How do I compute a large FFT using a smaller FFT?
If I have an FFT implementation of a certain size M (cardinality 2), how can I calculate the FFT of a set of sizes P = k * M, where k is also a power of 2?
#define M 256
#define P 1024
complex float x[P];
complex float X[P];
// Use FFT_M(y) to calculate X = FFT_P(x) here
[The question is fully expressed in a general sense. I know that FFT calculation is a huge field and many architectural optimizations have been researched and developed, but I'm trying to figure out how this can be done at a more abstract level. Please note that I am not an FFT expert (or DFT for that matter), so if the explanation can be phrased in simple terms that will be appreciated]
source to share
Here is an algorithm for computing a P-size FFT using two smaller FFT functions, sizes M and N (the original question begs the sizes M and k).
Inputs:
P is the size of the large FFT you want to compute.
M, N are chosen such that MN = P.
x [0 ... P-1] - input data.
Setting:
U is a two-dimensional array with M rows and N columns.
y is a vector of length P that will contain the FFT x.
Algorithm:
Step 1. Fill U from x with columns so that U looks like this: x(0) x(M) ... x(P-M)
x(1) x(M+1) ... x(P-M+1)
x(2) x(M+2) ... x(P-M+2)
... ... ... ...
x(M-1) x(2M-1) ... x(P-1)
Step 2. Replace each row of U with your own FFT (length N).
Step 3. Multiply each element U (m, n) by exp (-2 * pi * j * m * n / P).
Step 4. Replace each U column with your own FFT (length M).
Step 5. Consider the elements of U as strings in y, for example: y(0) y(1) ... y(N-1)
y(N) y(N+1) ... y(2N-1)
y(2N) y(2N+1) ... y(3N-1)
... ... ... ...
y(P-N) y(P-N-1) ... y(P-1)
Here is the MATLAB code that implements this algorithm. You can check it by typingfft_decomposition(randn(256,1), 8);
function y = fft_decomposition(x, M)
% y = fft_decomposition(x, M)
% Computes FFT by decomposing into smaller FFTs.
%
% Inputs:
% x is a 1D array of the input data.
% M is the size of one of the FFTs to use.
%
% Outputs:
% y is the FFT of x. It has been computed using FFTs of size M and
% length(x)/M.
%
% Note that this implementation doesn't explicitly use the 2D array U; it
% works on samples of x in-place.
q = 1; % Offset because MATLAB starts at one. Set to 0 for C code.
x_original = x;
P = length(x);
if mod(P,M)~=0, error('Invalid block size.'); end;
N = P/M;
% step 2: FFT-N on rows of U.
for m = 0 : M-1
x(q+(m:M:P-1)) = fft(x(q+(m:M:P-1)));
end;
% step 3: Twiddle factors.
for m = 0 : M-1
for n = 0 : N-1
x(m+n*M+q) = x(m+n*M+q) * exp(-2*pi*j*m*n/P);
end;
end;
% step 4: FFT-M on columns of U.
for n = 0 : N-1
x(q+n*M+(0:M-1)) = fft(x(q+n*M+(0:M-1)));
end;
% step 5: Re-arrange samples for output.
y = zeros(size(x));
for m = 0 : M-1
for n = 0 : N-1
y(m*N+n+q) = x(m+n*M+q);
end;
end;
err = max(abs(y-fft(x_original)));
fprintf( 1, 'The largest error amplitude is %g\n', err);
return;
% End of fft_decomposition().
source to share
Well, FFT is basically a recursive type of Fourier transform. He relies on the fact that, as wikipedia writes,
The most famous FFT algorithms depend on the factorization of N, but there are FFTs with complexity O (N log N) for> all N, even for prime N. Many FFT algorithms depend only on the fact that e ^ (-2pi * i / N) is The Nth primitive root of unity, and> can thus be applied to similar transformations over any finite field, such as number-theoretic transformations. Since the> inverse DFT is the same as the DFT, but with the opposite sign in the exponent and a factor of 1 / N, any FFT algorithm> can be easily adapted for it.
So this is already done in FFT. If you're talking about getting longer signals from your transform, you're better off doing DFT on frequency-limited datasets. There might be a way to do this from the frequency domain, but IDK if anyone actually did it. You could be the first !!!! :)
source to share
Kevin_o's answer worked pretty well. I took his code and eliminated the loops using some basic Matlab tricks. It is functionally identical to its version.
function y = fft_decomposition(x, M)
% y = fft_decomposition(x, M)
% Computes FFT by decomposing into smaller FFTs.
%
% Inputs:
% x is a 1D array of the input data.
% M is the size of one of the FFTs to use.
%
% Outputs:
% y is the FFT of x. It has been computed using FFTs of size M and
% length(x)/M.
%
% Note that this implementation doesn't explicitly use the 2D array U; it
% works on samples of x in-place.
q = 1; % Offset because MATLAB starts at one. Set to 0 for C code.
x_original = x;
P = length(x);
if mod(P,M)~=0, error('Invalid block size.'); end;
N = P/M;
% step 2: FFT-N on rows of U.
X=fft(reshape(x,M,N),[],2);
% step 3: Twiddle factors.
X=X.*exp(-j*2*pi*(0:M-1)'*(0:N-1)/P);
% step 4: FFT-M on columns of U.
X=fft(X);
% step 5: Re-arrange samples for output.
x_twiddle=bsxfun(@plus,M*(0:N-1)',(0:M-1))+q;
y=X(x_twiddle(:));
% err = max(abs(y-fft(x_original)));
% fprintf( 1, 'The largest error amplitude is %g\n', err);
return;
% End of fft_decomposition()
source to share