實(shí)驗(yàn)三離散傅立葉變換

實(shí)驗(yàn)三離散傅立葉變換第一頁，共四十三頁，2022年，8月28日一、實(shí)驗(yàn)?zāi)康?/p>

加深對(duì)離散傅立葉變換(DFT)的理解。掌握利用MATLAB語言進(jìn)行離散傅立葉變換和逆變換的方法。加深對(duì)離散傅立葉變換基本性質(zhì)的理解。掌握離散傅立葉變換快速算法的應(yīng)用。第二頁，共四十三頁，2022年，8月28日二、實(shí)驗(yàn)原理及方法

建立以時(shí)間t為自變量的“信號(hào)”與以頻率f為自變量的“頻率函數(shù)”(頻譜)之間的某種變換關(guān)系。所以“時(shí)間”或“頻率”取連續(xù)還是離散值,就形成各種不同形式的傅里葉變換對(duì)。傅里葉變換第三頁，共四十三頁，2022年，8月28日四種不同傅里葉變換對(duì)傅里葉級(jí)數(shù)(FS)：連續(xù)時(shí)間,離散頻率的傅里葉變換。周期連續(xù)時(shí)間信號(hào)傅里葉級(jí)數(shù)(FS)得到非周期離散頻譜密度函數(shù)。傅里葉變換(FT)：連續(xù)時(shí)間,連續(xù)頻率的傅里葉變換。非周期連續(xù)時(shí)間信號(hào)通過連續(xù)付里葉變換(FT)得到非周期連續(xù)頻譜密度函數(shù)。離散時(shí)間的傅里葉變換(DTFT):離散時(shí)間,連續(xù)頻率的傅里葉變換。非周期離散的時(shí)間信號(hào)(單位園上的Z變換(DTFT))得到周期性連續(xù)的頻率函數(shù)。離散傅里葉變換(DFT):離散時(shí)間,離散頻率的傅里葉變換。第四頁，共四十三頁，2022年，8月28日上面討論的前三種傅里葉變換對(duì),都不適用在計(jì)算機(jī)上運(yùn)算,因?yàn)橹辽僭谝粋€(gè)域(時(shí)域或頻域)中,函數(shù)是連續(xù)的。因?yàn)閺臄?shù)字計(jì)算角度我們感興趣的是時(shí)域及頻域都是離散的情況,這就是第四種離散傅里葉變換。第五頁，共四十三頁，2022年，8月28日離散傅里葉級(jí)數(shù)(DFS)離散時(shí)間序列x(n)滿足x(n)=x(n+rN)，稱為離散周期序列，其中N為周期，x(n)為主值序列。由傅立葉分析知道周期函數(shù)可由復(fù)指數(shù)的線性組合疊加得到。其頻率為基本頻率的倍數(shù)。從離散時(shí)間傅立葉變換的頻率周期性，我們知道諧波次數(shù)是有限的，其頻率為周期序列可表示成：第六頁，共四十三頁，2022年，8月28日其中叫做離散傅立葉級(jí)數(shù)系數(shù)，也稱為周期序列的頻譜，可由下式表示注意也是一個(gè)基本周期為N的周期序列。上面兩式稱為周期序列的傅立葉級(jí)數(shù)變換對(duì)。令表示復(fù)指數(shù)，可以得到以下：第七頁，共四十三頁，2022年，8月28日例：求出下面周期序列的DFSx(n)={……,0,1,2,3,0,1,2,3,0,1,2,3,……}基本周期為N=4，WN=W4=-j，

因而第八頁，共四十三頁，2022年，8月28日MATLAB實(shí)現(xiàn)矩陣-向量相乘運(yùn)算來實(shí)現(xiàn)。由于和均為周期函數(shù)，周期為N，可設(shè)和代表序列和的主值區(qū)間序列，則前面的兩個(gè)表達(dá)式可寫成：式中，矩陣WN為方陣——DFS矩陣。第九頁，共四十三頁，2022年，8月28日利用MATLAB實(shí)現(xiàn)傅立葉級(jí)數(shù)計(jì)算編寫函數(shù)實(shí)現(xiàn)DFS計(jì)算functionxk=dfs(xn,N)n=[0:1:N-1];%n的行向量k=n;%k的行向量WN=exp(-j*2*pi/N);%WN因子nk=n’*k;%產(chǎn)生一個(gè)含nk值的N乘N維矩陣WNnk=WN.^nk;%DFS矩陣xk=xn*WNnk;%DFS系數(shù)行向量第十頁，共四十三頁，2022年，8月28日例：xn=[0,1,2,3]，N=4xn=[0,1,2,3];N=4;xk=dfs(xn,N)’第十一頁，共四十三頁，2022年，8月28日逆運(yùn)算IDFSfunctionxn=idfs(xk,N)n=[0:1:N-1];k=n;WN=exp(-j*2*pi/N);nk=n’*k;WNnk=WN.^(-nk);xn=（xk*WNnk）/N;第十二頁，共四十三頁，2022年，8月28日xn=idfs(xk',4)x=xn'第十三頁，共四十三頁，2022年，8月28日周期重復(fù)次數(shù)對(duì)序列頻譜的影響理論上講，周期序列不滿足絕對(duì)可積條件，要對(duì)周期序列進(jìn)行分析，可以先取K個(gè)周期進(jìn)行處理，然后讓K無限增大，研究其極限情況。這樣可以觀察信號(hào)序列由非周期到周期變換時(shí)，頻譜由連續(xù)譜逐漸向離散譜過渡的過程。第十四頁，共四十三頁，2022年，8月28日例：已知一個(gè)矩形序列的脈沖寬度占整個(gè)周期的1/2，一個(gè)周期的采樣點(diǎn)數(shù)為10，用傅立葉級(jí)數(shù)變換求信號(hào)的重復(fù)周期數(shù)分別為1、4、7、10時(shí)的幅度頻譜。MATLAB程序：xn=[ones(1,5),zeros(1,5)];Nx=length(xn);Nw=1000;dw=2*pi/Nw;k=floor((-Nw/2+0.5):(Nw/2+0.5));forr=0:3;K=3*r+1;nx=0:(K*Nx-1);x=xn(mod(nx,Nx)+1);Xk=x*(exp(-j*dw*nx'*k))/K;subplot(4,2,2*r+1);stem(nx,x)axis([0,K*Nx-1,0,1.1]);ylabel('x(n)');subplot(4,2,2*r+2);plot(k*dw,abs(Xk))axis([-4,4,0,1.1*max(abs(Xk))]);ylabel('X(k)');end第十五頁，共四十三頁，2022年，8月28日從上圖可以看出，信號(hào)序列的周期數(shù)越多，則頻譜越是向幾個(gè)頻點(diǎn)集中，當(dāng)信號(hào)周期數(shù)趨于無窮大時(shí)，頻譜轉(zhuǎn)化為離散譜。第十六頁，共四十三頁，2022年，8月28日離散傅立葉變換（DFT）有限長序列x(n)表示為x(n)是非周期序列，但可以理解為周期序列的主值序列。由離散傅立葉級(jí)數(shù)DFS和IDFS引出有限長序列的離散傅立葉正、逆變換關(guān)系式第十七頁，共四十三頁，2022年，8月28日有限長序列傅立葉變換定義式為：

比較正、逆變換的定義式可以看出，只要把DFT公式中的系數(shù)改為，并最后乘以1/N，那么，DFT的計(jì)算程序就可以用來計(jì)算IDFT。第十八頁，共四十三頁，2022年，8月28日DFT與DFS的關(guān)系比較兩者的變換對(duì)，可以看出兩者的區(qū)別僅僅是將周期序列換成了有限長序列。有限長序列x(n)可以看作是周期序列的一個(gè)周期；反之周期序列可以看作是有限長序列x(n)以N為周期的周期延拓。由于公式非常相似，在程序編寫上也基本一致。第十九頁，共四十三頁，2022年，8月28日例：已知序列x(n)=[0,1,2,3,4,5,6,7]，求x(n)的DFT和IDFT，畫出序列傅立葉變換的幅度和相位圖，并將原圖像與逆變換圖像進(jìn)行比較。N=8;xn=0:N-1;n=0:N-1;xk=dft(xn,N);x=idft(xk,N);subplot(2,2,1);stem(n,xn)subplot(2,2,2);stem(n,abs(x))subplot(2,2,3);stem(n,abs(xk))subplot(2,2,4);stem(n,angle(xk))第二十頁，共四十三頁，2022年，8月28日第二十一頁，共四十三頁，2022年，8月28日三、快速傅立葉變換有限長序列通過離散傅里葉變換(DFT)將其頻域離散化成有限長序列.但其計(jì)算量太大(與N的平方成正比）,很難實(shí)時(shí)地處理問題,因此引出了快速傅里葉變換(FFT)。FFT并不是一種新的變換形式,它只是DFT的一種快速算法.并且根據(jù)對(duì)序列分解與選取方法的不同而產(chǎn)生了FFT的多種算法.第二十二頁，共四十三頁，2022年，8月28日DFT的快速算法—FFT是數(shù)字信號(hào)處理的基本方法和基本技術(shù)，是必須牢牢掌握的。時(shí)間抽選FFT算法的理論推導(dǎo)和流圖詳見《數(shù)字信號(hào)處理》教材。該算法遵循兩條準(zhǔn)則：（1）對(duì)時(shí)間奇偶分；（2）對(duì)頻率前后分。這種算法的流圖特點(diǎn)是：（1）基本運(yùn)算單元都是蝶形

任何一個(gè)長度為N=2M的序列，總可通過M次分解最后成為2點(diǎn)的DFT計(jì)算。如圖所示：第二十三頁，共四十三頁，2022年，8月28日WNk稱為旋轉(zhuǎn)因子計(jì)算方程如下：Xm+1(p)=Xm(p)+WNkXm(q)Xm+1(q)=Xm(p)-WNkXm(q)第二十四頁，共四十三頁，2022年，8月28日（2）同址（原位）計(jì)算這是由蝶形運(yùn)算帶來的好處，每一級(jí)蝶形運(yùn)算的結(jié)果Xm+1(p)無須另外存儲(chǔ)，只要再存入Xm(p)中即可，Xm+1(q)亦然。這樣將大大節(jié)省存儲(chǔ)單元。（3）變址計(jì)算輸入為“混序”（碼位倒置）排列，輸出按自然序排列，因而對(duì)輸入要進(jìn)行“變址”計(jì)算（即碼位倒置計(jì)算）。“變址”實(shí)際上是一種“整序”的行為，目的是保證“同址”。第二十五頁，共四十三頁，2022年，8月28日FFT的應(yīng)用凡是利用付里葉變換來進(jìn)行分析、綜合、變換的地方，都可以利用FFT算法來減少其計(jì)算量。FFT主要應(yīng)用在1、快速卷積2、快速相關(guān)3、頻譜分析第二十六頁，共四十三頁，2022年，8月28日快速傅立葉變換的MATLAB實(shí)現(xiàn)提供fft函數(shù)計(jì)算DFT格式

X=fft(x)

X=fft(x,N)如果x的長度小于N，則在其后填零使其成為N點(diǎn)序列，反之對(duì)x進(jìn)行截?cái)?，若省略變量N，則DFT的長度即為x的長度。如果N為2的冪，則得到高速的基-2FFT算法；若N不是2的乘方，則為較慢的混合算法。如果x是矩陣，則X是對(duì)矩陣的每一列向量作FFT。第二十七頁，共四十三頁，2022年，8月28日快速傅立葉逆變換（IFFT）函數(shù)調(diào)用格式

y=ifft(x)y=ifft(x,N)當(dāng)N小于x長度時(shí)，對(duì)x進(jìn)行截?cái)?，?dāng)N大于x長度時(shí)，對(duì)x進(jìn)行補(bǔ)零。第二十八頁，共四十三頁，2022年，8月28日fftshift函數(shù)功能：對(duì)fft的輸出進(jìn)行重新排列，將零頻分量移到頻譜的中心。調(diào)用格式

y=fftshift(x)

當(dāng)x為向量時(shí)，fftshift(x)直接將x中左右兩半交換而產(chǎn)生y。當(dāng)x為矩陣時(shí)，fftshift(x)直接將x中左右、上下進(jìn)行交換而產(chǎn)生y。第二十九頁，共四十三頁，2022年，8月28日由題目可得x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)fs=100N=128/1024例：已知信號(hào)由15Hz幅值0.5的正弦信號(hào)和40Hz幅值2的正弦信號(hào)組成，數(shù)據(jù)采樣頻率為100Hz,試?yán)L制N=128點(diǎn)DFT的幅頻圖。第三十頁，共四十三頁，2022年，8月28日fs=100;N=128;n=0:N-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,N);f=(0:length(y)-1)'*fs/length(y);mag=abs(y);stem(f,mag);title(‘N=128點(diǎn)’)第三十一頁，共四十三頁，2022年，8月28日第三十二頁，共四十三頁，2022年，8月28日利用FFT進(jìn)行功率譜的噪聲分析已知帶有測(cè)量噪聲信號(hào)其中f1=50Hz,f2=120Hz,為均值為零、方差為1的隨機(jī)信號(hào)，采樣頻率為1000Hz，數(shù)據(jù)點(diǎn)數(shù)N=512。試?yán)L制信號(hào)的功率譜圖。第三十三頁，共四十三頁，2022年，8月28日t=0:0.001:0.6;x=sin(2*pi*50*t)+sin(2*pi*120*t);y=x+2*randn(1,length(t));Y=fft(y,512);P=Y.*conj(Y)/512;%求功率f=1000*(0:255)/512;subplot(2,1,1);plot(y);subplot(2,1,2);plot(f,P(1:256));第三十四頁，共四十三頁，2022年，8月28日第三十五頁，共四十三頁，2022年，8月28日序列長度和FFT的長度對(duì)信號(hào)頻譜的影響。已知信號(hào)其中f1=15Hz,f2=40Hz,采樣頻率為100Hz.

在下列情況下繪制其幅頻譜。

Ndata=32,Nfft=32;Ndata=32,Nfft=128;第三十六頁，共四十三頁，2022年，8月28日fs=100;Ndata=32;Nfft=32;n=0:Ndata-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,Nfft);mag=abs(y);f=(0:length(y)-1)’*fs/length(y);subplot(2,1,1)plot(f(1:Nfft/2),mag(1:Nfft/2))title(‘Ndata=32,Nfft=32’)第三十七頁，共四十三頁，2022年，8月28日Nfft=128;n=0:Ndata-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,Nfft);mag=abs(y);f=(0:length(y)-1)’*fs/length(y);subplot(2,1,2)plot(f(1:Nfft/2),mag(1:Nfft/2))title(‘Ndata=32,Nfft=128’)第三十八頁，共四十三頁，2022年，8月28日第三十九頁，共四十三頁，2022年，8月28日線性卷積的FFT算法在MATLAB實(shí)現(xiàn)卷積的函數(shù)為CONV，對(duì)于N值較小的向量，這是十分有效的。對(duì)于N值較大的向量卷積可用FFT加快計(jì)算速度。由DFT性質(zhì)可知，若DFT[x1(n)]=X1