级联FFT（超采样FFT架构）的MATLAB代码及原理

一、原理

  该算法首先对输入时间序列的数据进行抽样，然后对抽样后数组内的数据进行 FFT 运算处理，然后进行交叉项的补偿，再对 FFT 之后不同数组间相同位置上的数据进行第 2 次 FFT 处理，从而达到一次 FFT 运算能够得到的效果。
  例如下图，做N点FFT，可以先将输入折叠成N=L*M矩阵形式，然后先对每行进行FFT（长度为M），得到的结果乘以补偿因子后，再对每列进行FFT（长度为L）。这样，便通过两级短FFT实现了一级长FFT同样的效果。

二、适用场景

  此方法适用于FFT长度较长而使得硬件的DSP资源无法满足的场景，其他情况下使用此方法并不能明显节约硬件资源。

三、MATLAB代码

  SuperSamplingFFT.m

function XN1N2 = SuperSamplingFFT(Xin,SubLength)

if(~isvector(Xin))
    error('Xin should be a vector');
end

row = isrow(Xin);

N = length(Xin);
N1 = SubLength;
N2 = ceil(N/N1); % Length of each sub-series

if(N1*N2~=N)
    error('SubLength should be a divider of the length of Xin');
end

i = sqrt(-1);
WN1 = exp(-2*i*pi/N1);
WN2 = exp(-2*i*pi/N2);
WN1N2 = exp(-2*i*pi/(N1*N2));


xN1N2 = reshape(Xin,N1,N2);

% FFT by row
YN1N2 = zeros(N1,N2);
for n1=1:N1
    YN1N2(n1,:) = fft(xN1N2(n1,:));
end

% Weighting 
YYN1N2 = zeros(N1,N2);
for n2=1:N2
    Y = YN1N2(:,n2);  % vector extraction
    Weight = WN1N2.^([0:N1-1]'*(n2-1));
    YN1N2(:,n2) = Y.*Weight;
end

% Second FFT
for n2=1:N2
    YYN1N2(:,n2) = fft(YN1N2(:,n2));
end

% Transpose and back to vector
if(row)
    XN1N2 = reshape(YYN1N2.',1,N);
else
    XN1N2 = reshape(YYN1N2.',N,1);
end

end

  Test_SuperSamplingFFT.m

N = 2048;
N1 = 512; % Number of sub-series to extract from original time-series

xN = rand(N,1)+i*rand(N,1);
XN = fft(xN);
XN1N2 = SuperSamplingFFT(xN,N1);

figure(1);
clf;
subplot(2,2,1);
plot([1:N],real(XN),'b',[1:N],real(XN1N2),'r');
title(['Real Part:  max error ' num2str(max(abs(real(XN-XN1N2))))]);

subplot(2,2,2);
plot([1:N],imag(XN),'b',[1:N],imag(XN1N2),'r');
title(['Imag Part:  max error ' num2str(max(abs(imag(XN-XN1N2))))]);

disp(['Maximum Modulus Error: ' num2str(max(abs(XN-XN1N2)))]);

xN = rand(1,N)+i*rand(1,N);
XN = fft(xN);
XN1N2 = SuperSamplingFFT(xN,N1);

subplot(2,2,3);
plot([1:N],real(XN),'b',[1:N],real(XN1N2),'r');
title(['Real Part:  max error ' num2str(max(abs(real(XN-XN1N2))))]);

subplot(2,2,4);
plot([1:N],imag(XN),'b',[1:N],imag(XN1N2),'r');
title(['Imag Part:  max error ' num2str(max(abs(imag(XN-XN1N2))))]);

disp(['Maximum Modulus Error: ' num2str(max(abs(XN-XN1N2)))]);

四、参考来源

  参考一：一种新的级联 FFT 算法, 张大炜, (中国电子进出口总公司，北京 100037)
  参考二：https://www.dsprelated.com/thread/3440/cascaded-ffts（这是Xilinx写的实现方法）
  下载链接：参考二中可以下载Xilinx写的实现文档，就是一个pdf文件，下载不了的话可以在我CSDN中下载，https://download.csdn.net/download/qq_35809085/90533630