現(xiàn)代信號(hào)處理

1、 現(xiàn)代信號(hào)處理課程設(shè)計(jì)實(shí)驗(yàn)報(bào)告 實(shí)驗(yàn)課題：現(xiàn)代信號(hào)處理 專(zhuān)業(yè)班級(jí)： 學(xué)生姓名： 學(xué)生學(xué)號(hào)： 指導(dǎo)老師： 完成時(shí)間： 目錄1 前言-22 課程設(shè)計(jì)內(nèi)容要求及題目-33 設(shè)計(jì)思想和系統(tǒng)功能結(jié)構(gòu)及功能說(shuō)明-44 關(guān)鍵部分的詳細(xì)描述和介紹，流程圖描述關(guān)鍵模塊和設(shè)計(jì)思想-75 問(wèn)題分析及心得體會(huì)-206 參考文獻(xiàn)-217 附錄：程序源代碼清單-21 一、前言 數(shù)字濾波在通信、圖像編碼、語(yǔ)音編碼、雷達(dá)等許多領(lǐng)域中有著十分廣泛的應(yīng)用。目前，數(shù)字信號(hào)濾波器的設(shè)計(jì)在圖像處理、數(shù)據(jù)壓縮等方面的應(yīng)用取得了令人矚目的進(jìn)展和成就。它是數(shù)字信號(hào)處理理論的一部分。數(shù)字信號(hào)處理主要是研究用數(shù)字或符號(hào)的序列來(lái)表示信號(hào)

2、波形，并用數(shù)字的方式去處理這些序列，以便估計(jì)信號(hào)的特征參量，或削弱信號(hào)中的多余分量和增強(qiáng)信號(hào)中的有用分量。具體來(lái)說(shuō)，凡是用數(shù)字方式對(duì)信號(hào)進(jìn)行濾波、變換、調(diào)制、解調(diào)、均衡、增強(qiáng)、壓縮、固定、識(shí)別、產(chǎn)生等加工處理，都可納入數(shù)字信號(hào)處理領(lǐng)域。 數(shù)字信號(hào)處理學(xué)科的一項(xiàng)重大進(jìn)展是關(guān)于數(shù)字濾波器設(shè)計(jì)方法的研究。關(guān)于數(shù)字濾波器，早在上世紀(jì)40年代末期就有人討論設(shè)計(jì)它的可能性問(wèn)題，在50年代也有人討論過(guò)數(shù)字濾波器，但直到60年代中期，才開(kāi)始形成關(guān)于數(shù)字濾波器的一整套完整的正規(guī)理論。在這一時(shí)期，提出了各種各樣的數(shù)字濾波器結(jié)構(gòu)，有的以運(yùn)算誤差最小為特點(diǎn)，有的則以運(yùn)算速度高見(jiàn)長(zhǎng)，而有的則二者兼而有之。出現(xiàn)

3、了數(shù)字濾波器的各種實(shí)現(xiàn)方法，對(duì)遞歸和非遞歸兩類(lèi)濾波器作了全面的比較，統(tǒng)一了數(shù)字濾波器的基本概念和理論。  數(shù)字濾波器與模擬濾波器相比，具有精度高、穩(wěn)定、體積小、重量輕、靈活、不要求阻抗匹配以及能實(shí)現(xiàn)模擬濾波器無(wú)法進(jìn)行的特殊濾波等優(yōu)點(diǎn)。上學(xué)期學(xué)習(xí)了數(shù)字信號(hào)處理這門(mén)課，這學(xué)期的課程設(shè)計(jì)使我更加形象具體的掌握這門(mén)課程，并且可以熟練的運(yùn)用MATLAB進(jìn)行編程，經(jīng)過(guò)此次課程設(shè)計(jì)達(dá)到了其設(shè)計(jì)內(nèi)容的基本要求，實(shí)現(xiàn)了基本功能，并對(duì)程序進(jìn)行了優(yōu)化，以達(dá)到最佳效果。在設(shè)計(jì)前制定一個(gè)設(shè)計(jì)方案，明確設(shè)計(jì)目標(biāo)，方法和過(guò)程，參閱相關(guān)文獻(xiàn)，咨詢老師和同學(xué)，嚴(yán)格按照方案進(jìn)行課程設(shè)計(jì)，最終解決問(wèn)題，完成任務(wù)要求。 二

4、、課程設(shè)計(jì)內(nèi)容要求及題目1、信號(hào)發(fā)生器用戶根據(jù)測(cè)試需要，可任選以下兩種方式之一生成測(cè)試信號(hào)：（1）直接輸入（或從文件讀?。y(cè)試序列；（2）輸入由多個(gè)不同頻率正弦信號(hào)疊加組合而成的模擬信號(hào)公式（如式所示）、采樣頻率（Hz）、采樣點(diǎn)數(shù)，動(dòng)態(tài)生成該信號(hào)的采樣序列，作為測(cè)試信號(hào)。 px=100*sin(2*pi*15*t)+100*sin(2*pi*20*t)+100*sin(2*pi*30*t) 2、頻譜分析 使用FFT對(duì)產(chǎn)生的測(cè)試信號(hào)進(jìn)行頻譜分析并展示其幅頻特性與相頻特性，指定需要濾除的頻帶，通過(guò)選擇濾波器類(lèi)型（IIR/FIR），確定對(duì)應(yīng)的濾波器（低通、高通）技術(shù)指標(biāo)。3、濾波器設(shè)計(jì) 根據(jù)以上技術(shù)

5、指標(biāo)（通帶截止頻率、通帶最大衰減、阻帶截止頻率、阻帶最小衰減），設(shè)計(jì)數(shù)字濾波器，生成相應(yīng)的濾波器系數(shù)，并畫(huà)出對(duì)應(yīng)的濾波器幅頻特性與相頻特性。（1）IIR DF設(shè)計(jì)：可選擇濾波器基型（巴特沃斯或切比雪夫型）；（2）FIR DF設(shè)計(jì)：使用窗口法（可選擇窗口類(lèi)型，并比較分析基于不同窗口、不同階數(shù)所設(shè)計(jì)數(shù)字濾波器的特點(diǎn)）。4、數(shù)字濾波 根據(jù)設(shè)計(jì)的濾波器系數(shù)，對(duì)測(cè)試信號(hào)進(jìn)行數(shù)字濾波，展示濾波后信號(hào)的幅頻特性與相頻特性，分析是否滿足濾波要求（對(duì)同一濾波要求，對(duì)比分析各類(lèi)濾波器的差異）。（1）IIRDF：要求通過(guò)差分方程迭代實(shí)現(xiàn)濾波（未知初值置零處理）；（2）FIRDF：要求通過(guò)快速卷積實(shí)現(xiàn)濾波（對(duì)于長(zhǎng)序列

6、，可以選擇使用重疊相加或重疊保留法進(jìn)行卷積運(yùn)算）。5、選做內(nèi)容將一段語(yǔ)音作為測(cè)試信號(hào)，通過(guò)頻譜展示和語(yǔ)音播放，對(duì)比分析濾波前后語(yǔ)音信號(hào)的變化，進(jìn)一步加深對(duì)數(shù)字信號(hào)處理的理解。三、設(shè)計(jì)思想和系統(tǒng)功能結(jié)構(gòu)及功能說(shuō)明1、設(shè)計(jì)思想 首先設(shè)定由多個(gè)不同頻率正弦信號(hào)疊加組合而成的模擬信號(hào)函數(shù)、決定采樣頻率（Hz）和采樣點(diǎn)數(shù)，然后使用FFT對(duì)產(chǎn)生的測(cè)試信號(hào)進(jìn)行頻譜分析并展示其幅頻特性與相頻特性，指定需要濾除的頻帶，接下來(lái)使用漢寧窗低通濾波器與巴特沃斯高通濾波器進(jìn)行濾波，分析濾波后的信號(hào)的幅頻特性以及相頻特性，最后設(shè)計(jì)語(yǔ)音濾波，檢查是否滿足濾波要求。2、系統(tǒng)結(jié)構(gòu)功能說(shuō)明（1）信號(hào)發(fā)生器 定義一個(gè)由4個(gè)不同頻率

7、疊加的正弦函數(shù)，輸出該函數(shù)，觀察其相頻以及幅頻特性。（2）頻譜分析 使用FFT對(duì)產(chǎn)生的測(cè)試信號(hào)進(jìn)行頻譜分析并展示其幅頻特性與相頻特性，指定需要濾除的頻帶，通過(guò)選擇濾波器類(lèi)型漢寧窗和巴特沃斯濾波器，確定對(duì)應(yīng)的濾波器（低通、高通）技術(shù)指標(biāo)Fp、Fc、Rp、Rs。然后對(duì)其濾波后信號(hào)進(jìn)行頻譜分析。（3）數(shù)字濾波器設(shè)計(jì)1、FIR 數(shù)字濾波器設(shè)計(jì)原理 FIR 數(shù)字濾波器可以用下面的差分方程來(lái)描述其中x（n）是輸入序列，y（n）是輸出序列，h（n）是系統(tǒng)的單位脈沖響應(yīng)。由于其采用卷積的方法得到輸出，故可以采用快速算法FFT進(jìn)行卷積計(jì)算，提高信號(hào)處理速度。FIR 數(shù)字濾波器滿足線性相位的條件是其乘法器的系數(shù)關(guān)

8、于中點(diǎn)對(duì)稱(chēng)，又因?yàn)槠浣Y(jié)構(gòu)無(wú)反饋的回路，所以是無(wú)條件的穩(wěn)定系統(tǒng)。2、 窗函數(shù)法設(shè)計(jì)及MATLAB 實(shí)現(xiàn)FIR 濾波器目前常用的設(shè)計(jì)方法有窗函數(shù)法和頻率采樣法，窗函數(shù)法是從時(shí)域進(jìn)行設(shè)計(jì)，而頻率采樣法是從頻域進(jìn)行設(shè)計(jì)。窗函數(shù)法由于簡(jiǎn)單、物理意義清晰，因而得到了較為廣泛的應(yīng)用。窗函數(shù)法設(shè)計(jì)的基本思想是：首先根據(jù)技術(shù)指標(biāo)要求，選取合適的階數(shù)N和窗函數(shù)的類(lèi)型，使其幅頻特性逼近理想濾波器幅頻特性。其次，因?yàn)槔硐霝V波器的是無(wú)限長(zhǎng)的，所以需要對(duì)進(jìn)行截?cái)?，?shù)學(xué)上稱(chēng)這種方法為窗函數(shù)法1。其中是希望逼近的濾波器的單位脈沖響應(yīng)，是希望逼近的濾波器的幅頻特性，是所設(shè)計(jì)的濾波器的單位脈沖響應(yīng)，是窗函數(shù)。在截?cái)嗟倪^(guò)程中要保證

9、FIR濾波器的穩(wěn)定性與線性相位的特性。由于截?cái)嘈?yīng)使所設(shè)計(jì)濾波器的幅頻特性存在誤差，稱(chēng)之為吉布斯效應(yīng)。最后，驗(yàn)證所設(shè)計(jì)的濾波器是否滿足要求，若不滿足，重新設(shè)計(jì)。常用的窗函數(shù)有矩形窗、三角窗、漢寧窗和凱塞窗等。窗函數(shù)的選擇一般要滿足以下三個(gè)條件：（1）具有較低的旁瓣幅度，尤其是第一旁瓣幅度。（2）旁瓣幅度下降速度要快，以利于增加阻帶衰減。（3）主瓣寬度要窄，以獲得較陡的過(guò)渡帶。FIR 濾波器的設(shè)計(jì)任務(wù)是選擇有限長(zhǎng)度的 h(n), 使函數(shù)H(ej)滿足技術(shù)要求。我選擇的是漢寧窗低通濾波。3、 IIR數(shù)字濾波器的設(shè)計(jì)IIR數(shù)字濾波器沖激響應(yīng)無(wú)限寬，與模擬濾波器相匹配。因此可采取在模擬濾波器設(shè)計(jì)的基礎(chǔ)

10、上進(jìn)一步變換的方法設(shè)計(jì)IIR濾波器，即通過(guò)已知的模擬濾波器系統(tǒng)的系統(tǒng)函數(shù)H(s)來(lái)設(shè)計(jì)數(shù)字濾波器的系統(tǒng)函數(shù)H(z)，主要是通過(guò)脈沖響應(yīng)不變法或雙線性變換法完成s平面到Z平面的變換。通過(guò)典型的模擬濾波器（如巴特沃斯濾波器、切比雪夫?yàn)V波器等）可實(shí)現(xiàn)一定參數(shù)要求的數(shù)字濾波器。IIR數(shù)字濾波器的設(shè)計(jì)具體步驟：<1> 按照既定規(guī)則將數(shù)字濾波器技術(shù)指標(biāo)轉(zhuǎn)為模擬低通濾波器技術(shù)指標(biāo)。 <2> 以轉(zhuǎn)換后的技術(shù)指標(biāo)為依據(jù)設(shè)計(jì)模擬低通濾波器H(s)（傳遞函數(shù)) 。<3> 再按照一定規(guī)則將H(s)轉(zhuǎn)換成H(z)（數(shù)字濾波器的傳遞函數(shù)）。若設(shè)計(jì)的是高通、帶通或帶阻濾波器，還需進(jìn)行以下

11、步驟：將高通、帶通或帶阻數(shù)字濾波器的技術(shù)指標(biāo)轉(zhuǎn)換為低通模擬濾波器的技術(shù)指標(biāo), 然后設(shè)計(jì)出低通H(s), 再將H(s)轉(zhuǎn)換為H(z)。我選擇的是巴特沃斯高通濾波器。（4）語(yǔ)音濾波1、語(yǔ)音信號(hào)的采集 MATLAB軟件平臺(tái)下，利用函數(shù)wavread對(duì)語(yǔ)音信號(hào)進(jìn)行采樣，記住采樣頻率和采樣點(diǎn)數(shù)。通過(guò)wavread函數(shù)的使用，我們很快理解了采樣頻率、采樣位數(shù)等概念。 這里我直接采用了一段現(xiàn)成的wav格式的語(yǔ)音信號(hào)。2、語(yǔ)音信號(hào)的頻譜分析 首先畫(huà)出語(yǔ)音信號(hào)的時(shí)域波形；然后對(duì)語(yǔ)音號(hào)進(jìn)行快速傅里葉變換，得到信號(hào)的頻譜特性。四、關(guān)鍵部分的詳細(xì)描述和介紹，流程圖描述關(guān)鍵模塊和設(shè)計(jì)思想1、信號(hào)發(fā)生器（1）輸入信號(hào)程序

12、如下：N=512;T=0.01;Fs=1/T;t=0:pi/256:2*pi;f1=7;f2=37;f3=39;f4=43;x=100*sin(2*pi*f1*t)+100*sin(2*pi*f2*t)+100*sin(2*pi*f3*t)+100*sin(2*pi*f4*t);plot(t,x,'g');grid;title('原始輸入信號(hào)');xlabel('時(shí)間/t');ylabel('振幅');X=fft(x,N);plot(X);grid;title('輸入信號(hào)頻譜');設(shè)計(jì)輸入信號(hào)并對(duì)其進(jìn)行傅里葉變換。

13、輸入信號(hào)如下圖（2）輸入信號(hào)頻譜分析幅頻特性Fs=100;len=512;f=Fs*(0:len/2-1)/len;plot(f,abs(X(1:len/2);gridtitle('原始信號(hào)幅頻特性曲線');相頻特性Fs=100;len=512;f=Fs*(0:len/2-1)/len;plot(f,angle(X(1:len/2);grid;title('相頻特性曲線');幅頻特性相頻特性2、巴特沃斯高通濾波器 高通濾波器采用Butterworth濾波器, 因?yàn)樗哂休^平坦的幅頻特性, 而且其元件值比較合乎實(shí)際、對(duì)精度要求稍寬。用N,Wn=buttord(Wp

14、,Ws,Rp,Rs,'s')語(yǔ)句實(shí)現(xiàn)Butterworth濾波器階數(shù)和邊緣頻率的確定；b,a=butter(N,Wn,'s')語(yǔ)句來(lái)實(shí)現(xiàn)Butterworth濾波器的設(shè)計(jì)。其中Wp為通帶截止頻率；Rp為通帶截止頻率處的衰減；Ws為阻帶截止頻率；Rs為阻帶截止頻率處的衰減； b,a是濾波器傳遞函數(shù)多項(xiàng)式的系數(shù)；n為濾波器的階數(shù)；Wn為截止頻率。該通帶截止頻率為14HZ。（1）程序如下：T=0.01;Fs=1/T;wp=0.4*pi;ws=0.3*pi;Rp=1;Rs=15;Nn=512;%采樣間隔Wp=wp/T;Ws=ws/T;N,Wn=buttord(Wp,Ws

15、,Rp,Rs,'s');%模擬濾波器階數(shù)z,p,k=buttap(N);%數(shù)字濾波器設(shè)計(jì)Bap,Aap=zp2tf(z,p,k);b,a=lp2hp(Bap,Aap,Wn);bz,az=bilinear(b,a,1/T);%脈沖響應(yīng)不變法figurefreqz(bz,az,Nn);H,F=freqz(bz,az,Nn,1/T);plot(F,20*log10(abs(H),'r');grid;title('巴特沃斯濾波器幅頻特性');xlabel('頻率/HZ');ylabel('振幅/dB');巴特沃斯幅頻特性（

16、2） 濾波后信號(hào)巴特沃斯高通濾波器濾除了7HZ的信號(hào)，保留了高于14HZ的信號(hào)。濾波后信號(hào)幅頻特性相頻特性3、漢寧窗低通濾波器通帶截止頻率為16HZ，濾除16HZ以上的信號(hào)，只保留7HZ的信號(hào)。程序如下：wp=0.3*pi;ws=0.45*pi;wdelta=ws-wp;N=ceil(8*pi/wdelta);if rem(N,2)=0 N=N+1;end;Nw=N;wc=(wp+ws)/2;n=0:N-1;alpha=(N-1)/2;%相位延遲m=n-alpha+0.00001;%0.00001為MATALB系統(tǒng)的精度hd=sin(wc*m)./(pi*m);%理性濾波器脈沖響應(yīng)win=ha

17、nn(Nw);h=hd.*win'%時(shí)域乘積對(duì)應(yīng)頻域卷積b=h;H,f=freqz(b,1,512,100);%100HZ的采樣頻率plot(f,20*log10(abs(H),'r');grid;xlabel('頻率/HZ');ylabel('振幅/dB');title('漢寧窗幅頻特性曲線');漢寧窗幅頻特性濾波后信號(hào)幅頻特性相頻特性4、 語(yǔ)音濾波器利用函數(shù)wavread對(duì)語(yǔ)音信號(hào)進(jìn)行采樣，利用巴特沃斯濾波器進(jìn)行濾波。程序如下：fs=32768;x1=wavread('C:UsersloryDesktopwo.

18、wav');t=0:1/32768:(size(x1)-1)/32768;Au=0.03;d=Au*cos(2*pi*10000*t)'x2=x1+d;wp=0.25*pi;ws=0.3*pi;Rp=1;Rs=15;Fs=32768;Ts=1/Fs;wp1=2/Ts*tan(wp/2); %將模擬指標(biāo)轉(zhuǎn)換成數(shù)字指標(biāo)ws1=2/Ts*tan(ws/2); N,Wn=buttord(wp1,ws1,Rp,Rs,'s');%選擇濾波器的最小階數(shù)Z,P,K=buttap(N); %創(chuàng)建butterworth模擬濾波器Bap,Aap=zp2tf(Z,P,K);b,a=lp

19、2lp(Bap,Aap,Wn); bz,az=bilinear(b,a,Fs); %用雙線性變換法實(shí)現(xiàn)模擬濾波器到數(shù)字濾波器的轉(zhuǎn)換H,W=freqz(bz,az); %繪制頻率響應(yīng)曲線figure(1)plot(W*Fs/(2*pi),abs(H)gridxlabel('頻率Hz')ylabel('頻率響應(yīng)幅度')title('Butterworth')f1=filter(bz,az,x2);figure(2)subplot(2,1,1)plot(t,x2) %畫(huà)出濾波前的時(shí)域圖title('濾波前的時(shí)域波形');subplot(

20、2,1,2)plot(t,f1); %畫(huà)出濾波后的時(shí)域圖title('濾波后的時(shí)域波形');sound(f1,32768); %播放濾波后的信號(hào)F0=fft(f1,1024);f=fs*(0:511)/1024;figure(3)y2=fft(x2,1024);subplot(2,1,1);plot(f,abs(y2(1:512); %畫(huà)出濾波前的頻譜圖title('濾波前的頻譜')xlabel('Hz');ylabel('fuzhi');subplot(2,1,2)F1=plot(f,abs(F0(1:512); %畫(huà)出濾波后的

21、頻譜圖title('濾波后的頻譜')xlabel('Hz');ylabel('fuzhi');巴特沃斯幅頻特性 五、問(wèn)題分析及心得體會(huì) 經(jīng)過(guò)這次課程設(shè)計(jì)，我的收獲良多，也遇到了很多問(wèn)題。在這之前我們并沒(méi)有接觸過(guò)MATLAB，所以在課程設(shè)計(jì)時(shí)查閱了大量有關(guān)MATLAB的書(shū)籍以及與課程相關(guān)的設(shè)計(jì)案例。但是編程過(guò)程中總是出現(xiàn)一些錯(cuò)誤，變量的定義，參數(shù)的傳遞，函數(shù)的調(diào)用等各方面大大小小的問(wèn)題層出不窮，為此，只有仔細(xì)認(rèn)真、不急不躁才能解決這些小問(wèn)題，這次課程設(shè)計(jì)我分成了四個(gè)模塊進(jìn)行編程，信號(hào)發(fā)生器，巴特沃斯高通濾波器，漢寧窗低通濾波器，還有語(yǔ)音濾波部分。在

22、巴特沃斯沃斯濾波時(shí)由于參數(shù)取的不好，有濾波殘留，改一下參數(shù)，濾波效果很好。開(kāi)始時(shí)幅頻特性與相頻特性是對(duì)稱(chēng)分布的，通過(guò)加入幾行代碼，讓它只顯示一半。做語(yǔ)音濾波時(shí)，讀不到我的語(yǔ)音文件，后來(lái)通過(guò)音頻轉(zhuǎn)換器把它的格式轉(zhuǎn)換成了wav格式的文件后，就可以正常運(yùn)行了。最后，在做GUI界面時(shí)，因?yàn)樗枰龅膱D太多，有很多個(gè)按鈕，后來(lái)我把各個(gè)模塊分成了四部分，再建立一個(gè)總的調(diào)用函數(shù)來(lái)調(diào)出各個(gè)模塊。形成了主界面與子界面的模式。畫(huà)面看著直觀，清晰。這次的課程設(shè)計(jì)，使我對(duì)MATLAB軟件運(yùn)用能力的有所提高，并且將所學(xué)的數(shù)字信號(hào)處理知識(shí)應(yīng)用在課程設(shè)計(jì)之中，提升了對(duì)數(shù)字信號(hào)處理技術(shù)的掌握。參數(shù)的使用，函數(shù)的調(diào)用，信號(hào)的產(chǎn)

23、生，頻譜分析，濾波器參數(shù)的設(shè)定及不同濾波器的設(shè)計(jì)都要注意，避免細(xì)節(jié)錯(cuò)誤。經(jīng)過(guò)這次課程設(shè)計(jì)，提高了我對(duì)問(wèn)題的處理能力，有了問(wèn)題，在經(jīng)過(guò)不斷的努力和老師與同學(xué)的幫助下成功的解決了。我不僅在數(shù)字信號(hào)處理知識(shí)的理解、設(shè)計(jì)思想的構(gòu)建有了進(jìn)一步的提升，也對(duì)MATLAB數(shù)學(xué)處理軟件有了一定的認(rèn)識(shí)和掌握，并認(rèn)識(shí)到編程習(xí)慣的重要性。鍛煉了我的動(dòng)手能力，豐富了我的思維模式。我們要盡最大努力完成任務(wù)。 六、參考文獻(xiàn)1 S.K.Mitra.Digital Signal Processing: A Computer Based Approach, 3rdEditionM,New York, USA:McGraw-Hil

24、l,20002 R.G.Lyons.Understanding Digital Signal Processing,2nd EditionM.New Jersey, USA:Prentice Hall,20053 程佩青.數(shù)字信號(hào)處理教程,第二版M.北京:清華大學(xué)出版社,20014 趙樹(shù)杰等.數(shù)字信號(hào)處理M.西安:西安電子科技大學(xué)出版社,19975 丁玉美等.數(shù)字信號(hào)處理時(shí)域離散隨機(jī)信號(hào)處理M.西安:西安電子科技大學(xué)出版社,20026 陳懷琛等.MATLAB及在電子信息課程中的應(yīng)用M,北京:電子工業(yè)出版社出版,2002 七、程序源代碼清單N=512;T=0.01;Fs=1/T;t=0:pi/2

25、56:2*pi;f1=7;f2=37;f3=39;f4=43;x=100*sin(2*pi*f1*t)+100*sin(2*pi*f2*t)+100*sin(2*pi*f3*t)+100*sin(2*pi*f4*t);figure;subplot(411)plot(t,x,'g');grid;title('原始輸入信號(hào)');xlabel('時(shí)間/t');ylabel('振幅');X=fft(x,N);subplot(412)plot(X);grid;title('輸入信號(hào)頻譜');Fs=100;len=512;f=

26、Fs*(0:len/2-1)/len;subplot(413)plot(f,angle(X(1:len/2);grid;title('相頻特性曲線');Fs=100;len=512;f=Fs*(0:len/2-1)/len;subplot(414)plot(f,abs(X(1:len/2);gridtitle('原始信號(hào)幅頻特性曲線');T=0.01;Fs=1/T;wp=0.4*pi;ws=0.3*pi;Rp=1;Rs=15;Nn=512;%采樣間隔Wp=wp/T;Ws=ws/T;N,Wn=buttord(Wp,Ws,Rp,Rs,'s');%模擬

27、濾波器階數(shù)z,p,k=buttap(N);%數(shù)字濾波器設(shè)計(jì)Bap,Aap=zp2tf(z,p,k);b,a=lp2hp(Bap,Aap,Wn);bz,az=bilinear(b,a,1/T);%脈沖響應(yīng)不變法figurefreqz(bz,az,Nn);H,F=freqz(bz,az,Nn,1/T);figure;subplot(411)plot(F,20*log10(abs(H),'r');grid;title('巴特沃斯濾波器幅頻特性');xlabel('頻率/HZ');ylabel('振幅/dB');N=512;T=0.01;

28、t=0:pi/256:2*pi;f1=7;f2=37;f3=39;f4=43;x=100*sin(2*pi*f1*t)+100*sin(2*pi*f2*t)+100*sin(2*pi*f3*t)+100*sin(2*pi*f4*t);y2=filtfilt(bz,az,x);subplot(412)plot(t,y2,'r');title('濾波后信號(hào)');xlabel('時(shí)間/s');grid;legend('filtfilt');Fs=100;len=200;Y1=fft(y2,len);f=Fs*(0:len/2-1)/le

29、n;subplot(413)plot(f,abs(Y1(1:len/2);grid;title('濾波后信號(hào)幅頻特性曲線');Fs=100;len=512;f=Fs*(0:len/2-1)/len;subplot(414)plot(f,abs(X(1:len/2);gridtitle('原始信號(hào)幅頻特性曲線');wp=0.3*pi;ws=0.45*pi;wdelta=ws-wp;N=ceil(8*pi/wdelta);if rem(N,2)=0 N=N+1;end;Nw=N;wc=(wp+ws)/2;n=0:N-1;alpha=(N-1)/2;%相位延遲m=n-

30、alpha+0.00001;%0.00001為MATALB系統(tǒng)的精度hd=sin(wc*m)./(pi*m);%理性濾波器脈沖響應(yīng)win=hann(Nw);h=hd.*win'%時(shí)域乘積對(duì)應(yīng)頻域卷積b=h;H,f=freqz(b,1,512,100);%100HZ的采樣頻率figure;subplot(411)plot(f,20*log10(abs(H),'r');grid;xlabel('頻率/HZ');ylabel('振幅/dB');title('漢寧窗幅頻特性曲線');N=512;T=0.01;t=0:pi/256:

31、2*pi;f1=7;f2=37;f3=39;f4=43;x=100*sin(2*pi*f1*t)+100*sin(2*pi*f2*t)+100*sin(2*pi*f3*t)+100*sin(2*pi*f4*t);y4=filtfilt(b,1,x);subplot(412)plot(t,y4,'r');grid;title('FIR濾波后輸出信號(hào)')legend('filtfilt');xlabel('時(shí)間/t');ylabel('振幅');Fs=100;len=200;Y2=fft(y4,len);f=Fs*(0:len/2-1)/len;subplot(413)plot(f,abs(Y2(1:len/2);grid;title('濾波后信號(hào)幅頻特性曲線');global Y2subplot(414)plot(angle(Y2);grid