頻譜峰值搜索

1、實驗名稱：頻譜峰值搜索實驗?zāi)康模?、熟練掌握離散信號的DFT實現(xiàn)方法；2、熟練掌握Matlab實現(xiàn)DFT的方法，提高編程實踐能力；3、增強(qiáng)自我學(xué)習(xí)能力，查閱文獻(xiàn)搜索能力；4、掌握離散信號的時域與頻域的對應(yīng)關(guān)系。實驗原理：1.離散復(fù)正弦信號的DFTNj2二nkX(kfi)=x(n)eN(1)n=02、Matlab主要函數(shù)一fft(signal,N);signal:輸入信號，N：fft的點(diǎn)數(shù)函數(shù)的作用是對輸入信號做N點(diǎn)的DFTfftshift(fft(signal,N);將零頻點(diǎn)移到頻譜的中間3譜峰搜索算法采用一維黃金分割精搜算法5取代分級搜索過程中的遞歸精搜。對于一維黃金分割精搜算法，若函數(shù)f(

2、x)有且僅有一個極大值位于區(qū)間(a,c)上，則有當(dāng)a<b<c時，f(b)af(a)且f(b)af(c)。此時若在區(qū)間(b,c)上選取一點(diǎn)x,當(dāng)f(b)Af(x)時，則f(b)Af(a)且f(b)>f(x),即極大值點(diǎn)將位于三元點(diǎn)組a<b<x對應(yīng)的區(qū)間(a,x)上；否則，當(dāng)f(b)<f(x)時，則f(x)Af(b)且f(x)>f(b)>f(c),即極大值點(diǎn)將位于三元點(diǎn)組b<x<c對應(yīng)的區(qū)間(b,c)上。在這些三元點(diǎn)組中，其中間點(diǎn)對應(yīng)的函數(shù)值都是每一輪求解過程中的最大值.這一過程下圖所示，繼續(xù)對三元點(diǎn)組所對應(yīng)的區(qū)間進(jìn)行劃分，直到區(qū)間足夠小

3、，小到以至于該區(qū)間上任何一點(diǎn)都可以表示函數(shù)的極大值點(diǎn)。下面是一維黃金分割搜索算法的基本原理。給定三元點(diǎn)組a<b<c，假設(shè)b是a,cN間的一個分割wo即bac-b-二w,=1wcac-a(2)再假設(shè)一個試探點(diǎn)X位于區(qū)間（b,c）之間，且有(3)通過選取試探點(diǎn)X后，可以將極大值點(diǎn)壓縮到相對長度為0+2的區(qū)間（21）?；蛘呤窍鄬﹂L度為1W的區(qū)間（b,c）上?？紤]到搜索極大值的最壞情況，應(yīng)該使得下式成立：w+z=1-wz=1-2w（4）可見試探點(diǎn)X應(yīng)該選為點(diǎn)b關(guān)于區(qū)間（a，c）的對稱點(diǎn)。同樣應(yīng)該保證X是（b，c）之間的一個分割w0x-bz二w二wc-b1-w（5）把式（4）代入到式（5）中

4、得到如下的二次方程2w-3w+1=0（6）解得其根為（考慮到0<w<1,舍棄另一根）3-、/5/、w=0.3819660（7）2所以對于給定的三元點(diǎn)組所對應(yīng)的區(qū)間上，每次選取試探點(diǎn)都是位于較大的一段子區(qū)間上，并且距離原來中間點(diǎn)0.3189660的位置上。也就是說通過一次這樣的試探點(diǎn)的選取可以使得極大值點(diǎn)將位于原區(qū)間0.6180340的更小的區(qū)間上。繼續(xù)這樣的計算，那么區(qū)間將變得越來越小，小到以至于該區(qū)間上任何一點(diǎn)都可以表示函數(shù)的極大值點(diǎn)?？赼實驗步驟：1、設(shè)置輸入信號的參數(shù)以及DFT變換的點(diǎn)數(shù)根據(jù)要求，輸入信號的模擬頻率為f1=0.111111111,那么采樣頻率滿足fs>2

5、f1即可，為方便觀察頻率最大值位置，取fs=2Hz。給定DFT點(diǎn)數(shù)為1024點(diǎn)，而為了使的被觀察的頻譜峰值在頻譜圖的中央，將抽樣時間取在NN1t=,一fsfsfs的區(qū)間，采樣間隔為Ts=1/fs。其中N=512,滿足采樣點(diǎn)數(shù)為512M2=1024點(diǎn)。這樣得到輸入信號的表達(dá)式為signal=ej2;rf1nTs(8)2、對信號進(jìn)行DFT,并畫出頻譜圖。(1)在MATLAB中應(yīng)用fft(signal,N)對信號signal做N點(diǎn)的FFT;(2)分別應(yīng)用函數(shù)fftshift、abs又tDFT結(jié)果調(diào)整和取絕對值；(3)設(shè)置橫坐標(biāo)。根據(jù)fs和N的對應(yīng)關(guān)系，得到很坐標(biāo)的取值范圍是-1,1o在MATLAB中

6、設(shè)置為f=(1:2*N)-N)*(fs/(2*N);3、運(yùn)用一維黃金分割方法找出頻譜峰值。根據(jù)頻譜的峰值范圍，以及分割法的原理，設(shè)置個參數(shù)為：a=0.5,c=0.5,w=0.318966。創(chuàng)建計算相應(yīng)頻率點(diǎn)的幅值計算函數(shù)。根據(jù)DFT的計算定義，f、處的頻譜值為N1j-nkN1j"n"nN"lxx(fx)=£x(nTs)eN=£x(nTs)es(9)nz0n-0根據(jù)計算精度，將MATLAB計算精度設(shè)置為formatlong。并設(shè)計計算迭代次數(shù)的變量iterations。通過判別b點(diǎn)和x點(diǎn)的幅值大小來更新參數(shù)，參數(shù)更新如下、x(b)<x(f)

7、則：bTa,xTb,cTc(10)x(b)x(x)則：aTa,bTb,xTc實驗結(jié)果：1、運(yùn)行程序（程序見附錄）,得到頻譜圖如圖1所示輸入信號的DFT70SO5040302G10Q-10X;0.1133¥60.03-08-06-0.4-02a020406ftSfHz圖1復(fù)正弦信號的頻譜由于圖上顯示精度的原因，直接找到的最大值不是我們所需要的最大值，通過峰值搜索函數(shù)得到最大值。2、得到搜索結(jié)果為:fmax=0.111111111474039iterations=53。實驗結(jié)果分析：1、由于DFT的點(diǎn)數(shù)1024比較多，而頻譜范圍較小，所以離散的頻譜在圖上顯得像連續(xù)譜一樣。2、因為輸入信號

8、的模擬頻率的值為i=0.111111111,所以圖上離散的點(diǎn)上沒有顯示最大值點(diǎn)。3、從搜索的結(jié)果看出，最大值在滿足精度要求的情況下是正確的，說明一維黃金搜索方法在本實驗中是可行的。4、在a,c取值距離最大值較遠(yuǎn)的情況下經(jīng)過53次循環(huán)迭代可以得到最大值。說明該算法收斂比較快。附錄：clearall;clc;formatlongN=512;設(shè)置采樣點(diǎn)數(shù)為2N=1024fs=2;%設(shè)置采樣頻率為2Hzt=-N/fs:1/fs:N/fs-1/fs;%采樣時間序列f1=0.111111111;輸入信號的頻率signal=exp(i*2*pi*f1*t);輸入信號的采樣序列signalDtf=abs(ff

9、tshift(fft(signal);%對信號進(jìn)行DFT%畫出波形f=(1:2*N)-N)*(fs/(2*N);設(shè)置橫坐標(biāo)plot(f,20*log10(signalDtf);holdon;gridxlabel('fHz');ylabel('20log10(幅度)');title('輸入信號的DFT');%采用一維黃金分割精度算法%吉合算法的特點(diǎn)選擇a=-0.2,c=-0.2,w=0.3189660;%初始化幅值ampX=2;ampB=0;a=-0.5;c=0.5;w=0.3189660;z=1-2*w;iterations=0;%初始化迭代次數(shù)

10、whileabs(ampX-ampB)>1.0000e-0011;%設(shè)置收斂目標(biāo)%定義算法參數(shù)b=w*(c-a)+a;x=z*(c-a)+b;iterations=iterations+1;迭代次數(shù)更新%通過DFT定義計算a,b,c,x點(diǎn)的幅值ampB=Amplitude(b);ampA=Amplitude(a);ampC=Amplitude(c);ampX=Amplitude(x);compairBX=ampB>ampX;%比較b,x點(diǎn)的幅值大小ifcompairBX=1;a=a;b=b;c=x;elsecompairBX=0;a=b;b=x;c=c;endendamplitudeMax=xiterationsfunctionampx=Amplitude(x)%頻率幅值對應(yīng)幅值計算函數(shù)%計算任意頻率對應(yīng)點(diǎn)的幅度值%變量x為模擬平頻率