基-4FFT和分裂基FFT的MATLAB仿真實現(xiàn),與Walsh變換比較

1、Harbin Institute of Technology數(shù)字信號處理作業(yè)院 系： 電信學(xué)院通信工程系姓 名：學(xué) 號：哈爾濱工業(yè)大學(xué)作業(yè)一：在MATLAB 上完成4096點的基-4FFT和分裂基 FFT的運算程序一、4096 點的基-4 FFT1、假設(shè)被分析信號為一個四頻信號x sin t2sinsin對該信號在時間t上每隔T 1s采樣一次，共取 4096個采樣點，得到長度為N 4096的序列x n ，對離散序列做基-4 FFT。2、基-4 FFT的MATLAB仿真結(jié)果在對基-4 FFT的結(jié)果進(jìn)行仿真時，為了驗證其正確性，使用了MATLAB中自帶的fft函數(shù)對分析信號進(jìn)行基-2 FFT變換，

2、與基-4 FFT變換的結(jié)果進(jìn)行對比。420-2-405001000150020002500300035004000分析信號的時域波形3000200010000matlab中基-2 FFT下 的頻譜300005001000150020002500300035004000基-4 FFT下的頻譜2000100005001000150020002500300035004000圖1-14096點的基-4 FFT的頻譜圖圖1-1是MATLAB中基-4 FFT的仿真結(jié)果，從圖中明顯可以看出，四頻信號 x的基-4 FFT計算出的頻譜圖與基-2 FFT計算出的頻譜完全相同，證明了基-4 FFT程序正確。3、基-

3、4 FFT 的 MATLAB 程序%DFT 的點數(shù)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 的分解級數(shù)， M=log4(N)Wn=exp(-2j*pi/N);%旋轉(zhuǎn)因子temp=zeros(1,N);%中間臨時數(shù)組，用于存儲各級蝶形運算計算結(jié)果%四進(jì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;%運算級之間的循環(huán) %第 L 級數(shù)據(jù)分組數(shù) %第 L-1 級數(shù)據(jù)分組數(shù)%第 L 級組間數(shù)據(jù)間隔個數(shù)， 也是組內(nèi)%第 L-1 級組間數(shù)據(jù)間隔個數(shù)，也是組%分組上限 %組內(nèi)數(shù)據(jù)上限for L=1:M group_cont_2=4A(M-L); group_cont_1=4A(M-L+1);group_interval_2=4AL;數(shù)據(jù)個數(shù)group_interval_1=4A(L-1);內(nèi)數(shù)據(jù)個數(shù)G=gro

6、up_cont_2-1;K=group_interval_1 -1 ;for g=0:G計算各組中的數(shù)據(jù)for k0=0:K%每一級中包含的組循環(huán)，遍歷每一組,%遍歷每一組中的所有數(shù)據(jù)，計算次級數(shù)據(jù)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;%每組數(shù)據(jù)中第一個數(shù)據(jù)序號%每一級所乘旋轉(zhuǎn)因子的指數(shù)因子 %每組數(shù)據(jù)中四個數(shù)據(jù)的序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;%遞推公式計算結(jié)果%將temp中臨時存儲的第L級計算結(jié)果賦值給x%作為次級運算的輸入end%繪制subplot(3,1,3); plot(abs(x);axis(0 4096 0 3000); title(基-4 FFT下的頻譜);4 FFT下的頻譜二、4096點的分裂基FFT1、假設(shè)被分析信號與前面相同，也是

8、四頻信號，234x sin t sin t sin t sin t2345對該信號每隔1秒取4096個采樣點，得到長度為N 4096的序列x n ，對離 散序列做分裂基FFT。2、分裂基FFT的MATLAB仿真結(jié)果在對分裂基FFT的結(jié)果進(jìn)行仿真時，為了驗證其正確性，使用 MATLAB中自帶的fft函數(shù)對分析信號進(jìn)行基-2 FFT變換，與分裂FFT變換的結(jié)果進(jìn)行對比分析信號的時域波形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的仿真結(jié)果，從圖中明顯可以看出，四頻 信號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結(jié)果的正確性subplot(3,1,2);plot(abs(y);axis(0 4096 0 3000); title(matlab 中基-2 FFT 下的頻譜);ticN=4096;% FFT 點數(shù)M=log(N)/log(4);% FFT 變換級數(shù)% 倒序運算 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% 進(jìn)行第一級的 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%進(jìn)行后面M-1級倒L形分裂基蝶形運算N2=2;for K=2:MN2=N2*2;% 每一級 DFT 點數(shù)N4=N2/4;% 每一級對應(yīng)的蝶形數(shù)目，一個基本碟形運算為倒 L 形for J=1:N4W_J_N2=exp(-j*2*pi*(J -1)/N2);% 旋轉(zhuǎn)因子；W_3J_N2=exp(-j*2*pi

12、*3*(J -1)/N2);% 旋轉(zhuǎn)因子；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作業(yè)二：比較基-2 FFT，基-4 FFT ,分裂基 FFT和 Walsh變換的運算時間FFT的復(fù)雜度主要來源于計算過程中與旋轉(zhuǎn)因子的乘法，基-2 FFT的每個基本蝶形中需乘1次旋轉(zhuǎn)因子，運算量為N log2 N ;基-4 FFT的每個基本蝶形中3N只需3次乘旋轉(zhuǎn)因子，運算量為 一 log4 N，分裂基FFT充分結(jié)合了基-2 FFT和4基-4 FFT的優(yōu)點，對

14、序列的偶序號項和奇數(shù)項項分別使用基 -2 FFT和基-4 FFT， 按照理論分析，分裂基FFT的運算量最小，基-4 FFT次之，基-2 FFT運算量最 大。變換基底后有得到 Walsh變換，可以由時域轉(zhuǎn)化到列率域來分析。為了比較幾種FFT和Walsh變換的運算量，統(tǒng)計幾種變換的運算時間并對比。 使用MATLAB中的tic/toc命令，tic用在程序的開始，作用是啟動一個計時器， 然后在程序尾部放一個toc，表示終止計時器，并返回tic啟動以來的總時間，單 位為秒(s)。1、FFT的運算時間由于軟件函數(shù)庫中自帶的基-2 fft函數(shù)和Walsh變換的fwht函數(shù)是預(yù)先編譯好 的內(nèi)部代碼，執(zhí)行效率非

15、常高，所以就會出現(xiàn)基-2 FFT的運算時間遠(yuǎn)少于基-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%進(jìn)行參數(shù)設(shè)置以及初始化N=4096;m=log2(N);n=0:N-1;nT=bin2dec(fliplr(dec2bin(n,m)+1;y=x(nT);%進(jìn)行基 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 運算結(jié)果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 程序計算出來，運行 一次程序得到以下結(jié)果：Elapsed time is 0.022429 seconds.2、Walsh變換的運算時間Walsh 變換的 MATLAB 程序如下：%Walsh變換快速算法的MATLAB程序clo

18、se all;clear all;clc;%生成四頻分析信號，點數(shù)為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變換結(jié)果圖subplot(2,1,2);plot(abs(y1);axis(0 4096 0 0.3);title(列率域 Walsh變換的結(jié)果);toc05001000150020002500300035004000圖1-3 Walsh變換的MATLAB仿真結(jié)果3、運算時間比較綜上所述，幾種變換的運算時間如下表所示變換類型（N=4096）運算時間基-2 FFTElapsed time is 0.129479 seco nds.