1、實(shí)驗(yàn)四 信號(hào)的頻域分析實(shí)驗(yàn)?zāi)康模?、掌握FFT和IFFT；2、掌握功率密度普；3、掌握數(shù)字濾波；一、實(shí)驗(yàn)講義信號(hào)的頻域分析在實(shí)際中有著廣泛的應(yīng)用，是LTI系統(tǒng)分析的基礎(chǔ)。對于連續(xù)信號(hào)的分析，如周期信號(hào)的傅里葉級(jí)數(shù)、非周期信號(hào)的傅里葉變換等，都需要利用其相應(yīng)的公式進(jìn)行分解或變換，以便獲得傅里葉級(jí)數(shù)和傅里葉變換的表達(dá)形式，然后利用繪圖工具繪制其頻譜。對于連續(xù)信號(hào)而言，MATLAB無法提供專門的函數(shù)進(jìn)行頻譜分析和處理。因此，本節(jié)只討論離散信號(hào)和離散系統(tǒng)的頻譜分析，并介紹頻域分析的應(yīng)用（以數(shù)字濾波為例）。1 離散傅里葉變換及其逆變換傅里葉變換可以將時(shí)域信號(hào)轉(zhuǎn)換為頻域信號(hào)，以便分析信號(hào)的頻域特性，其逆變

2、換則把頻域信號(hào)轉(zhuǎn)換為時(shí)域信號(hào)。傅里葉變換的原理是把一個(gè)時(shí)域信號(hào)分解成用不同頻率的正弦信號(hào)（或復(fù)信號(hào)）線性組合的形式，這樣時(shí)域信號(hào)所包含的頻率成分就一目了然。離散傅里葉變換DFT是數(shù)字信號(hào)分析的主要工具，快速傅里葉變換是對離散傅里葉變換DFT進(jìn)行快速計(jì)算的有效算法。在MATLAB中，利用函數(shù)fft和ifft分別計(jì)算一維信號(hào)的離散傅里葉變換和其逆變換。對于二維信號(hào)，使用的離散傅里葉變換和離散傅里葉變換逆變換函數(shù)則為fft2和ifft2。函數(shù)fft和ifft的調(diào)用格式為：y=fft(x,n)；Y=ifft(X,n)說明：（1）x和X分別為待變換的輸入向量，x為時(shí)域信號(hào)，X為頻域信號(hào)。（2）n表示進(jìn)行

3、變換的點(diǎn)數(shù)，可以缺省，若輸入的序列比n短，則fft和ifft用0填充序列，使其長度為n；若輸入的序列比n長，則截短輸入序列。（3）n缺省時(shí)，變換的點(diǎn)數(shù)為輸入序列的長度。（4）n為2的冪數(shù)時(shí)，計(jì)算速度最快。例1 對信號(hào)進(jìn)行傅里葉變換，然后對變換后的序列進(jìn)行傅里葉逆變換，并繪制它們的圖像。N=512; % N為采樣點(diǎn)數(shù)T=1; %T為采樣時(shí)間終點(diǎn)t=linspace(0,T,N); % 給出N個(gè)采樣時(shí)間ti(I=1:N)x=sin(2*pi*50*t)+5*cos(2*pi*120*t); % 求各采樣點(diǎn)樣本值xfigure(1),plot(x(1:30); %繪制輸入信號(hào)的部分波形dt=t(2)

4、-t(1); % 采樣周期f=1/dt; % 采樣頻率(Hz)X=fft(x); % 計(jì)算x的快速傅立葉變換XF=X(1:N/2+1); % F(k)=X(k)(k=1:N/2+1)f=f*(0:N/2)/N; % 使頻率軸f從零開始figure(2),plot(f,abs(F),'-*') % 繪制振幅-頻率圖xlabel('Frequency');ylabel('|F(k)|');z=ifft(X) ; %傅里葉逆變換Z=real(z); %逆變換之后的信號(hào)是復(fù)信號(hào)，去其實(shí)部即可figure(3),plot(Z(1:30); %繪制逆變換之后

5、的信號(hào)程序運(yùn)行結(jié)果如圖1、圖2和圖3所示。圖1 輸入信號(hào)的時(shí)域圖圖2 輸入信號(hào)的頻譜圖3 逆傅里葉變換后的信號(hào)可以看出，1和2是完全一致的，說明傅里葉變換和傅里葉逆變換是嚴(yán)格對應(yīng)的。從圖2可知，頻譜在50Hz和120Hz有兩個(gè)沖激，說明原時(shí)域信號(hào)中同時(shí)含有頻率為50Hz和120Hz的兩個(gè)周期信號(hào)。2 信號(hào)的功率密度譜信號(hào)的功率密度譜是信號(hào)功率譜密度值在頻率范圍內(nèi)的分布，具體表現(xiàn)就是不同頻段上波形幅值的變化。功率譜是進(jìn)行頻率分析的基礎(chǔ)，在許多工程方面有著重要的意義。在MATLAB中，求取信號(hào)的功率密度譜用函數(shù)psd來實(shí)現(xiàn)，該函數(shù)的調(diào)用格式為：Pxx,f=psd(Xn,nfft,Fs,window

6、,noverlap)說明：（1）Pxx為輸入信號(hào)Xn的功率譜密度數(shù)值序列，f為與Pxx對應(yīng)的頻率序列。（2）Xn為輸入的時(shí)域信號(hào)，nfft為計(jì)算FFT的單位寬度，F(xiàn)S為采樣頻率，window為聲明窗函數(shù)的類型，noverlap是處理Xn混疊的點(diǎn)數(shù)。 （3）使用psd函數(shù)計(jì)算信號(hào)的功率密度譜的基本原理是：把輸入時(shí)域信號(hào)Xn分成許多連續(xù)的區(qū)域，對每個(gè)區(qū)域加窗（窗函數(shù)的類型由window設(shè)置），然后做FFT計(jì)算，每兩個(gè)相鄰的區(qū)域之間有個(gè)重疊的區(qū)域（大小由noverlap設(shè)置）。把每個(gè)區(qū)域計(jì)算所得的序列先做模計(jì)算，然后對所有區(qū)域的序列進(jìn)行求和，最后用這個(gè)和除以計(jì)算過的功率分布范圍就獲得到了信號(hào)的功率密

7、度譜。例2 在時(shí)域信號(hào)摻入隨機(jī)噪聲，并繪制該信號(hào)的功率密度譜。t=0:0.001:1;x=sin(2*pi*60*t)+2*sin(2*pi*160*t)+randn(size(t) ; %信號(hào)加噪聲nfft=256; %設(shè)置nfft值Fs=1000; %設(shè)置采樣頻率window=hanning(nfft); %設(shè)置窗函數(shù)為漢寧窗，窗寬為256noverlap=128; %混疊寬度為128點(diǎn)Pxx,f=psd(x,nfft,Fs,window,noverlap); %計(jì)算功率密度譜plot(t(1:200),x(1:200); %繪制輸入信號(hào)figure,plot(f,Pxx) %繪制功率密度

8、譜程序運(yùn)行結(jié)果如圖4所示。圖4 信號(hào)的功率密度譜從圖4可以看出，經(jīng)過psd函數(shù)的計(jì)算，隨機(jī)噪聲的能量被大大地弱化，而周期信號(hào)的能量得到了加強(qiáng)。3 信號(hào)的互相關(guān)功率密度譜信號(hào)的相關(guān)分為自相關(guān)和互相關(guān)兩種類型，分別說明一個(gè)信號(hào)與自己或另外一個(gè)信號(hào)之間在頻域上的相似性。信號(hào)的互相關(guān)功率密度譜（如果兩個(gè)信號(hào)完全相同，則為自相關(guān)功率密度譜）在故障診斷和狀態(tài)預(yù)測等方面有廣泛的應(yīng)用，例如檢測振動(dòng)信號(hào)中是否有周期信號(hào)、檢測零部件裂縫的位置及孔洞的大小等。在MATLAB中，求取信號(hào)的互相關(guān)功率密度譜用函數(shù)csd來實(shí)現(xiàn)，該函數(shù)的調(diào)用格式為：Pxy,f=csd(x,y,nfft,Fs,window,noverlap

9、)說明：（1）x和y分別為兩個(gè)時(shí)域的輸入信號(hào)，nfft為計(jì)算FFT的單位寬度，F(xiàn)s為采樣頻率，window是指窗函數(shù)的類型，noverlap是處理Xn混疊的點(diǎn)數(shù)。（2）Pxy為互相關(guān)功率密度譜的數(shù)值序列，f為相應(yīng)的頻率軸。（3）該函數(shù)與功率密度譜函數(shù)psd的原理基本相似。（4）該函數(shù)的主要作用是用來辨認(rèn)信號(hào)中的周期成分。例3 求一個(gè)含有頻率為100Hz正弦波的隨機(jī)噪聲信號(hào)的自相關(guān)功率密度譜，并繪制其圖像。t=0:0.001:1;x=sin(2*pi*100*t)+randn(size(t); %含有頻率為100Hz正弦波的隨機(jī)信號(hào)nfft=256; %設(shè)置nfft值Fs=1000; %設(shè)置采樣

10、頻率window=hanning(256); %設(shè)置窗函數(shù)為漢寧窗，窗寬為256noverlap=128; %混疊寬度為128點(diǎn)Pxx,f=csd(x,x,nfft,Fs,window,noverlap); %計(jì)算自相關(guān)功率密度譜subplot(1,2,1),plot(t(1:200),x(1:200);title(含有正弦信號(hào)的隨機(jī)噪聲); %繪制輸入信號(hào)subplot(1,2,2),plot(f,Pxx); title(功率密度譜) %繪制功率密度譜程序運(yùn)行結(jié)果如圖5所示。 圖5含有正弦信號(hào)的隨機(jī)噪聲及其自相關(guān)功率密度譜由圖5可知，自相關(guān)功率密度值都大于0，事實(shí)上自相關(guān)功率譜還是偶函數(shù)。從

11、時(shí)域圖像中很難看出此輸入信號(hào)中含有周期成分，更不要說確定周期成分的周期，但是從自相關(guān)譜中很容易就看出這個(gè)雜亂無章的信號(hào)中含有頻率為100Hz的周期成分。例4 試?yán)L制兩個(gè)時(shí)域信號(hào)之間的互相關(guān)功率譜。t=0:0.001:1;x=sin(2*pi*300*t)+2*sin(2*pi*160*t)+randn(size(t);y=sin(2*pi*100*t)+2*sin(2*pi*180*t)+randn(size(t);nfft=256; %設(shè)置nfft值Fs=1000; %設(shè)置采樣頻率window=hanning(256); %設(shè)置窗函數(shù)為漢寧窗，窗寬為256noverlap=128; %混疊寬

12、度為128點(diǎn)Pxy,f=csd(x,y,nfft,Fs,window,noverlap); %計(jì)算互相關(guān)功率密度譜plot(f,Pxy) %繪制功率密度譜圖像程序運(yùn)行結(jié)果如圖6所示：圖6 兩個(gè)信號(hào)的互相關(guān)功率譜兩個(gè)信號(hào)分別含有頻率為160Hz、300Hz和100Hz、180Hz的正弦波，從信號(hào)的互相關(guān)功率譜可以明顯觀察到這些頻率的周期信號(hào)。從圖9-20中可知，互相關(guān)功率譜不全為正值。4 數(shù)字濾波濾波器設(shè)計(jì)的目的是對信號(hào)進(jìn)行數(shù)據(jù)序列的頻率變換，去掉不需要的頻率成分，變換成輸出序列。一般有低通（濾除高頻成分）、高通（濾除低頻成分）、帶通（濾除通頻帶兩端的頻率成分）和帶阻（濾除聲明的阻頻帶部分）等幾

13、種濾波器。一般來講，數(shù)字濾波器的輸出序列y(n)的Z變換Y(z)與濾波器的輸入序列x(n)的Z變換X(z)是相關(guān)聯(lián)的，常常表述如下： （1）其中，H(z)是濾波器的傳遞函數(shù)，常量a(i)和 b(i)是濾波器的系數(shù)，而濾波器的階次是m和n中的較大值，濾波器的系數(shù)分別以兩個(gè)行向量a和b的形式被存儲(chǔ)起來。在MATLAB中，濾波運(yùn)算用函數(shù)filter來實(shí)現(xiàn)，該函數(shù)的調(diào)用格式為：y=filter(b,a,x)其中，a、b是濾波器的系數(shù)向量。濾波器的輸出序列y和輸入序列x的長度（即采樣點(diǎn)數(shù)）是相等的。在MATLAB中，不同的濾波器類型是通過不同的函數(shù)計(jì)算來實(shí)現(xiàn)的。濾波器有很多類型，其對應(yīng)的函數(shù)也各不相同，

14、下面介紹其中的一類濾波器即Butterworth濾波器。當(dāng)計(jì)算出Butterworth濾波器的系數(shù)后，即可用之來實(shí)現(xiàn)信號(hào)的濾波。例5 在正弦信號(hào)中加入隨機(jī)噪聲，采樣頻率為1000Hz，請濾除其中頻率為30Hz以上的噪聲。t=0:0.001:0.5;x=sin(50*pi*t)+randn(size(t); %生成輸入序列b,a=butter(10,30/500); %計(jì)算濾波器系數(shù)y=filter(b,a,x); %進(jìn)行數(shù)字濾波subplot(1,2,1),plot(t,x); axis(0 0.5 -4 4) ;title(含噪信號(hào)); %繪制濾波前信號(hào)的圖像subplot(1,2,2),p

15、lot(t,y); axis(0 0.5 -2 2) ; title(濾波結(jié)果); %繪制濾波后信號(hào)的圖像程序運(yùn)行結(jié)果如圖7所示： 圖7 輸入序列及其濾波結(jié)果程序說明：（1）butter(10,30/500)是一個(gè)濾波器函數(shù)，它的第一個(gè)參數(shù)表示這個(gè)濾波器的階數(shù)是10，此處的階數(shù)為估計(jì)值。在不同的情況下應(yīng)當(dāng)選取什么樣的階次，MATLAB里面有專門的函數(shù)來計(jì)算這個(gè)數(shù)值。butter函數(shù)的階次估計(jì)函數(shù)為buttord，其具體的用法此處不予介紹。（2）butter(10,30/500)的第二個(gè)參數(shù)是濾波器的標(biāo)準(zhǔn)化頻率。在MATLAB中，包括butter在內(nèi)的所有濾波器函數(shù)的截止頻率都采用標(biāo)準(zhǔn)化頻率，即以Hz為單位的截止頻率除以采樣頻率的一半，所以標(biāo)準(zhǔn)頻率是一個(gè)分布在區(qū)間0,1的比值。（3）從圖7可知，信號(hào)中的高頻成分被濾除了，但是由于其中仍然含