fft - 如何使用较小的 FFT 计算较大的 FFT？

Question

如果我有一个特定大小 M(2 的幂)的 FFT 实现，我如何计算一组大小 P=k*M 的 FFT，其中 k 也是 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

[这个问题是故意在一般意义上表达的。我知道 FFT 计算是一个广阔的领域，研究和开发了许多特定于体系结构的优化，但我想了解的是这在更抽象的层面上是如何实现的。请注意，我不是 FFT(或 DFT，就此而言)专家，因此，如果可以用简单的术语进行解释，我们将不胜感激]

Answer 1

这是一个使用两个较小的 FFT 函数计算大小 P 的算法，大小为 M 和 N(原始问题称为大小 M 和 k)。

输入:
P 是您希望计算的大型 FFT 的大小。
选择 M, N 使得 MN=P。
x[0...P-1]为输入数据。

设置:
U 是一个 M 行 N 列的二维数组。
y 是长度为 P 的向量，它将保存 x 的 FFT。

算法:
步骤 1. 从 x 中按列填充 U，使 U 如下所示:
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)

步骤 2. 用自己的 FFT(长度为 N)替换 U 的每一行。
步骤 3. 将 U(m,n) 的每个元素乘以 exp(-2*pi*j*m*n/P)。
步骤 4. 用自己的 FFT(长度为 M)替换 U 的每一列。
step 5. 逐行读出U的元素到y中，像这样:

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)

这是实现该算法的 MATLAB 代码。您可以通过键入 fft_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().