基于MATLAB的頻譜分析儀設(shè)計(jì)

1、基于MATLAB的信號(hào)頻譜分析儀的實(shí)現(xiàn)一、概述信號(hào)處理幾乎涉及到所有的工程技術(shù)領(lǐng)域，而頻譜分析又是信號(hào)處理中一個(gè)非常重要的分析手段。一般的頻譜分析都依靠傳統(tǒng)頻譜分析儀來完成，價(jià)格昂貴，體積龐大，不便于工程技術(shù)人員的攜帶。虛擬頻譜分析儀改變了原有頻譜分析儀的整體設(shè)計(jì)思路，用軟件代替了硬件，使工程技術(shù)人員可以用一部筆記本電腦到現(xiàn)場就可輕松完成信號(hào)的采集、處理及頻譜分析。在工程領(lǐng)域中， MATLAB是一種倍受程序開發(fā)人員青睞的語言，對于一些需要做大量數(shù)據(jù)運(yùn)算處理的復(fù)雜應(yīng)用以及某些復(fù)雜的頻譜分析算法MATLAB顯得游刃有余。本文將重點(diǎn)介紹基于MATLAB的虛擬頻譜分析儀的設(shè)計(jì)。本文設(shè)計(jì)的虛擬頻譜分析儀

2、的功能包括：(1)音頻信號(hào)信號(hào)輸入。輸入的途徑包括從聲卡輸入、從WAV文件輸入、從信號(hào)發(fā)生器輸入；(2)信號(hào)波形分析。包括幅值、頻率、周期、相位的估計(jì)，并計(jì)算統(tǒng)計(jì)量的峰值、均值、均方值和方差等信息；(3)信號(hào)頻譜分析。頻率、周期的估計(jì)，圖形顯示幅值譜、相位譜和功率譜等信息的曲線。二、實(shí)驗(yàn)原理 2.1時(shí)域抽樣定理時(shí)域抽樣定理給出了連續(xù)信號(hào)抽樣過程中信號(hào)不失真的約束條件：對于基帶信號(hào)，信號(hào)抽樣頻率 大于等于2倍的信號(hào)最高頻率 ，即 。時(shí)域抽樣是把連續(xù)信號(hào) 變成適于數(shù)字系統(tǒng)處理的離散信號(hào) 。對連續(xù)信號(hào) 以間隔T抽樣，則可得到的離散序列為 。 圖1 連續(xù)信號(hào)抽樣的離散序列若 ，則信號(hào) 與 的頻譜之間存

3、在：其中： 的頻譜為 ， 的頻譜為 ?？梢姡盘?hào)時(shí)域抽樣導(dǎo)致信號(hào)頻譜的周期化。 (rad/s)為抽樣角頻率， 為抽樣頻率。數(shù)字角頻率與模擬角頻率的關(guān)系為：=T。2.2 離散傅立葉變換（DFT）有限長序列的離散傅立葉變換（DFT）為逆變換為2.3 快速傅立葉變換（FFT）在各種信號(hào)序列中，有限長序列占重要地位。對有限長序列可以利用離散傅立葉變換(DFT)進(jìn)行分析。DFT不但可以很好的反映序列的頻譜特性，而且易于用快速算法(FFT)在計(jì)算機(jī)上進(jìn)行分析。有限長序列的DFT是其z變換在單位圓上的等距離采樣，或者說是序列傅立葉的等距離采樣，因此可以用于序列的譜分析。FFT是DFT的一種快速算法，它是對變

4、換式進(jìn)行一次次分解，使其成為若干小數(shù)據(jù)點(diǎn)的組合，從而減少運(yùn)算量。 MATLAB為計(jì)算數(shù)據(jù)的離散快速傅立葉變換，提供了一系列豐富的數(shù)學(xué)函數(shù)，主要有Fft、Ifft、Fft2 、Ifft2, Fftn、ifftn和Fftshift、Ifftshift等。當(dāng)所處理的數(shù)據(jù)的長度為2的冪次時(shí)，采用基-2算法進(jìn)行計(jì)算，計(jì)算速度會(huì)顯著增加。所以，要盡可能使所要處理的數(shù)據(jù)長度為2的冪次或者用添零的方式來添補(bǔ)數(shù)據(jù)使之成為2的冪次。Fft函數(shù)調(diào)用方式：Y=fft(X)；Yfft(X,N)；Yfft(X,dim)或Yfft(X,N,dim)。函數(shù)Ifft的參數(shù)應(yīng)用與函數(shù)Fft完全相同。2.4 波形分析原理2.4.1

5、信號(hào)頻率、幅值和相位估計(jì)(1)頻率(周期)檢測對周期信號(hào)來說，可以用時(shí)域波形分析來確定信號(hào)的周期，也就是計(jì)算相鄰的兩個(gè)信號(hào)波峰的時(shí)間差、或過零點(diǎn)的時(shí)間差。這里采用過零點(diǎn)(ti)的時(shí)間差T(周期)。頻率即為f= 1/T，由于能夠求得多個(gè)T值(ti有多個(gè))，故采用它們的平均值作為周期的估計(jì)值。(2)幅值檢測在一個(gè)周期內(nèi)，求出信號(hào)最大值ymax與最小值ymin的差的一半，即A= (ymax-ymin)/2，同樣，也會(huì)求出多個(gè)A值，但第1個(gè)A值對應(yīng)的ymax和ymin不是在一個(gè)周期內(nèi)搜索得到的，故以除第1個(gè)以外的A值的平均作為幅值的估計(jì)值。(3)相位檢測采用過零法，即通過判斷與同頻零相位信號(hào)過零點(diǎn)時(shí)刻

6、，計(jì)算其時(shí)間差，然后換成相應(yīng)的相位差。=2(1-ti/T)，x表示x的小數(shù)部分，同樣，以的平均值作為相位的估計(jì)值。下圖是本設(shè)計(jì)利用過零檢測法估算信號(hào)周期、頻率、幅值和相位的流程圖。圖1 過零檢測法估算信號(hào)的周期、幅值和相位流程圖2.4.2數(shù)字信號(hào)統(tǒng)計(jì)量估計(jì)(1)峰值P的估計(jì)在樣本數(shù)據(jù)x中找出最大值與最小值，其差值為雙峰值，雙峰值的一半即為峰值。P=max(yi)-min(yi)/2(2)均值估計(jì)式中，N為樣本容量，下同。(3)均方值估計(jì)(4)方差估計(jì)2.5頻譜分析原理時(shí)域分析只能反映信號(hào)的幅值隨時(shí)間的變化情況，除單頻率分量的簡單波形外，很難明確提示信號(hào)的頻率組成和各頻率分量大小，而頻譜分析能很

7、好的解決此問題。由于從頻域能獲得的主要是頻率信息，所以本節(jié)主要介紹頻率(周期)的估計(jì)與頻譜圖的生成。2.5.1頻率、周期的估計(jì)對于Y(kf)，如果當(dāng)kf=f時(shí)，Y(kf)取最大值，則f為頻率的估計(jì)值，由于采樣間隔的誤差，f也存在誤差，其誤差最大為f/ 2。周期T=1/f。從原理上可以看出，如果在標(biāo)準(zhǔn)信號(hào)中混有噪聲，用上述方法仍能夠精確地估計(jì)出原標(biāo)準(zhǔn)信號(hào)的頻率和周期，這個(gè)將在下一章做出驗(yàn)證2.5.2頻譜圖為了直觀地表示信號(hào)的頻率特性，工程上常常將Fourier變換的結(jié)果用圖形的方式表示，即頻譜圖。以頻率f為橫坐標(biāo)，|Y(f)|為縱坐標(biāo)，可以得到幅值譜；以頻率f為橫坐標(biāo)，argY(f)為縱坐標(biāo)，可

8、以得到相位譜；以頻率f為橫坐標(biāo)，ReY(f)為縱坐標(biāo)，可以得到實(shí)頻譜；以頻率f為橫坐標(biāo)，ImY(f)為縱坐標(biāo)，可以得到虛頻譜。根據(jù)采樣定理，只有頻率不超過Fs/2的信號(hào)才能被正確采集，即Fourier變換的結(jié)果中頻率大于Fs/2的部分是不正確的部分，故不在頻譜圖中顯示。即橫坐標(biāo)f0,Fs/2三、程序設(shè)計(jì)及編寫采樣頻率Fs與采樣點(diǎn)數(shù)N是聲音信號(hào)輸入時(shí)共同需要作用的參數(shù)，故將其獨(dú)立出來。下面為別介紹三種輸入方式的實(shí)現(xiàn)。3.1 聲卡輸入這里聲卡輸入是指由麥克風(fēng)錄音得到的聲音信號(hào)的輸入，MATLAB提供了wavrecord函數(shù)，該函數(shù)能夠?qū)崿F(xiàn)讀取麥克風(fēng)錄音信號(hào)。以下是“開始錄音”按鈕的回調(diào)函數(shù)內(nèi)容。%

9、首先獲得設(shè)定的Fs值Fs=str2double(get(findobj(Tag,samplerate),String);%根據(jù)設(shè)定的錄音時(shí)長進(jìn)行錄音，將其存入handles.y中handles.y=wavrecord(str2double(get(handles.recordtime,String)*Fs, Fs,int16);%保存handles結(jié)構(gòu)體，使得handles.y在別的函數(shù)中也能使用guidata(hObject,handles);%在波形顯示區(qū)繪出波形plot(handles.time,handles.y);title(WAVE);%將所采到的點(diǎn)的數(shù)量輸出在“采樣點(diǎn)數(shù)”中ysi

10、ze=size(handles.y)set(handles.samplenum,String,num2str(ysize(1);3.2 WAV文件輸入MATLAB提供了wavread函數(shù)，該函數(shù)能夠方便的打開并讀取WAV文件中的聲音信息，并且同時(shí)讀取所有聲道。下面是“打開文件”按鈕回調(diào)函數(shù)的部分代碼。其它代碼與聲卡輸入的類似。%從WAV文件中讀取的聲音信息并臨時(shí)存放到temp變量中temp = wavread(get(findobj(Tag,filename),String);%獲得所選擇的聲道channel=str2double(get(handles.channel,String);%將指

11、定聲道的信息存放到handles.y中handles.y=temp(:,channel);3.3 信號(hào)發(fā)生器輸入MATLAB有產(chǎn)生標(biāo)準(zhǔn)信號(hào)的函數(shù)，如sawtooth能夠產(chǎn)生三角波或鉅齒波，首先利用get函數(shù)獲得波形soundtype，頻率frequency，幅值amp和相位phase，然后是以下代碼。switch soundtypecase 1%標(biāo)準(zhǔn)正弦波y=amp*sin(2*pi*x*frequency+phase);case 2%方波y=amp*sign(sin(2*pi*x*frequency+phase);case 3%三角波y=amp*sawtooth(2*pi*x*frequen

12、cy+phase,0.5);case 4%鉅齒波y=amp*sawtooth(2*pi*x*frequency+phase);case 5%白噪聲y=amp*(2*rand(size(x)-1);otherwiseerrordlg(Illegal wave type,Choose errer);endif get(handles.add,Value)=0.0handles.y=y;%若沒有勾選上“混迭”，則將生成的波形賦給handles.yelse%否則將生成的波形與原有波形疊加handles.y=handles.y+y;end3.4時(shí)域分析2.1.2節(jié)給出時(shí)域分析中的過零檢測算法流程，故這里

13、不給出過零檢測的代碼。MATLAB提供了mean，std函數(shù)，能夠方便地計(jì)算均值、標(biāo)準(zhǔn)差。下面是過零檢測之后的代碼，其中T為過零檢測得到的周期(向量)，amp為過零檢測得到的幅值(向量)，n為過零點(diǎn)數(shù)。freq=Fs/mean(T);%計(jì)算頻率set(handles.outt,String,1/freq);%輸出周期估計(jì)值set(handles.outfreq,String,num2str(freq);%輸出頻率估計(jì)值%計(jì)算并輸出幅值，以幅值均值作為其估計(jì)set(handles.outamp,String,num2str(mean(amp(2:n-1);%將待分析信號(hào)的過零點(diǎn)與標(biāo)準(zhǔn)信號(hào)的過零點(diǎn)

14、相比較，從而得出相位phase=2*pi*(1-(ti(1:n-1)-1)./T+floor(ti(1:n-1)-1)./T);set(handles.outphase,String,num2str(mean(phase);%最大值與最小值的一半即為峰值set(handles.outpeak,String,(max(handles.y(from:to)-min(handles.y(from:to)/2);%from，to即是界面中的“從第from點(diǎn)到第to點(diǎn)”%計(jì)算并輸出均值set(handles.outmean,String,mean(handles.y(from:to);%計(jì)算并輸出均方值

15、set(handles.outmeansquare,String,mean(handles.y(from:to).2);%計(jì)算半輸出方差set(handles.outs,String,std(handles.y(from:to)2);3.5頻域分析頻域分析需要作Fourier變換，MATLAB提供了fft函數(shù)，能夠方便地實(shí)現(xiàn)快速Fourier變換算法。以下代碼省去了從界面中獲得from、to、Fs的部分，也省去了繪圖后設(shè)置橫、縱坐標(biāo)軸的名稱的部分。%首先提取出待分析的樣本，將其存入sample中sample=handles.y(from:to);%生成離散化的頻率點(diǎn)，以采樣頻率作為離散化的間隔

16、f=linspace(0,Fs/2,(to-from+1)/2);%對樣本作快速Fourier變換，變換結(jié)果存入Y中Y=fft(sample,to-from+1);C,I=max(abs(Y);%獲得幅值最大的點(diǎn)及其所對應(yīng)的下標(biāo)值I%則f(I)為最大的幅值所對應(yīng)的頻率，即信號(hào)頻率的估計(jì)值set(handles.foutt,String,1/f(I);%計(jì)算并輸出周期的估計(jì)值set(handles.foutfreq,String,f(I);%輸出頻率的估計(jì)值Y=Y(1:(to-from+1)/2);%為與f對應(yīng)，只取Y的前半部分plot(handles.plot1,f,2*sqrt(Y.*con

17、j(Y);%繪制幅值譜曲線plot(handles.plot2,f,angle(Y);%繪制相位譜曲線plot(handles.plot3,f,real(Y);%繪制實(shí)頻譜曲線plot(handles.plot4,f,imag(Y);%繪制虛頻譜曲線plot(handles.plot5,f,abs(Y).2);%繪制功率譜曲線四、軟件運(yùn)行及結(jié)果分析4.1標(biāo)準(zhǔn)正弦信號(hào)的頻率估計(jì)用信號(hào)發(fā)生器生成標(biāo)準(zhǔn)正弦信號(hào)，然后分別進(jìn)行時(shí)域分析與頻域分析，得到的結(jié)果如圖4所示。從圖中可以看出，時(shí)域分析的結(jié)果為f= 400Hz，頻域分析的結(jié)果為f= 400.2Hz，而標(biāo)準(zhǔn)信號(hào)的頻率為400Hz，從而對于標(biāo)準(zhǔn)信號(hào)時(shí)域分析的精度遠(yuǎn)高于頻域分析的精度。圖3 標(biāo)準(zhǔn)正弦波分析4.2非標(biāo)準(zhǔn)正弦信號(hào)的頻率估計(jì)先成生幅值600的標(biāo)準(zhǔn)正弦信號(hào)，再將幅值300的白噪聲信號(hào)與其混迭，對最終得到的信號(hào)進(jìn)行時(shí)域分析與頻域分析，結(jié)果如圖5所示，可以看出，時(shí)域分析的結(jié)果為f= 388.68Hz，頻域分析的結(jié)果為f= 400.21Hz，而標(biāo)準(zhǔn)信號(hào)的頻率為400Hz，從而對于帶噪聲的正弦信號(hào)頻域分析的精度