基-4FFT和分裂基FFT的MATLAB仿真实现,与Walsh变换比较

1、Harbin Institute of Technology数字信号处理作业院 系： 电信学院通信工程系姓 名：学 号：哈尔滨工业大学作业一：在MATLAB 上完成4096点的基-4FFT和分裂基 FFT的运算程序一、4096 点的基-4 FFT1、假设被分析信号为一个四频信号x sin t2sinsin对该信号在时间t上每隔T 1s采样一次，共取 4096个采样点，得到长度为N 4096的序列x n ，对离散序列做基-4 FFT。2、基-4 FFT的MATLAB仿真结果在对基-4 FFT的结果进行仿真时，为了验证其正确性，使用了MATLAB中自带的fft函数对分析信号进行基-2 FFT变换，

2、与基-4 FFT变换的结果进行对比。420-2-405001000150020002500300035004000分析信号的时域波形3000200010000matlab中基-2 FFT下 的频谱300005001000150020002500300035004000基-4 FFT下的频谱2000100005001000150020002500300035004000图1-14096点的基-4 FFT的频谱图图1-1是MATLAB中基-4 FFT的仿真结果，从图中明显可以看出，四频信号 x的基-4 FFT计算出的频谱图与基-2 FFT计算出的频谱完全相同，证明了基-4 FFT程序正确。3、基-

3、4 FFT 的 MATLAB 程序%DFT 的点数clc; clear; N=4096;%生成四频的分析信号t=0:4095;x=sin(1*pi/2*t)+sin(2*pi/3*t)+sin(3*pi/4*t)+sin(4*pi/5*t);%四频的信号subplot(3,1,1);plot(x);axis(0 4096 -4 4);title( 分析信号的时域波形 );subplot(3,1,2);%绘制基-2 FFT的频谱图，用于比较plot(abs(fft(x);axis(0 4096 0 3000);title(matlab 中基-2 FFT 下的频谱);M=log(N)/log(4)

4、;%基 4FFT 的分解级数， M=log4(N)Wn=exp(-2j*pi/N);%旋转因子temp=zeros(1,N);%中间临时数组，用于存储各级蝶形运算计算结果%四进制逆序n=0:N- 1 ;screen=ones(1,N);n=bitor(bitand(n,screen*hex2dec(cccc)/4,bitand(n,screen*hex2dec(3333)*4);n=bitor(bitand(n,screen*hex2dec(f0f0)/16,bitand(n,screen*hex2dec(0f0f)*16);n=bitor(bitand(n,screen*hex2dec(ff

5、00)/256,bitand(n,screen*hex2dec(00ff)*256); n=n /4A(8-M)+1;for n1=1:Ntemp(n(n1)=x(n1);end x=temp;%运算级之间的循环 %第 L 级数据分组数 %第 L-1 级数据分组数%第 L 级组间数据间隔个数， 也是组内%第 L-1 级组间数据间隔个数，也是组%分组上限 %组内数据上限for L=1:M group_cont_2=4A(M-L); group_cont_1=4A(M-L+1);group_interval_2=4AL;数据个数group_interval_1=4A(L-1);内数据个数G=gro

6、up_cont_2-1;K=group_interval_1 -1 ;for g=0:G计算各组中的数据for k0=0:K%每一级中包含的组循环，遍历每一组,%遍历每一组中的所有数据，计算次级数据k=kO+g*group_i nterval_2+1; m=group_c on t_2*k0;k1=k;k2=k1+group_i nterval_1; k3=k2+group_i nterval_1;k4 =k3+group_i nterval_1;%每组数据中第一个数据序号%每一级所乘旋转因子的指数因子 %每组数据中四个数据的序end end x=temp;X1=x(k1);X2=WnAm*x

7、(k2);X3=WnA(2*m)*x(k3);X4=WnA(3*m)*x(k4);temp(k1)=X1+X2+X3+X4; temp(k2)=X1-1j*X2 -X3+1j*X4; temp(k3)=X1 -X2+X3 -X4; temp(k4)=X1+1j*X2 -X3-1j*X4;%递推公式计算结果%将temp中临时存储的第L级计算结果赋值给x%作为次级运算的输入end%绘制subplot(3,1,3); plot(abs(x);axis(0 4096 0 3000); title(基-4 FFT下的频谱);4 FFT下的频谱二、4096点的分裂基FFT1、假设被分析信号与前面相同，也是

8、四频信号，234x sin t sin t sin t sin t2345对该信号每隔1秒取4096个采样点，得到长度为N 4096的序列x n ，对离 散序列做分裂基FFT。2、分裂基FFT的MATLAB仿真结果在对分裂基FFT的结果进行仿真时，为了验证其正确性，使用 MATLAB中自带的fft函数对分析信号进行基-2 FFT变换，与分裂FFT变换的结果进行对比分析信号的时域波形matlab中基-2 FFT下 的频谱3000LILLIILL 12000 -1000 -0丄 IJ.05001000150020002500300035004000分裂基FFT下的频谱3000 1111cc1r-T

9、2000 -1000 -0 iIIi_.05001000150020002500300035004000图1-24096点的分裂基FFT的频谱图图1-2是MATLAB中分裂基FFT的仿真结果，从图中明显可以看出，四频 信号x的分裂基FFT计算出的频谱图与基-2 FFT计算出的频谱完全相同，证明 了分裂基FFT程序正确。3、分裂基FFT的MATLAB程序%时域抽取法分裂基 FFT算法clear;clc;t=0:(N-1);x=si n(1*pi/2*t+si n(2*pi/3*t+si n(3*pi/4*t+si n(4*pi/5*t;% 四频输入信号subplot(3,1,1);plot(x)

10、;axis(0 4096 -3 3),title(分析信号的时域波形);y=fft(x,4096);%绘制基-2 FFT的频谱图,用于验证基-4 FFT结果的正确性subplot(3,1,2);plot(abs(y);axis(0 4096 0 3000); title(matlab 中基-2 FFT 下的频谱);ticN=4096;% FFT 点数M=log(N)/log(4);% FFT 变换级数% 倒序运算 J=1;N1=N-1;for I=1:N1if I JT=x(I);x(I)=x(J);x(J)=T;endK=N/2;while K JJ=J-K;K=K/2;endJ=J+K;e

11、nd% 进行第一级的 2 点 FFT 运算 IS=1;step=4;while ISNfor n=IS:step:N %2 点 DFTa=x(n)+x(n+1);b=x(n)-x(n+1);x(n)=a;%原位计算x(n+1)=b;endIS=2*step-1;step=4*step;end%进行后面M-1级倒L形分裂基蝶形运算N2=2;for K=2:MN2=N2*2;% 每一级 DFT 点数N4=N2/4;% 每一级对应的蝶形数目，一个基本碟形运算为倒 L 形for J=1:N4W_J_N2=exp(-j*2*pi*(J -1)/N2);% 旋转因子；W_3J_N2=exp(-j*2*pi

12、*3*(J -1)/N2);% 旋转因子；IS=J;step=2*N2;while ISNfor n=IS:step:N-1 a=x(n)+x(n+2*N4)*W_J_N2+x(n+3*N4)*W_3J_N2; b=x(n+N4)-j*(x(n+2*N4)*W_J_N2 -x(n+3*N4)*W_3J_N2); c=x(n)-x(n+2*N4)*W_J_N2 -x(n+3*N4)*W_3J_N2; d=x(n+N4)+j*(x(n+2*N4)*W_J_N2 -x(n+3*N4)*W_3J_N2); x(n)=a;x(n+N4)=b; x(n+2*N4)=c; x(n+3*N4)=d;endIS

13、=2*step-N2+J; step=4*step;endendend% 绘制分裂基 FFT 频谱图 subplot(3,1,3);plot(abs(x);axis(0 4096 0 3000);title(分裂基FFT下的频谱);toc作业二：比较基-2 FFT，基-4 FFT ,分裂基 FFT和 Walsh变换的运算时间FFT的复杂度主要来源于计算过程中与旋转因子的乘法，基-2 FFT的每个基本蝶形中需乘1次旋转因子，运算量为N log2 N ;基-4 FFT的每个基本蝶形中3N只需3次乘旋转因子，运算量为 一 log4 N，分裂基FFT充分结合了基-2 FFT和4基-4 FFT的优点，对

14、序列的偶序号项和奇数项项分别使用基 -2 FFT和基-4 FFT， 按照理论分析，分裂基FFT的运算量最小，基-4 FFT次之，基-2 FFT运算量最 大。变换基底后有得到 Walsh变换，可以由时域转化到列率域来分析。为了比较几种FFT和Walsh变换的运算量，统计几种变换的运算时间并对比。 使用MATLAB中的tic/toc命令，tic用在程序的开始，作用是启动一个计时器， 然后在程序尾部放一个toc，表示终止计时器，并返回tic启动以来的总时间，单 位为秒(s)。1、FFT的运算时间由于软件函数库中自带的基-2 fft函数和Walsh变换的fwht函数是预先编译好 的内部代码，执行效率非

15、常高，所以就会出现基-2 FFT的运算时间远少于基-4FFT。例如，在做4096点基-4 FFT的仿真时，4096点基-2 FFT和基-4 FFT的运 算时间分别是：Elapsed time is 0.013441 seco nds.(基-2 FFT) Elapsed time is 0.032516 seco nds.(基-4 FFT)为了能够公平地比较几种变换的运算时间，在这里给出了一个自编的基-2FFT程序，如下所示：%基2FFT算法的MATLAB程序close all;clear all;clc;%生成分析信号t=0:4095;x=si n(1*pi/2*t+si n(2*pi/3*t

16、+si n(3*pi/4*t+si n(4*pi/5*t;% 四频输入信号subplot(2,1,1),plot(x);axis(0 4095 -2 2);title(分析信号的时域波形);tic%进行参数设置以及初始化N=4096;m=log2(N);n=0:N-1;nT=bin2dec(fliplr(dec2bin(n,m)+1;y=x(nT);%进行基 2FFT 变换运算for mm=1:m;n z=2Amm;u=1;wn=exp(-i*2*pi/nz);for j=1:nz/2;for k=j:nz:N;kp=k+nz/2;t=y(kp)*u;y(kp)=y(k)-t;y(k)=y(k

17、)+t;endu=u*wn;endend%绘制基 2FFT 运算结果subplot(2,1,2);plot(abs(y);axis(0 4096 0 3000);title(基2FFT计算出的频谱);toc由以上程序得到的基-2 FFT运算时间为：Elapsed time is 0.129479 seconds.分裂基 FFT 的运算时间可以直接由前面的分裂基 FFT 程序计算出来，运行 一次程序得到以下结果：Elapsed time is 0.022429 seconds.2、Walsh变换的运算时间Walsh 变换的 MATLAB 程序如下：%Walsh变换快速算法的MATLAB程序clo

18、se all;clear all;clc;%生成四频分析信号，点数为4096a=0:4095;x=sin(1*pi/2*a)+sin(2*pi/3*a)+sin(3*pi/4*a)+sin(4*pi/5*a); subplot(2,1,1);plot(x);axis(O 4095 -4 4);title(分析信号的时域波形);ticy1= fWht(x,hadamard);%绘制 Walsh变换结果图subplot(2,1,2);plot(abs(y1);axis(0 4096 0 0.3);title(列率域 Walsh变换的结果);toc05001000150020002500300035004000图1-3 Walsh变换的MATLAB仿真结果3、运算时间比较综上所述，几种变换的运算时间如下表所示变换类型（N=4096）运算时间基-2 FFTElapsed time is 0.129479 seco nds.