如何使用较小的FFT计算大尺寸FFT？

Question

如何使用较小的FFT计算大尺寸FFT？

3

如果我有一个大小为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）专家，因此如果能用简单的术语来解释，那将不胜感激。

- ysap

1

只需执行一次基数为4的Cooley-Tukey FFT算法。这将把FFT分解成四个大小为1/4的FFT。 - Mysticial

4个回答

1

kevin_o的回答效果很好。我采用了他的代码，并使用一些基本的Matlab技巧消除了循环。它在功能上与他的版本完全相同。

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()

- Ed_B

这样写更易读了。但是第五步仍然过于晦涩难懂。其实很简单：y = reshape(X.', [], 1);，就这样。这个函数不过是标准的频率抽取FFT方法而已。在YouTube上有数百个相关视频。 - Harry

0

嗯，FFT基本上是一种递归类型的傅里叶变换。它依赖于维基百科所述的事实：

最著名的FFT算法依赖于N的因式分解，但是有O（N log N）复杂度的FFT适用于所有N，即使是素数N。许多FFT算法仅依赖于e^（-2pi * i / N）是第N个本原单位根这一事实，因此可以应用于任何有限域上的类似变换，例如数论变换。由于反向DFT与DFT相同，但指数中的符号相反并且有1 / N因子，因此任何FFT算法都可以轻松地为其进行调整。

因此，在FFT中已经基本完成了这个过程。如果您想从转换中获取更长时间的信号，则最好在有限频率的数据集上执行DFT。也许有一种方法可以从频域中完成它，但我不知道是否有人真正做到了。你可以成为第一个尝试的人!!!! :)

- Kevin R.

0

你可以只使用基数为2的FFT的最后log2(k)次传递，前提是前面的FFT结果来自适当交错的数据子集。

- hotpaw2

谢谢。这个假设是什么意思？我假设你的意图是在执行FFT_M()步骤后，重新组织这些变换的输出，然后像往常一样继续应用传递。 - ysap

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- kevin_o · Accepted Answer

这是一种计算大小为P的FFT的算法，使用两个较小的FFT函数，它们的大小分别为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的每一列。
步骤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().