快速傅里葉變換基2時(shí)間抽取FFT算法

7.4實(shí)驗(yàn)4：快速傅里葉變換-基2時(shí)間抽取FFT算法matlab實(shí)現(xiàn)實(shí)驗(yàn)?zāi)康?.練習(xí)利用matlab6.5中工具箱中的信號(hào)處理函數(shù)2.熟悉快速傅里葉變換的根本原理3.熟悉基2DIT-FFT運(yùn)算的MATLAB程序并運(yùn)用涉及函數(shù)信號(hào)處理函數(shù)X=fft(x)或者X=fft(x,N):自定義功能函數(shù)functiony=myfft(xr,n)實(shí)驗(yàn)原理與方法1DIT-FFT算法的根本原理有限長序列x(n)的N點(diǎn)DFT定義為：，式中，其整數(shù)次冪簡(jiǎn)稱為旋轉(zhuǎn)因子。N符合2的整數(shù)冪，N為2的幾次冪，那么需要進(jìn)行幾次分解。碟形運(yùn)算流圖符號(hào)如下：2DIT-FFT算法的運(yùn)算規(guī)律及編程思想為了編寫DIT-FFT算法的運(yùn)算程序，首先要分析其運(yùn)算規(guī)律，總結(jié)編程思想并繪出程序框圖。由右圖可知，DIT-FFT算法的運(yùn)算過程很有規(guī)律。2.1原位計(jì)算對(duì)點(diǎn)的FFT共進(jìn)行M級(jí)運(yùn)算，每級(jí)由N/2個(gè)蝶形運(yùn)算組成。在同一級(jí)中，每個(gè)蝶的輸入數(shù)據(jù)只對(duì)本蝶有用，且輸出節(jié)點(diǎn)與輸入節(jié)點(diǎn)在同一水平線上，這就意味著每算完一個(gè)蝶后，所得數(shù)據(jù)可立即存入原輸入數(shù)據(jù)所占用的數(shù)組元素(存儲(chǔ)單元)，這種原位(址)計(jì)算的方法可節(jié)省大量內(nèi)存。2.2蝶形運(yùn)算實(shí)現(xiàn)FFT運(yùn)算的核心是蝶形運(yùn)算，找出蝶形運(yùn)算的規(guī)律是編程的根底。formm=1:m%將DFT做m次基2分解，從左到右，對(duì)每次分解作DFT運(yùn)算Nmr=2^mm;u=1;%旋轉(zhuǎn)因子u初始化WN=exp(-j*2*pi/Nmr);%本次分解的根本DFT因子WN＝exp(-i*2*pi/Nmr)forn=1:Nmr/2%本次跨越間隔內(nèi)的各次碟形運(yùn)算fork=n:Nmr:N%本次碟形運(yùn)算的跨越間隔為Nmr=2^mmkp=k+Nmr/2;%確定碟形運(yùn)算的對(duì)應(yīng)單元下標(biāo)(對(duì)稱性)t=x(kp)*u;%碟形運(yùn)算的乘積項(xiàng)x(kp)=x(k)-t;%碟形運(yùn)算的加法項(xiàng)x(k)=x(k)+t;endu=u*WN;%修改旋轉(zhuǎn)因子，多乘一個(gè)根本DFT因子WN2.3序列倒序?yàn)榱吮ＷC運(yùn)算輸出的X(k)按順序排列，要求序列x(n)倒序輸入，即在運(yùn)算前要先對(duì)輸入的序列進(jìn)行位序顛倒。如果總點(diǎn)數(shù)為的x(n)的順序數(shù)是用M位二進(jìn)制數(shù)表示，那么倒序數(shù)只需將順序數(shù)的二進(jìn)制位倒置即可，按照這一規(guī)律用硬件電路和匯編語言很容易產(chǎn)生倒序數(shù)。3MATLAB程序?qū)崿F(xiàn)MATLAB提供的fft函數(shù)是一個(gè)計(jì)算DFT的智能程數(shù)據(jù)倒序程序框圖序，能自動(dòng)選擇快速算法進(jìn)行DFT運(yùn)算，由于它是一個(gè)內(nèi)建函數(shù)，用type命令看不到程序代碼。MATLAB等高級(jí)語言實(shí)現(xiàn)倒序時(shí)，直接倒置二進(jìn)制數(shù)位的方法不可取，還須找出產(chǎn)生倒序的十進(jìn)制規(guī)律。將十進(jìn)制順序數(shù)用I表示，與之對(duì)應(yīng)的二進(jìn)制數(shù)用IB表示。十進(jìn)制倒序數(shù)用J表示，與之對(duì)應(yīng)的二進(jìn)制數(shù)用JB表示。JB是IB的位倒置結(jié)果，十進(jìn)制順序數(shù)I增加1，相當(dāng)于IB最低位加1且逢2向高位進(jìn)1，即相當(dāng)于JB最高位加1且逢2向低位進(jìn)1。JB的變化規(guī)律反映到J的變化分二種情況：如果JB的最高位是0，那么直接由加1得到下一個(gè)倒序值；如果JB的最高位是1，那么要先將最高位變0，再在次高位加1。但次高位加1時(shí)，同樣要判斷0、1值，如果是0，那么直接加1，否那么要先將次高位變0，再判斷下一位。依此類推，直到完成最高位加1，逢2向右進(jìn)位的運(yùn)算。利用這一算法可按順序數(shù)I的遞增順序，依次求得與之對(duì)應(yīng)的倒序數(shù)J。為了節(jié)省內(nèi)存，數(shù)據(jù)倒序可原址進(jìn)行，當(dāng)I=J時(shí)不需要交換，當(dāng)I≠J時(shí)需要交換數(shù)據(jù)。另外，為了防止再次調(diào)換前面已經(jīng)調(diào)換過的一對(duì)數(shù)據(jù)，只對(duì)I<J的情況進(jìn)行數(shù)據(jù)交換即可實(shí)現(xiàn)數(shù)據(jù)倒序操作。數(shù)據(jù)倒序程序框圖實(shí)驗(yàn)內(nèi)容及步驟%基2DIT-FFT運(yùn)算的MATLAB程序clc;closeall;clear;formatcompact;%輸入數(shù)據(jù)并計(jì)算常量xn=[0,1,2,3,4,5,6,7];%可取任意序列M=nextpow2(length(xn)),N=2^M,form=0:N/2-1;%旋轉(zhuǎn)因子指數(shù)范圍WN(m+1)=exp(-j*2*pi/N)^m;%計(jì)算旋轉(zhuǎn)因子endA=[xn,zeros(1,N-length(xn))];%數(shù)據(jù)輸入disp('輸入到各存儲(chǔ)單元的數(shù)據(jù):'),disp(A);%數(shù)據(jù)倒序操作J=0;%給倒序數(shù)賦初值forI=0:N-1;%按序交換數(shù)據(jù)和算倒序數(shù)ifI<J;%條件判斷及數(shù)據(jù)交換T=A(I+1);A(I+1)=A(J+1);A(J+1)=T;end%算下一個(gè)倒序數(shù)K=N/2;whileJ>=K;J=J-K;K=K/2;endJ=J+K;enddisp('倒序后各存儲(chǔ)單元的數(shù)據(jù):'),disp(A);%分級(jí)按序依次進(jìn)行蝶形運(yùn)算forL=1:M;%分級(jí)計(jì)算disp('運(yùn)算級(jí)次:'),disp(L);B=2^(L-1);forR=0:B-1;%各級(jí)按序蝶算P=2^(M-L)*R;forK=R:2^L:N-2;%每序依次計(jì)算T=A(K+1)+A(K+B+1)*WN(P+1);A(K+B+1)=A(K+1)-A(K+B+1)*WN(P+1);A(K+1)=T;endenddisp('本級(jí)運(yùn)算后各存儲(chǔ)單元的數(shù)據(jù):'),disp(A);enddisp('輸出各存儲(chǔ)單元的數(shù)據(jù):'),Xk=A,disp('調(diào)用fft函數(shù)運(yùn)算的結(jié)果:'),fftxn=fft(xn,N),作業(yè)快速傅里葉變換-基2時(shí)間抽取FFT算法matlab實(shí)現(xiàn)要求:輸入信號(hào)x(t)=0.6sin(200πt)+sin(400πt)+0.3sin(800πt)。1、實(shí)現(xiàn)基2時(shí)間抽取的FFT算法的1024點(diǎn)，4096點(diǎn)變換，并與matlab自帶函數(shù)進(jìn)行比擬，包括運(yùn)算時(shí)間和精度實(shí)驗(yàn)報(bào)告要求1、給出快速傅里葉變換程序；2、對(duì)給定信號(hào)進(jìn)行頻譜分析；3、給出與matlab自帶函數(shù)比擬結(jié)果實(shí)驗(yàn)四：1.快速傅里葉變換程序代碼：functiony=myfft(xr,n)p=0:n-1;nu=log2(n);p1=p;b=zeros(1,n);fort=1:nu;p2=floor(p1/2);b=b*2+(p1-2*p2);p1=p2;end;yr(p+1)=xr(b+1);xr=yr;t=0:n/2-1;forv=0:n/2-1;w=exp(-2*i*pi*t/n);end;yr(p+1)=xr(b+1);xr=yr;t=0:n/2-1;forv=0:n/2-1;w=exp(-2*i*pi*t/n);end;form=1:nu;h=2^(m-1);k=1;while(k<n+1)fort=1:h;y=bitshift(k-1,nu-m,nu)+1;xch(k)=xr(k)+w(y)*xr(k+h);k=k+1;end;fort=1:h;y=bitshift(k-1-h,nu-m,nu)+1;xch(k)=xr(k-h)-xr(k)*w(y);k=k+1;end;end;xr=xch;end;y=xr2.輸入信號(hào)x(t)=0.6sin(200πt)+sin(400πt)+0.3sin(800πt)的頻譜分析N=1024的傅里葉變換fs=100;N=1024;n=0:N-1;t=n/fs;x=0.6*sin(200*pi*t)+sin(400*pi*t)+0.3*sin(800*pi*t);y=myfft(x,N);mag=abs(y);f=n*fs/N;subplot(2,2,1),plot(f,mag);xlabel('頻率/Hz');ylabel('振幅');title('N=128');gridon;Subplot(2,2,2);plot(f,angle(y));xlabel('頻率/Hz');ylabel('相位');title('N=128');gridon;N=1024的傅里葉自帶函數(shù)fs=100;N=1024;n=0:N-1;t=n/fs;x=0.6*sin(200*pi*t)+sin(400*pi*t)+0.3*sin(800*pi*t);y=fft(x,N);mag=abs(y);f=(0:length(y)-1)'*fs/length(y);subplot(221);plot(f,mag);xlabel('Frequence(Hz)');ylabel('Magnitude');title('N=1024');grid;subplot(222);plot(f,angle(y));xlabel('Frequence(Hz)');ylabel('Magnitude');title('N=1024');grid;N=4096的傅里葉變換fs=100;N=4096;n=0:N-1;t=n/fs;x=0.6*sin(200*pi*t)+sin(400*pi*t)+0.3*sin(800*pi*t);y=myfft(x,N);mag=abs(y);f=n*fs/N;subplot(2,2,1),plot(f,mag);xlabel('頻率/Hz');ylabel('振幅');title('N=4096');gridon;Subplot(2,2,2);plot(f,angle(y));xlabel('頻率/Hz');ylabel('相位');title('N=4096');gridon;N=4096的傅里葉變換自帶函數(shù)fs=100;N=4096;n=0:N-1;t=n/fs;x=0.6*sin(200*pi*t)+sin(400*pi*t)+0.3*sin(800*pi*t);y=fft(x,N);mag=abs(y);f=(0:length(y)-1)'*fs/length(y);subplot(221);plot(f,mag);xlabel('Frequence(Hz)');ylabel('Magnitude');title('N=4096');grid;subplot(222);plot(f,angle(y));xlabel('Frequence(Hz)');ylabel('Magnitude');title('N=4096');grid;3.與matlab自帶函數(shù)進(jìn)行比擬結(jié)果(運(yùn)算時(shí)間和精度):myfft()函數(shù)精度與matlab自帶的函數(shù)相比所用的時(shí)