Non-square FFT

Hello,

I wonder if one can perform a non-square DFT efficiently. Specifically, let's say I've got a grid function with N = 2 * K points, and I want to find the lower half (by frequency magnitude) of its DFT only. An obvious way would be to perform an FFT of the input grid function and then throw away the K higher frequency components. I wonder if there is a faster way to obtain the same result.

Thanks,

Dmitry

[Walter Roberson]

y = nufft(x,t,f) maybe ?? I have no idea what the performace is of nufft()

[Dmitry Medvedev]

Thanks! I'll try to check what exactly nufft does.

Answer 1

[Matt J]

There's probably no significant efficiency to be gained if you just want to get rid of half of the frequencies, as in the example you've mentioned. For a highly sparse FFT, for example if you want just want to keep the first 3 out of 1000 frequencies, it makes sense to pre-compute a non-square DFT matrix,

n=1000;

M=exp(-1j *(2*pi/n)* (0:n-1).*(0:2).'); % 3xn DFT matrix

and multiply it with your signal.

[Dmitry Medvedev]

Thanks! Yes, I have figured out that for a very small number of frequencies a direct summation would work. So far it seems that there is no efficient algorithm in a more general case. I hoped that there might be an algorithm precisely for the case that I described, where exactly half of the frequencies is of interest, since this is the case I am most likely to deal with directly. But for now it looks that there is no such algorithm.

Thanks,

Dmitry

[Matt J]

Here's a timing comparison. Basically, you need to be discarding the vast majority of the frequencies before FFT-and-discard becomes sub-optimal.

n=1e6;
x=rand(1,n);
t=0:numel(x)-1;
f=0:80;

M=exp(-1j *(2*pi/n)* (0:n-1).*f.'); % 3xn DFT matrix


tic; y=fft(x);y=y(f+1); toc  %Use FFT and discard

Elapsed time is 0.042146 seconds.

tic; M*x(:);toc              %Use pre-computed matrix

Elapsed time is 0.101585 seconds.

tic; nufft(x,t,f); toc       %Use NUFFT

Elapsed time is 0.206988 seconds.

[Dmitry Medvedev]

Thanks for the time comparison! Do I understand correctly that you keep only 81 out of the million of frequencies?

Thanks,

Dmitry

[Matt J]

Yes, the first 81.

[Walter Roberson]

This does not surprise me. fft() calls into high performance libraries that are coded in C or Fortran

[Chris Turnes]

Another option is to use the Chirp-Z Transform to compute a partial FFT, though for simple cases FFT + throw away is still probably more efficient.

However, if you have a simple division like wanting only the first half of the frequency vector, you could always apply the same recursion rules that the FFT itself does and just do two smaller FFTs, summing the result.

For instance, for a length-100 signal, if you only want the first 50 coefficients of the FFT you could do

x = randn(100,1);
% FFT + throw away.
y1 = fft(x);
y1 = y1(1:50);
% Manually compute two FFTs.
y2 = fft(x(1:2:end)) + exp(-1j*2*pi/100*(0:49)') .* fft(x(2:2:end));
norm(y1-y2)

ans = 1.9119e-14

[Dmitry Medvedev]

Thanks, Chris! Yes, I have come up with this algorithm, too. I will look into chirp-z transform, too.