Search code examples
fftlarge-datadft

MATLAB Dftmtx for Huge "N"


I have a vector of size M (say 500), which I up-sample by a factor of MM=500, so that my new vector is now size N=500 x 500=250000. I am using an optimisation algorithm, and need to carryout the fft/dft of the up-sampled vector of size N using the DFT Matrix, and not the inbuilt function.

However, this becomes prohibitive due to memory constraints. Is there any way to go about it? I have seen a similar question here Huge Fourier matrix - MATLAB but this is regarding a Huge Matrix, where the solution is to break the matrix into columns and do the operation column by column. In my case, the vector has 250000 rows.

Would it be wise to split the rows into pieces, say 500 each and iterate the same thing 500 times, and concatenate the results in end ?


Solution

  • If using the FFT is an option, the matrix of twiddle factors does not appear explicitly, so the actual memory requirements are on the order of O(N).

    If you must use the explicit DFT matrix, then it is possible to break down the computations using submatrices of the larger DFT matrix. Given an input x of length N, and assuming we wish to divide the large DFT matrix into BlockSize x BlockSize submatrices, this can be done with the following matlab code:

    y = zeros(size(x));
    Imax = ceil(N / BlockSize); % divide the rows into Imax chunks
    Jmax = ceil(N / BlockSize); % divide the columns into Jmax chunks
    
    % iterate over the blocks
    for i=0:Imax-1
      imin = i*BlockSize;
      imax = min(i*BlockSize+BlockSize-1,N-1);
      for j=0:Jmax-1
        jmin = j*BlockSize;
        jmax = min(j*BlockSize+BlockSize-1,N-1);
        [XX,YY] = meshgrid(jmin:jmax, imin:imax);
    
        % compute the DFT submatrix
        W = exp(-2* pi * 1i * XX .* YY / N);
    
        % apply the DFT submatrix on a chunk of the input and add to the output
        y([imin:imax] + 1) = y([imin:imax] + 1) + W * x([jmin:jmax] + 1);
      end
    end
    

    If needed it would be fairly easy to adapt the above code to use different block size along the rows than along the columns.