如何使用较小尺寸的FFT计算大尺寸FFT？
如果我有一个特定大小为M(2的幂)的快速傅立叶变换实现，我如何计算一组大小为P=k*M的快速傅立叶变换，其中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，就这一点而言)专家，所以如果可以用简单的术语来解释，将不胜感激
回答 4
Stack Overflow用户
发布于 2014-03-29 07:12:20
这里有一个算法，使用两个大小为M和N的较小的FFT函数来计算大小为P的FFT (原始问题称为大小为M和k)。
输入：
P是要计算的大型FFT的大小。
M、N被选择成使得MN=P。
x0...P-1是输入数据。
设置：
U是具有M行和N列的2D阵列。
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.将U的每一行替换为它自己的FFT (长度为N)。
步骤3.将U(m，n)的每个元素乘以exp(-2*pi*j*m*n/P)。
步骤4.将U的每一列替换为它自己的FFT (长度为M)。
步骤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().Stack Overflow用户
发布于 2017-09-12 10:00:42
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()Stack Overflow用户
发布于 2012-03-30 20:45:22
您可以只使用基数-2FFT的最后log2(k)遍，假设之前的FFT结果来自适当交错的数据子集。
https://stackoverflow.com/questions/9940192
复制相似问题
腾讯云开发者
