近似核FFT的实现及其仿真
参考《基于单比特和近似核FFT的瞬时测频算法及FPGA实现_袁俊榆_2.2_3.1》
参考《一种基于高阶近似核DFT的快速实现算法_吕远》
参考《数字信号处理_胡广书_4.2.2》
近日在查单比特测频方面的论文时，看到了使用近似核FFT来简化DFT运算的方式，颇感兴趣，研究一下。
近似核FFT
利用圆的外切正方形上的点，近似相同角度的Wn系数值，
对应到64点FFT的话，则系数如下图
代码实现
参照DIF-FFT的实现方式，函数fft_pipe使用多级蝶形运算级联的方式进行运算。值得注意的是，compensation的使用，当需要对Wn系数进行放缩的时候（诸如使用硬件实现），为了保证在运算过程中，各级联部分的上下两路数据处于同一量级，需要在不含Wn的运算式中对运算结果进行放大。
function [result] = fft_pipe(Xn,Wn,compensation)
pipe1 = zeros(1,64);
pipe2 = zeros(1,64);
pipe3 = zeros(1,64);
pipe4 = zeros(1,64);
pipe5 = zeros(1,64);
pipe6 = zeros(1,64);
%-- 第一级蝶形运算
% disp('----第一级蝶形运算');
for n=0:2:1
for m=1:32
pipe1(1,m+(n+0)*32) = compensation*(Xn(1,m+(n+0)*32) + Xn(1,m+(n+1)*32));
% index_0 = m+(n+0)*32-1;
% index_1 = m+(n+0)*32-1;
% index_2 = m+(n+1)*32-1;
% disp(['Y',num2str(index_0),'=(X',num2str(index_1),'+X',num2str(index_2),')']);
end
for m=1:32
pipe1(1,m+(n+1)*32) = (Xn(1,m+(n+0)*32) - Xn(1,m+(n+1)*32)) * Wn(1,(m-1)*1+1);
% index_0 = m+(n+1)*32-1;
% index_1 = m+(n+0)*32-1;
% index_2 = m+(n+1)*32-1;
% index_3 = (m-1)*1+1-1;
% disp(['Y',num2str(index_0),'=(X',num2str(index_1),'-X',num2str(index_2),') *Wn',num2str(index_3)]);
end
end
%-- 第二级蝶形运算
for n=0:2:3
for m=1:16
pipe2(1,m+(n+0)*16) = compensation*(pipe1(1,m+(n+0)*16) + pipe1(1,m+(n+1)*16));
end
for m=1:16
pipe2(1,m+(n+1)*16) = (pipe1(1,m+(n+0)*16) - pipe1(1,m+(n+1)*16)) * Wn(1,(m-1)*2+1);
end
end
%-- 第三级蝶形运算
for n=0:2:7
for m=1:8
pipe3(1,m+(n+0)*8) = compensation*(pipe2(1,m+(n+0)*8) + pipe2(1,m+(n+1)*8));
end
for m=1:8
pipe3(1,m+(n+1)*8) = (pipe2(1,m+(n+0)*8) - pipe2(1,m+(n+1)*8)) * Wn(1,(m-1)*4+1);
end
end
%-- 第四级蝶形运算
for n=0:2:15
for m=1:4
pipe4(1,m+(n+0)*4) = compensation*(pipe3(1,m+(n+0)*4) + pipe3(1,m+(n+1)*4));
end
for m=1:4
pipe4(1,m+(n+1)*4) = (pipe3(1,m+(n+0)*4) - pipe3(1,m+(n+1)*4)) * Wn(1,(m-1)*8+1);
end
end
%-- 第五级蝶形运算
for n=0:2:31
for m=1:2
pipe5(1,m+(n+0)*2) = compensation*(pipe4(1,m+(n+0)*2) + pipe4(1,m+(n+1)*2));
end
for m=1:2
pipe5(1,m+(n+1)*2) = (pipe4(1,m+(n+0)*2) - pipe4(1,m+(n+1)*2)) * Wn(1,(m-1)*16+1);
end
end
%-- 第六级蝶形运算
for n=0:2:63
for m=1:1
pipe6(1,m+(n+0)*1) = compensation*(pipe5(1,m+(n+0)*1) + pipe5(1,m+(n+1)*1));
end
for m=1:1
pipe6(1,m+(n+1)*1) = (pipe5(1,m+(n+0)*1) - pipe5(1,m+(n+1)*1)) * Wn(1,(m-1)*32+1);
end
end
%-- 逆序取数
result = zeros(1,64);
for m = 0:63
index = bin2dec(reverse(dec2bin(m,6)))+1;
result(1,index) = pipe6(1,m+1);
end
end
%清空一切
clc;clear all;close all;
%获得 近似核系数
W00 = complex(+4.0,-0.0); W08 = complex(4.0,-4); W16 = complex(-0.0,-4); W24 = complex(-4,-4.0);
W01 = complex(+4.0,-0.5); W09 = complex(3.5,-4); W17 = complex(-0.5,-4); W25 = complex(-4,-3.5);
W02 = complex(+4.0,-1.0); W10 = complex(3.0,-4); W18 = complex(-1.0,-4); W26 = complex(-4,-3.0);
W03 = complex(+4.0,-1.5); W11 = complex(2.5,-4); W19 = complex(-1.5,-4); W27 = complex(-4,-2.5);
W04 = complex(+4.0,-2.0); W12 = complex(2.0,-4); W20 = complex(-2.0,-4); W28 = complex(-4,-2.0);
W05 = complex(+4.0,-2.5); W13 = complex(1.5,-4); W21 = complex(-2.5,-4); W29 = complex(-4,-1.5);
W06 = complex(+4.0,-3.0); W14 = complex(1.0,-4); W22 = complex(-3.0,-4); W30 = complex(-4,-1.0);
W07 = complex(+4.0,-3.5); W15 = complex(0.5,-4); W23 = complex(-3.5,-4); W31 = complex(-4,-0.5);
W32 = complex(-4.0,0.0); W40 = complex(-4.0,4); W48 = complex(0.0,4); W56 = complex(4,4.0);
W33 = complex(-4.0,0.5); W41 = complex(-3.5,4); W49 = complex(0.5,4); W57 = complex(4,3.5);
W34 = complex(-4.0,1.0); W42 = complex(-3.0,4); W50 = complex(1.0,4); W58 = complex(4,3.0);
W35 = complex(-4.0,1.5); W43 = complex(-2.5,4); W51 = complex(1.5,4); W59 = complex(4,2.5);
W36 = complex(-4.0,2.0); W44 = complex(-2.0,4); W52 = complex(2.0,4); W60 = complex(4,2.0);
W37 = complex(-4.0,2.5); W45 = complex(-1.5,4); W53 = complex(2.5,4); W61 = complex(4,1.5);
W38 = complex(-4.0,3.0); W46 = complex(-1.0,4); W54 = complex(3.0,4); W62 = complex(4,1.0);
W39 = complex(-4.0,3.5); W47 = complex(-0.5,4); W55 = complex(3.5,4); W63 = complex(4,0.5);
W_approx = zeros(1,64);
W_approx(1,0*8+1:1*8) = [W00 W01 W02 W03 W04 W05 W06 W07];
W_approx(1,1*8+1:2*8) = [W08 W09 W10 W11 W12 W13 W14 W15];
W_approx(1,2*8+1:3*8) = [W16 W17 W18 W19 W20 W21 W22 W23];
W_approx(1,3*8+1:4*8) = [W24 W25 W26 W27 W28 W29 W30 W31];
W_approx(1,4*8+1:5*8) = [W32 W33 W34 W35 W36 W37 W38 W39];
W_approx(1,5*8+1:6*8) = [W40 W41 W42 W43 W44 W45 W46 W47];
W_approx(1,6*8+1:7*8) = [W48 W49 W50 W51 W52 W53 W54 W55];
W_approx(1,7*8+1:8*8) = [W56 W57 W58 W59 W60 W61 W62 W63];
% 获得精确的Wn系数
for m=0:63
W_exact(1,m+1) = exp(-j*2*pi*m/64);
end
%生成单频信号
N = 64;
t = (0:N-1);
sig = sin(2*pi*0.2*t);
Xn = 25*sig;
%近似核算法运算
Wn = W_approx*2; %W_approx约是W_exact的4倍
compensation = 8;
[result] = fft_pipe(Xn,Wn,compensation);
figure;
plot(20*log10(abs(result)/max(abs(result))));
title(['近似核多级蝴蝶运算结果']);
%精确FFT运算
Wn = floor(1024*W_exact);
Xn = 25*sig;
compensation = 1024;
[result] = fft_pipe(Xn,Wn,compensation);
figure;
plot(20*log10(abs(result)/max(abs(result))));
title(['精确的多级蝴蝶运算结果']);
%内置FFT函数运算
figure;
plot(20*log10(abs(fft(Xn,64))/max(abs(fft(Xn,64)))));
title('内置FFT运算结果');
W_approx是近似核的系数，W_exact是精确的Wn系数，仿真过程很简单，产生单频信号，改变幅度、W系数的幅度。使用近似核FFT的话，动态会差一些。