數(shù)字信號處理論文

1、數(shù)字信號處理實驗報告班級：1402014姓名：陳曉霞學(xué)號系方式安電子科技大學(xué)電子工程學(xué)院摘要本試驗報告按照數(shù)字信號處理課程實驗2017的內(nèi)容及要求而進行的MATLAB仿真實驗。共有三個實驗，每個實驗又有若干實驗要求。其中第一個實驗主要是實現(xiàn)常見序列的產(chǎn)生，序列運算的應(yīng)用，以及復(fù)雜序列的產(chǎn)生和序列的線性卷積卷積；實驗二主要是通過編程實現(xiàn)DFT ,IDFT,DIT-FFT,DIF-FFT,IFFT并利用DFT對連續(xù)信號分析；實驗三則要求設(shè)計IIR以及FIR低通，帶通，高通濾波器并對所給音頻信號進行相應(yīng)處理。報告主要包含解決上述內(nèi)容的程序源代碼，m

2、atlab處理圖像，以及實驗結(jié)論三部分。同時在報告所處的文件夾中含有相應(yīng).m文件以及處理后的音頻信號。緒論 數(shù)字信號處理起源于十八世紀的數(shù)學(xué)，隨著信息科學(xué)和計算機技術(shù)的迅速發(fā)展，數(shù)字信號處理的理論與應(yīng)用得到迅速發(fā)展，形成一門極其重要的學(xué)科。當今數(shù)字信號處理的理論和方法已經(jīng)得到長足的發(fā)展，成為數(shù)字化時代的重要支撐，其在各個學(xué)科和技術(shù)領(lǐng)域中的應(yīng)用具有悠久的歷史，已經(jīng)滲透到我們生活和工作的各個方面。 數(shù)字信號處理相對于模擬信號處理具有許多優(yōu)點，比如靈活性好，數(shù)字信號處理系統(tǒng)的性能取決于系統(tǒng)參數(shù)，這些參數(shù)很容易修改，并且數(shù)字系統(tǒng)可以分時復(fù)用，用一套數(shù)字系統(tǒng)可以分是處理多路信號；高精度和高穩(wěn)定性，數(shù)字系

3、統(tǒng)的運算字符有足夠高的精度，同時數(shù)字系統(tǒng)不會隨使用環(huán)境的變化而變化，尤其使用了超大規(guī)模集成的DSP芯片，簡化了設(shè)備，更提高了系統(tǒng)穩(wěn)定性和可靠性；便于開發(fā)和升級，由于軟件可以方便傳送，復(fù)制和升級，系統(tǒng)的性能可以得到不斷地改善；功能強，數(shù)字信號處理不僅能夠完成一維信號的處理，還可以試下安多維信號的處理；便于大規(guī)模集成，數(shù)字部件具有高度的規(guī)范性，對電路參數(shù)要求不嚴格，容易大規(guī)模集成和生產(chǎn)。 從20世紀90年代末以來，數(shù)字信號處理課程幾乎在國內(nèi)所有大學(xué)的電氣信息類等學(xué)科專業(yè)的本科生和研究生中開設(shè)，且是本科生的必修課和研究生的學(xué)位課。數(shù)字信號處理課程在經(jīng)歷了不斷豐富發(fā)展的過程 后日臻完善。課程體系主要有

4、信號分析與處理，離散系統(tǒng)設(shè)計構(gòu)成，同時增加了近代信號處理的理論和方法，并將MATLAB與數(shù)字信號處理有機結(jié)合，作為信號處理的仿真分析手段，從而將理論分析與計算機仿真融為一體。 信號處理在生物醫(yī)學(xué)工程，地震學(xué)，聲吶，雷達，通信，控制等領(lǐng)域都日益顯示其重要作用。例如在醫(yī)學(xué)信號或地震信號分析中，我們需要提取某些重要的特征參數(shù)，在雷達和通信信號處理中，我們希望提出信號中的噪聲或干擾。 數(shù)字信號處理用途廣泛，對其進行一系列學(xué)習(xí)與研究也是非常必要的。同時數(shù)字信號處理技術(shù)正在不斷地克服自身的一些缺點，以驚人的速度發(fā)展，不斷擴大應(yīng)用領(lǐng)域，毫無疑問其前景非常好。實驗一實驗?zāi)康募由顚π蛄谢局R的掌握理解實驗原理

5、與方法1.幾種常見的典型序列：2.序列運算的應(yīng)用：數(shù)字信號處理中經(jīng)常需要將被加性噪聲污染的信號中移除噪聲，假定信號 s(n)被噪聲d(n)所污染，得到了一個含噪聲的信號。我們的目的是對 x(n)進行運算，產(chǎn)生一個合理逼近 s(n)的信號 y(n)，從而實現(xiàn)去噪的目的。通常做法是對時刻 n 的樣本附近的一些 樣本進行求平均，產(chǎn)生輸出信號，這是一種簡單有效的辦法，例如三點滑動平均算法的表達式如下：。3.復(fù)雜信號的產(chǎn)生：調(diào)幅廣播（AM）是一種重要廣播形式，使用的是振幅調(diào)制信號，其產(chǎn)生可用低頻調(diào)制信號（如聲音信號）來調(diào)制高頻正弦信號（如載 波信號），不妨假設(shè)低頻調(diào)制信號和高頻正弦信號均為正弦信號形式，

6、具體如下： ，振幅調(diào)制信號為： 4.線性卷積：實驗內(nèi)容及步驟1.序列的產(chǎn)生：根據(jù)所學(xué)的知識產(chǎn)生單位階躍序列，復(fù)指數(shù)序列，實指數(shù)序列，正弦序列，隨機信號序列。 2.序列運算：設(shè)計一個混有白噪聲的正弦信號，并通過三點滑動平均算法實現(xiàn)去噪，并通過時域圖對比原正弦信號與去噪后獲得的結(jié)果，觀察去噪效果。3.復(fù)雜信號的產(chǎn)生：AM1323kHz是陜西交通廣播，張旭東正在發(fā)出一個頻率為 800Hz 的聲音，試以上述信息為基礎(chǔ)，利用matlab 工具產(chǎn)生此時交通廣播臺此時發(fā)出的廣播信號。4.序列線性卷積驗證實驗：已知：試利用對位相乘法或者你熟悉的方法手動計算 y(n)，并通過matlab對手動計算結(jié)果進行驗證，

7、畫圖說明。實驗結(jié)果分析及結(jié)論總結(jié) 1序列的產(chǎn)生：a) 單位階躍序列：n0=0;n1=-10;n2=10;n=n1:n2; %建立序列自變量向量x=(n-n0)>=0; %建立單位階躍序列向量stem(n,x); %繪制單位階躍序列曲線axis(-10 10 -1 1.5); %固定x與y軸刻度最大值及最小值xlabel('n');ylabel('x(n)'); %標注x與y軸代表參數(shù)title('單位階躍序列') %添加曲線名稱 圖1.1 單位階躍序列 b) 復(fù)指數(shù)序列：n=-10:1:10;A=1;alpha=-0.1+0.3j; x=A

8、*exp(alpha*n); %建立復(fù)指數(shù)序列向量 subplot(221);stem(n,real(x);title('實部');xlabel('n') subplot(222);stem(n,imag(x);title('虛部');xlabel('n')subplot(223);stem(n,abs(x);title('振幅 ');xlabel('n')subplot(224);stem(n,x);title('復(fù)指數(shù)序列');xlabel('n') %依次繪制復(fù)

9、指數(shù)序列的實部，虛部，振幅以及自身曲線圖1.2 復(fù)指數(shù)序列c) 實指數(shù)序列：clear all;close all;clc;n=0:30;a1=0.8;a2=-0.8;a3=1.5;a4=-1.5; %取不同底數(shù)x1=a1.n;x2=a2.n;x3=a3.n;x4=a4.n; %建立實指數(shù)序列subplot(221);stem(n,x1);title('x=0.8n');xlabel('n') %繪制曲線subplot(222);stem(n,x2);title('x=(-0.8)n');xlabel('n')subplot(22

10、3);stem(n,x3);title('x=1.5n');xlabel('n')subplot(224);stem(n,x4);title('x=(-1.5)n');xla 圖1.3 實指數(shù)序列 d) 正弦序列：n=-50:50;A=5;B=0;w=2*pi*0.01;x=A*sin(w*n+B); %建立正弦序列stem(n,x); %繪制曲線xlabel('n');ylabel('x(n)');title('正弦序列')圖1.4 正弦序列e) 隨機信號序列：n=1:20;A=1;x=A*ran

11、d(1,20); %建立隨機序列stem(n,x); %繪制曲線xlabel('n');ylabel('x(n)');title('隨機序列') 圖1.5 隨機信號序列2.序列運算：clear all；close all；clcR=51; d=2*(rand(1,R)-0.5); %建立加性噪聲向量n=0:R-1; s=sin(pi/6)*n); %建立原始離散信號x=s+d; %建立被噪聲影響后的信號subplot(2,1,1); plot(n,d,'r-',n,s,'g-',n,x,'b-.')

12、; %繪制加性噪聲曲線xlabel('Time index n');ylabel('Amplitude');legend('dn','sn','xn'); %添加圖標d添加標注x1=0 0 x;x2=0 x 0;x3=x 0 0;y=(x1+x2+x3)/3; %三點滑動平均算法subplot(2,1,2);plot(n,y(2:R+1),'r-',n,s,'g-'); %繪制濾出噪聲后的信號legend('yn','sn');xlabel('

13、Time index n');ylabel('Amplitude');圖1.6 除噪聲前后信號 3.復(fù)雜信號的產(chǎn)生：clear allclose allclcR=30;n=1:R-1;A=2;m=0.5;x_l=cos(2*pi*0.8*n); %建立調(diào)制信號x_h=cos(2*pi*1323*n); %建立高頻載波信號y(n)=A*(1+m.*x_l).*x_h; %建立調(diào)幅信號stem(n,y(n) %繪制調(diào)幅曲線title('調(diào)幅廣播')圖1.7 調(diào)幅廣播4.序列線性卷積：clear all;N=5;M=6;L=M+N-1;x=1,2,3,4,5;

14、h=6,2,3,6,4,2; %建立卷積序列y=conv(x,h); %進行線性卷積運算nx=0:N-1;nh=0:M-1;ny=0:L-1; subplot(131); stem(nx,x);xlabel('n');ylabel('x(n)');grid on;title('x(n)');subplot(132);stem(nh,h);xlabel('n');ylabel('h(n)');grid on;title('h(n)');subplot(133);stem(ny,y);xlabel(&#

15、39;n');ylabel('y(n)');grid on;title('y(n)=x(n)*h(n)');圖1.8 線性卷積驗證結(jié)果：計算結(jié)果與對位相乘法計算所得結(jié)果一致。 思考題1) 可以通過什么參數(shù)控制序列的增長率與衰減率？哪個參數(shù)可以控制序列的幅值？圖1.9 復(fù)指數(shù)序列振幅控制因素圖1.10 復(fù)指數(shù)序列變化率控制因素通過MATLAB繪制圖形結(jié)果得出：對于復(fù)指數(shù)序列，可以控制序列的幅值，控制序列的增長率與衰減率。2) 復(fù)指數(shù)序列的實部和虛部分別是什么？圖1.11 復(fù)指數(shù)序列的實部與虛部特點 通過MATLAB繪制圖形結(jié)果得出： 當時，復(fù)指數(shù)序列的實部

16、和虛部分別是按指數(shù)規(guī)律增長的正弦振蕩序列；當時，復(fù)指數(shù)序列的實部和虛部分別是按指數(shù)規(guī)律衰減的正弦振蕩序列；當時，復(fù)指數(shù)序列即為虛指數(shù)序列，其實部和虛部分別是按照等幅的正弦振蕩序列。3) 正弦序列的頻率是如何控制的？如何可以改變其頻率？圖1.12 正弦序列頻率控制因素 通過MATLAB繪制圖形結(jié)果可以得出： 正弦序列的頻率是由參數(shù)控制的，并且可通過修改的取值來改變正弦序列的頻率。4) 對于一個幅值為5v,頻率為20Hz,初始相位為60度正弦連續(xù)信號，如何實現(xiàn)對其不失真采樣？ 圖1.13 正弦連續(xù)信號采樣 通過MATLAB繪制圖形結(jié)果可以看出：實現(xiàn)采樣不失真的條件是采樣周期滿足：。（以上代碼見qu

17、es1.m,ques12.m,ques2.m,ques3.m,ques4.m）實驗二實驗?zāi)康?加深對DFT和FFT的了解實驗原理與方法:1. DFT的定義： 2.利用DFT對連續(xù)時間信號進行頻譜分析： 連續(xù)時間信號的傅里葉變換所得信號的頻譜函數(shù)是頻率的連續(xù)函數(shù)，而序列的離散傅里葉變換是一種時域與頻域均離散化了的變換，適合做數(shù)值運算，成為分析離散時間信號和系統(tǒng)的有力工具。而由于DFT具有選頻特性，所以常用它對連續(xù)時間信號進行頻譜分析。 利用DFT對連續(xù)時間信號進行頻譜分析的過程如下圖所示：圖2.1 利用DFT對連續(xù)時間信號進行頻譜分析的過程 從圖中可以看出，利用DFT對連續(xù)時間信號進行頻譜分析是

18、近似的，其近似程度與信號帶寬，采樣頻率和截取長度等參數(shù)有關(guān)。 此外，利用DFT對連續(xù)時間信號進行頻譜分析過程中會出現(xiàn)混疊現(xiàn)象，截斷效應(yīng)，柵欄效應(yīng)以及信號的頻譜分辨力等問題，與之對應(yīng)分別有相應(yīng)的解決方法，對此需要進行學(xué)習(xí)和理解。 3.FFT快速離散傅里葉變換： 序列的FFT算法就是把長序列的DFT逐步分解成幾個短序列的DFT，并利用的周期性和對稱性來減少DFT的運算次數(shù)。最常用的是基2FFT算法，它又分為時域抽取法和頻域抽取法。 實驗內(nèi)容及步驟現(xiàn)有一信號s(t)，有四個頻率分量1050KHz、1100KHz、1200KHz、1500KHz，請完成以下實驗。注意：（1）s(t)請用實信號表示；（2

19、）不允許使用matlab中自帶的DFT及相關(guān)函數(shù)，除非需要對所編寫函數(shù)進行必要的驗證。根據(jù)DFT的定義編寫一個DFT函數(shù)和IDFT函數(shù)，并對s(n)進行DFT運算，對s(k)，即DFTs(n)進行IDFT運算，驗證所編寫函數(shù)的正確性。要求：s(k)與s(t)的頻譜形式一致，即在-范圍內(nèi)對表示。對信號s(t)進行數(shù)字信號分析時，需要對DSP系統(tǒng)參數(shù)盡心謹慎的選擇。（1）如果我們只檢測到了5ms的信號，通過DFT操作是否可以將四個頻率分量全部檢測到，如果可以檢測到，請作圖說明，如果未檢測到，有漏掉的頻率分量，請繼續(xù)完成以下問題；（2）如果你決定進行第2步實驗，說明我們需要在參數(shù)選擇上進行做出改變了

20、，通常我們會有以下措施：a）補零；b）改變采樣率；c）增加被分析信號的有效長度。通過實驗說明上述三種措施，哪些可以將被觀測信號s(t)中的頻率分量全部分析出來，哪些不可以。要提供理論分析和相關(guān)的實驗仿真。1965年，J.W.Cooley與J.W.Tukey提出了DFT的快速算法，即FFT，大幅度的降低了DFT的運算量。（1）請自行編寫一個DIT-基2FFT函數(shù)和DIF-基2FFT函數(shù)，對自定義的一個N=8的序列x(n)進行FFT操作，并對比兩個函數(shù)的輸出結(jié)果，最后通過matlab自帶的FFT函數(shù)對自己所編寫的兩個函數(shù)進行驗證。（2）基于（1）中的兩個FFT函數(shù)，分別通過兩個方法實現(xiàn)IFFT操作

21、，并通過與原x(n)的對比驗證IFFT操作的正確性。實驗結(jié)果分析及結(jié)論總結(jié) 據(jù)DFT的定義編寫一個DFT函數(shù)和IDFT函數(shù)，并對s(n)進行DFT運算，對DFTs(n)進行IDFT運算，驗證所編寫函數(shù)的正確性。A. 主函數(shù)：clear all;close all;clc;dt=0.00001;tf=300;t=0:dt:tf; %tf決定t的取值范圍，可更改st=sin(2*pi*1050*t/20000)+sin(2*pi*1100*t/20000)+sin(2*pi*1200*t/20000)+sin(2*pi*1500*t/20000); %連續(xù)信號s(t)subplot(221)plo

22、t(t,st);title('連續(xù)信號s(t)');xlabel('t');ylabel('s(t)'); %繪制連續(xù)信號圖像Ts=1; %采樣周期，可更改n=0:1:tf/Ts; %采樣點n的向量sn=sin(2*pi*1050*n/20000)+sin(2*pi*1100*n/20000)+sin(2*pi*1200*n/20000)+sin(2*pi*1500*n/20000); %采樣信號s(n)subplot(222)stem(n,sn);title('采樣信號s(n)');xlabel('n');yla

23、bel('s(n)'); %繪制采樣信號圖像Sk=dft(sn); %s(n)做DFT得到S(k)s1n=idft(Sk); %S(k)做IDFT得到s1(n)B. DFT函數(shù)：function Xk= dft(xn)%計算離散傅里葉變換%Xk=dft(xn)%Xk為在n在0,N-1區(qū)間的DFT系數(shù)數(shù)組%xn為N點有限持續(xù)時間序列N=length(xn);disp('輸出N');disp(N); %有效長度N的輸出n=0:1:N-1; %n的行向量k=0:1:N-1; %k的行向量W=exp(-j*2*pi/N); %Wn因子nk=n'*k; %產(chǎn)生一個

24、含nk值的N乘N維矩陣Xk=xn*(W.nk); %DFT系數(shù)的行向量disp('輸出Xk');disp(Xk); %Sk的輸出subplot(223)stem(k,Xk);title('s(n)做DFT所得S(k)');xlabel('k');ylabel('S(k)');%DFT的圖像繪制endC. IDFT函數(shù)：function xn= idft(Xk)%計算逆離散傅里葉變換% xn=idft(Xk,N)%xn為n在0,N-1區(qū)間的N點有限持續(xù)時間序列%Xk為k在0,N-1區(qū)間的DFT系數(shù)數(shù)組N=length(Xk);dis

25、p('序列有效長度');disp(N);n=0:1:N-1; %n的行向量k=0:1:N-1; %k的行向量W=exp(-j*2*pi/N); %Wn因子nk=n'*k; %產(chǎn)生一個含nk值的N乘N維矩陣xn=Xk*(W.-nk)/N; %IDFT系數(shù)的行向量disp('輸出xn');disp(xn); %sn的輸出subplot(224)stem(n,xn);title('S(k)做IDFT所得s1(n)');xlabel('n');ylabel('s1(n)'); %IDFT圖像繪制endD. 實驗結(jié)果

26、：圖2.2 DFT和IDFT圖像E. 由圖像可以看出對s(n)進行DFT得到的S（k）經(jīng)過IDFT運算后所得s1(n)與序列s(n)一樣，所以所編寫的程序正確。（1） 只檢測到了5ms的信號，通過DFT操作是否可以將四個頻率分量全部檢測到i. 代碼 clear all;close all;clc;dt=0.00001;tf=0.005;t=0:dt:tf; %tf決定t的取值范圍，可更改st=sin(2*pi*1050*t/20000)+sin(2*pi*1100*t/20000)+sin(2*pi*1200*t/20000)+sin(2*pi*1500*t/20000); %連續(xù)信號s(t)

27、subplot(221)plot(t,st);title('連續(xù)信號s(t)');xlabel('t');ylabel('s(t)'); %繪制連續(xù)信號圖像Ts=1; %采樣周期，可更改n=0:1:tf/Ts; %采樣點n的向量sn=sin(2*pi*1050*n/20000)+sin(2*pi*1100*n/20000)+sin(2*pi*1200*n/20000)+sin(2*pi*1500*n/20000); %采樣信號s(n)subplot(222)stem(n,sn);title('采樣信號s(n)');xlabel(&

28、#39;n');ylabel('s(n)'); %繪制采樣信號圖像Sk=dft(sn); %s(n)做DFT得到S(k)ii. 圖像圖2.3 只檢測到5ms信號的s(t)的DFTiii. 結(jié)論：不能夠檢測到，需進行第二問，（2） a）補零；b）改變采樣率；c）增加被分析信號的有效長度。通過實驗說明上述三種措施。a) 補零：clear all;close all;clc;dt=0.00001;tf=0.005;t=0:dt:tf; %tf決定t的取值范圍，可更改st=sin(2*pi*1050*t/20000)+sin(2*pi*1100*t/20000)+sin(2*p

29、i*1200*t/20000)+sin(2*pi*1500*t/20000); %連續(xù)信號s(t)subplot(221)plot(t,st);title('連續(xù)信號s(t)');xlabel('t');ylabel('s(t)'); %繪制連續(xù)信號圖像Ts=1; %采樣周期，可更改n=0:1:tf/Ts; %采樣點n的向量sn=sin(2*pi*1050*n/20000)+sin(2*pi*1100*n/20000)+sin(2*pi*1200*n/20000)+sin(2*pi*1500*n/20000); %采樣信號s(n)n1=0:1:t

30、f/Ts+1000;sn1=sn,zeros(1,1000);subplot(222)stem(n1,sn1);title('采樣信號s(n)');xlabel('n');ylabel('s(n)'); %繪制采樣信號圖像Sk=dft(sn1); %s(n)做DFT得到S(k)title('補零所得圖像')圖2.4 補零所得圖像b) 改變采樣頻率：clear all;close all;clc;dt=0.00001;tf=0.005;t=0:dt:tf; %tf決定t的取值范圍，可更改st=sin(2*pi*1050*t/2000

31、0)+sin(2*pi*1100*t/20000)+sin(2*pi*1200*t/20000)+sin(2*pi*1500*t/20000); %連續(xù)信號s(t)subplot(221)plot(t,st);title('連續(xù)信號s(t)');xlabel('t');ylabel('s(t)'); %繪制連續(xù)信號圖像Ts=0.00001; %采樣周期，可更改n=0:1:tf/Ts; %采樣點n的向量sn=sin(2*pi*1050*n/20000)+sin(2*pi*1100*n/20000)+sin(2*pi*1200*n/20000)+si

32、n(2*pi*1500*n/20000); %采樣信號s(n)subplot(222)stem(n,sn);title('采樣信號s(n)');xlabel('n');ylabel('s(n)'); %繪制采樣信號圖像Sk=dft(sn); %s(n)做DFT得到S(k)title('改變采樣頻率所得圖像')圖2.5 改變采樣頻率所得圖像c) 增加被分析信號有效長度：clear all;close all;clc;dt=0.00001;tf=0.005;t=0:dt:tf; %tf決定t的取值范圍，可更改st=sin(2*pi*1

33、050*t/20000)+sin(2*pi*1100*t/20000)+sin(2*pi*1200*t/20000)+sin(2*pi*1500*t/20000); %連續(xù)信號s(t)subplot(221)plot(t,st);title('連續(xù)信號s(t)');xlabel('t');ylabel('s(t)'); %繪制連續(xù)信號圖像Ts=1; %采樣周期，可更改n=0:1:tf/Ts; %采樣點n的向量sn=sin(2*pi*1050*n/20000)+sin(2*pi*1100*n/20000)+sin(2*pi*1200*n/20000

34、)+sin(2*pi*1500*n/20000); %采樣信號s(n)subplot(222)stem(n,sn);title('采樣信號s(n)');xlabel('n');ylabel('s(n)'); %繪制采樣信號圖像N=length(sn)+1000;n=0:1:N-1; %n的行向量k=0:1:N-1; %k的行向量W=exp(-j*2*pi/N); %Wn因子nk=n'*k; %產(chǎn)生一個含nk值的N乘N維矩陣Sk=sn*(W.nk); %DFT系數(shù)的行向量disp('輸出Xk');disp(Sk); %Sk的

35、輸出subplot(223)stem(k,Sk);title('s(n)做DFT所得S(k)');xlabel('k');ylabel('S(k)');%DFT的圖像繪制title('增加被分析信號有效長度所得圖像') 圖2.6 增加被分析信號有效長度所得圖像（1） 編寫一個DIT-基2FFT函數(shù)和DIF-基2FFT函數(shù)，對自定義的一個N=8的序列x(n)進行FFT操作，并對比兩個函數(shù)的輸出結(jié)果，最后通過matlab自帶的FFT函數(shù)對自己所編寫的兩個函數(shù)進行驗證。i. DIT-基2FFT函數(shù)a) 代碼：clc;close all;

36、clear;format compact; xn=1 3 5 9 6 7 0 3; %初始序列xnM=nextpow2(length(xn), N=2M, for m=0:N/2-1; %旋轉(zhuǎn)因子指數(shù)范圍 WN(m+1)=exp(-j*2*pi/N)m; %計算旋轉(zhuǎn)因子 end A=xn,zeros(1,N-length(xn); %數(shù)據(jù)輸入，長度不夠補零disp('輸入到各存儲單元的數(shù)據(jù):'),disp(A); J=N/2; %給倒序數(shù)賦初值 for I=1:N/2-1; %按序交換數(shù)據(jù)和算倒序數(shù) if I<J; %條件判斷及數(shù)據(jù)交換 T=A(I+1);A(I+1)=A

37、(J+1);A(J+1)=T; end K=N/2; while J>=K; J=J-K;K=K/2; end J=J+K; end %disp('倒序后各存儲單元的數(shù)據(jù):'),disp(A); for L=1:M; %分級計算 % disp('運算級次:'),disp(L); B=2(L-1); for R=0:B-1; %各級按序蝶算 P=2(M-L)*R; for K=R:2L:N-2; %每序依次計算 T=A(K+1)+A(K+B+1)*WN(P+1); A(K+B+1)=A(K+1)-A(K+B+1)*WN(P+1); A(K+1)=T; end

38、 end %disp('本級運算后各存儲單元的數(shù)據(jù):'),disp(A); end %disp('輸出各存儲單元的數(shù)據(jù):')Xk1=A;subplot(211)stem(Xk1); title('8點DITFFT計算X(k)')%disp('調(diào)用fft函數(shù)運算的結(jié)果:'),Xk2=fft(xn,N);subplot(212)stem(Xk2);title('對結(jié)果進行FFT驗證X(k)')b) 圖像結(jié)果：圖2.7 DIT-基2FFTc) 結(jié)論：兩者計算結(jié)果一致，故程序正確。ii. DIF-基2FFT函數(shù)a) 代碼：

39、clear allclose allclcxn=1 3 5 9 6 7 0 3; %初始序列xnN=8;n=0:N-1;m=log2(N); %蝶形運算的次數(shù)d=bin2dec(fliplr(dec2bin(0:N-1,m)+1; %將其倒置排序Xk = xn(d); %將到序排列的x作為Xk初始值for M=1:m % 將DFT作m次基2分解,從左到右，對每次分解作DFT運算 Nm=2M;u=1; % 旋轉(zhuǎn)因子u初始化為WN0=1 WN=exp(-i*2*pi/Nm); % 本次分解的基本DFT因子WN=exp(-i*2*pi/Nmr) for p=1:Nm/2 % 本次跨越間隔內(nèi)的各次蝶形

40、運算 for k=p:Nm:N % 本次蝶形運算的跨越間隔為Nmr=2mm kp=k+Nm/2; % 確定蝶形運算的對應(yīng)單元下標 a=Xk(kp)*u; % 蝶形運算的乘積項 Xk(kp)=Xk(k)-a; % 蝶形運算 Xk(k)=Xk(k)+a; % 蝶形運算 end u=u*WN; % 修改旋轉(zhuǎn)因子,多乘一個基本DFT因子WN end end y=Xk;ditfft=abs(y);subplot(211);stem(n,ditfft,'b'); %繪制傅里葉變換xlabel('k');ylabel('X(k)');title('8點

41、DIFFFT計算X(k)')grid on;check=fft(xn,N);fft_abs=abs(check);subplot(212)stem(n,fft_abs,'b'); %進行快速傅里葉驗證xlabel('k');ylabel('X(k)');title('對結(jié)果進行FFT驗證X(k)')grid on;b) 圖像結(jié)果：圖2.8 DIF-基2FFTc) 結(jié)論：兩者計算結(jié)果一致，故程序正確。（2）基于（1）中的兩個FFT函數(shù)，分別通過兩個方法實現(xiàn)IFFT操作，并通過與原x(n)的對比驗證IFFT操作的正確性。i.

42、旋轉(zhuǎn)因子指數(shù)變極性法a) 代碼：clc;close all;clear;format compact; Xk=1 3 5 9 6 7 0 3; %初始序列XkM=nextpow2(length(Xk), N=2M, for m=0:N/2-1; %旋轉(zhuǎn)因子指數(shù)范圍 WN(m+1)=exp(j*2*pi/N)m; %計算旋轉(zhuǎn)因子(改變極性) end A=Xk,zeros(1,N-length(Xk); %數(shù)據(jù)輸入，長度不夠補零disp('輸入到各存儲單元的數(shù)據(jù):'),disp(A); J=N/2; %給倒序數(shù)賦初值 for I=1:N/2-1; %按序交換數(shù)據(jù)和算倒序數(shù) if I

43、<J; %條件判斷及數(shù)據(jù)交換 T=A(I+1);A(I+1)=A(J+1);A(J+1)=T; end K=N/2; while J>=K; J=J-K;K=K/2; end J=J+K; end %disp('倒序后各存儲單元的數(shù)據(jù):'),disp(A); for L=1:M; %分級計算 % disp('運算級次:'),disp(L); B=2(L-1); for R=0:B-1; %各級按序蝶算 P=2(M-L)*R; for K=R:2L:N-2; %每序依次計算 T=A(K+1)+A(K+B+1)*WN(P+1); A(K+B+1)=A(K

44、+1)-A(K+B+1)*WN(P+1); A(K+1)=T; end end %disp('本級運算后各存儲單元的數(shù)據(jù):'),disp(A); end %disp('輸出各存儲單元的數(shù)據(jù):')xn1=A*1/N; %結(jié)果乘以1/Nsubplot(211)stem(xn1); title('旋轉(zhuǎn)因子指數(shù)變極性法實現(xiàn)IFFT計算xn')%disp('調(diào)用fft函數(shù)運算的結(jié)果:'),xn2=ifft(Xk,N);subplot(212)stem(xn2);title('對結(jié)果進行IFFT驗證xn')b) 圖像結(jié)果：圖2

45、.9 旋轉(zhuǎn)因子指數(shù)變極性法c) 結(jié)論：通過圖像可以看出程序編寫無誤。ii. 直接調(diào)用FFT子程序法：a) 代碼：clc;close all;clear;format compact; Xk=1 3 5 9 6 7 0 3; %初始序列XkXk1=conj(Xk); %取初始序列的共軛%調(diào)用FFT子程序計算M=nextpow2(length(Xk1), N=2M, for m=0:N/2-1; %旋轉(zhuǎn)因子指數(shù)范圍 WN(m+1)=exp(-j*2*pi/N)m; %計算旋轉(zhuǎn)因子 end A=Xk1,zeros(1,N-length(Xk1); %數(shù)據(jù)輸入，長度不夠補零disp('輸入到各

46、存儲單元的數(shù)據(jù):'),disp(A); J=N/2; %給倒序數(shù)賦初值 for I=1:N/2-1; %按序交換數(shù)據(jù)和算倒序數(shù) if I<J; %條件判斷及數(shù)據(jù)交換 T=A(I+1);A(I+1)=A(J+1);A(J+1)=T; end K=N/2; while J>=K; J=J-K;K=K/2; end J=J+K; end %disp('倒序后各存儲單元的數(shù)據(jù):'),disp(A); for L=1:M; %分級計算 % disp('運算級次:'),disp(L); B=2(L-1); for R=0:B-1; %各級按序蝶算 P=2

47、(M-L)*R; for K=R:2L:N-2; %每序依次計算 T=A(K+1)+A(K+B+1)*WN(P+1); A(K+B+1)=A(K+1)-A(K+B+1)*WN(P+1); A(K+1)=T; end end %disp('本級運算后各存儲單元的數(shù)據(jù):'),disp(A); end %disp('輸出各存儲單元的數(shù)據(jù):')xn1=conj(A)*1/N; %共軛并乘以1/Nsubplot(211)stem(xn1); title('直接調(diào)用FFT子程序法實現(xiàn)IFFT計算xn')%disp('調(diào)用fft函數(shù)運算的結(jié)果:'

48、;),xn2=ifft(Xk,N);subplot(212)stem(xn2);title('對結(jié)果進行IFFT驗證xn')b) 圖像結(jié)果： 圖2.10 直接調(diào)用FFT子程序法c) 結(jié)論：由MATLAB圖像對比可知，程序編寫無誤。實驗三實驗?zāi)康模?加深對FIR濾波器和IIR濾波器的認識與理解。實驗原理與方法 數(shù)字濾波器基于單位脈沖響應(yīng)的特性分為無限長單位脈沖響應(yīng)數(shù)字濾波器（IIR）和有限長單位脈沖響應(yīng)數(shù)字濾波器（FIR）。 IIR濾波器可以利用成熟的模擬濾波器進行設(shè)計，但是是非線性相位，常用方法有脈沖響應(yīng)不變法和雙極性不變法；FIR濾波器可嚴格線性相位并能夠任意幅度特性，且為因

49、果穩(wěn)定系統(tǒng)，可用FFT計算，但是階次比IIR高的多，常用方法窗函數(shù)法。實驗內(nèi)容及步驟1、利用設(shè)計IIR低通、帶通、高通三個濾波器，對音頻信號中的低頻部分、中頻部分、高頻部分進行提??；2、利用設(shè)計FIR低通、帶通、高通三個濾波器，對音頻信號中的低頻部分、中頻部分、高頻部分進行提?。?、對FIR濾波器和IIR濾波器的處理效果進行對比，包括：感性對比和理性分析。4、試將音頻中所提取的各部分的能量進行改變，再將各部分重組，比較一下重組后的音頻與原始音頻的變化。實驗結(jié)果分析及結(jié)論總結(jié)1. 利用設(shè)計IIR低通、帶通、高通三個濾波器，對音頻信號中的低頻部分、中頻部分、高頻部分進行提??；i. 代碼：a) 主函

50、數(shù)：y,Fs=audioread('E曉霞m文件實驗3大王叫我來巡山.wav');%讀取音頻，獲得采樣頻率N=length(y); %采樣點數(shù)T=N/Fs; %采樣時間turn=(0:N-1)/Fs; %定義時間軸furn=(-N/2:N/2-1)/N*Fs; %定義頻率軸% 未處理前信號波形figure(1)subplot(2,1,1)plot(turn,y(:,1) %畫出時域波形圖xlabel('t'),ylabel('s')title('音頻信號時域波形圖')yf=fftshift(fft(y(:

51、,1); %對信號做快速傅里葉變換subplot(2,1,2)plot(furn,abs(yf) %畫出頻域波形圖xlabel('f'),ylabel('s')title('音頻信號頻域波形圖')ylow=IIR_LP(y,Fs,furn); %低頻部分ybp=IIR_BP(y,Fs,furn); %中頻部分yhigh=IIR_HP(y,Fs,furn); %高頻部分b) 低通：function ylow=IIR_LP(y,Fs,furn)%橢圓型低通濾波器設(shè)計Wp=2*1000/Fs; %通帶截止頻率Ws=2*1110/Fs; %阻帶截止頻率Rp=3; %通帶內(nèi)最大衰減系數(shù)Rs=60; %阻帶內(nèi)最大衰減系數(shù)No,Wc=ellipord(Wp,Ws,Rp,Rs); %設(shè)計濾波器的階數(shù)以及3dB截止頻率Bz,Az=ellip(No,Rp,Rs,Wc,'low'); %計算濾波器系統(tǒng)函數(shù)分子分母多項式 ylow=filter(Bz,Az,y); % 進行低通濾波ylowf1=fftshift(fft(ylow(:,1); % 對信號做FFT變換 figure(2)plot(furn,abs(ylowf1)xlabel('f'),ylabel('s')title('信號低