數字信號處理論文

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

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

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

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

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

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

7、畫圖說明。實驗結果分析及結論總結 1序列的產生：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軸代表參數title('單位階躍序列') %添加曲線名稱 圖1.1 單位階躍序列 b) 復指數序列：n=-10:1:10;A=1;alpha=-0.1+0.3j; x=A

8、*exp(alpha*n); %建立復指數序列向量 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('復指數序列');xlabel('n') %依次繪制復

9、指數序列的實部，虛部，振幅以及自身曲線圖1.2 復指數序列c) 實指數序列：clear all;close all;clc;n=0:30;a1=0.8;a2=-0.8;a3=1.5;a4=-1.5; %取不同底數x1=a1.n;x2=a2.n;x3=a3.n;x4=a4.n; %建立實指數序列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 實指數序列 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.復雜信號的產生：clear allclose allclcR=30;n=1:R-1;A=2;m=0.5;x_l=cos(2*pi*0.8*n); %建立調制信號x_h=cos(2*pi*1323*n); %建立高頻載波信號y(n)=A*(1+m.*x_l).*x_h; %建立調幅信號stem(n,y(n) %繪制調幅曲線title('調幅廣播')圖1.7 調幅廣播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 線性卷積驗證結果：計算結果與對位相乘法計算所得結果一致。 思考題1) 可以通過什么參數控制序列的增長率與衰減率？哪個參數可以控制序列的幅值？圖1.9 復指數序列振幅控制因素圖1.10 復指數序列變化率控制因素通過MATLAB繪制圖形結果得出：對于復指數序列，可以控制序列的幅值，控制序列的增長率與衰減率。2) 復指數序列的實部和虛部分別是什么？圖1.11 復指數序列的實部與虛部特點 通過MATLAB繪制圖形結果得出： 當時，復指數序列的實部

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

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

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

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

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

21、，并通過與原x(n)的對比驗證IFFT操作的正確性。實驗結果分析及結論總結 據DFT的定義編寫一個DFT函數和IDFT函數，并對s(n)進行DFT運算，對DFTs(n)進行IDFT運算，驗證所編寫函數的正確性。A. 主函數：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函數：function Xk= dft(xn)%計算離散傅里葉變換%Xk=dft(xn)%Xk為在n在0,N-1區(qū)間的DFT系數數組%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; %產生一個

24、含nk值的N乘N維矩陣Xk=xn*(W.nk); %DFT系數的行向量disp('輸出Xk');disp(Xk); %Sk的輸出subplot(223)stem(k,Xk);title('s(n)做DFT所得S(k)');xlabel('k');ylabel('S(k)');%DFT的圖像繪制endC. IDFT函數：function xn= idft(Xk)%計算逆離散傅里葉變換% xn=idft(Xk,N)%xn為n在0,N-1區(qū)間的N點有限持續(xù)時間序列%Xk為k在0,N-1區(qū)間的DFT系數數組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; %產生一個含nk值的N乘N維矩陣xn=Xk*(W.-nk)/N; %IDFT系數的行向量disp('輸出xn');disp(xn); %sn的輸出subplot(224)stem(n,xn);title('S(k)做IDFT所得s1(n)');xlabel('n');ylabel('s1(n)'); %IDFT圖像繪制endD. 實驗結果

26、：圖2.2 DFT和IDFT圖像E. 由圖像可以看出對s(n)進行DFT得到的S（k）經過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. 結論：不能夠檢測到，需進行第二問，（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; %產生一個含nk值的N乘N維矩陣Sk=sn*(W.nk); %DFT系數的行向量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函數和DIF-基2FFT函數，對自定義的一個N=8的序列x(n)進行FFT操作，并對比兩個函數的輸出結果，最后通過matlab自帶的FFT函數對自己所編寫的兩個函數進行驗證。i. DIT-基2FFT函數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; %旋轉因子指數范圍 WN(m+1)=exp(-j*2*pi/N)m; %計算旋轉因子 end A=xn,zeros(1,N-length(xn); %數據輸入，長度不夠補零disp('輸入到各存儲單元的數據:'),disp(A); J=N/2; %給倒序數賦初值 for I=1:N/2-1; %按序交換數據和算倒序數 if I<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('倒序后各存儲單元的數據:'),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('本級運算后各存儲單元的數據:'),disp(A); end %disp('輸出各存儲單元的數據:')Xk1=A;subplot(211)stem(Xk1); title('8點DITFFT計算X(k)')%disp('調用fft函數運算的結果:'),Xk2=fft(xn,N);subplot(212)stem(Xk2);title('對結果進行FFT驗證X(k)')b) 圖像結果：圖2.7 DIT-基2FFTc) 結論：兩者計算結果一致，故程序正確。ii. DIF-基2FFT函數a) 代碼：

39、clear allclose allclcxn=1 3 5 9 6 7 0 3; %初始序列xnN=8;n=0:N-1;m=log2(N); %蝶形運算的次數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; % 旋轉因子u初始化為WN0=1 WN=exp(-i*2*pi/Nm); % 本次分解的基本DFT因子WN=exp(-i*2*pi/Nmr) for p=1:Nm/2 % 本次跨越間隔內的各次蝶形

40、運算 for k=p:Nm:N % 本次蝶形運算的跨越間隔為Nmr=2mm kp=k+Nm/2; % 確定蝶形運算的對應單元下標 a=Xk(kp)*u; % 蝶形運算的乘積項 Xk(kp)=Xk(k)-a; % 蝶形運算 Xk(k)=Xk(k)+a; % 蝶形運算 end u=u*WN; % 修改旋轉因子,多乘一個基本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('對結果進行FFT驗證X(k)')grid on;b) 圖像結果：圖2.8 DIF-基2FFTc) 結論：兩者計算結果一致，故程序正確。（2）基于（1）中的兩個FFT函數，分別通過兩個方法實現IFFT操作，并通過與原x(n)的對比驗證IFFT操作的正確性。i.

42、旋轉因子指數變極性法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; %旋轉因子指數范圍 WN(m+1)=exp(j*2*pi/N)m; %計算旋轉因子(改變極性) end A=Xk,zeros(1,N-length(Xk); %數據輸入，長度不夠補零disp('輸入到各存儲單元的數據:'),disp(A); J=N/2; %給倒序數賦初值 for I=1:N/2-1; %按序交換數據和算倒序數 if I

43、<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('倒序后各存儲單元的數據:'),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('本級運算后各存儲單元的數據:'),disp(A); end %disp('輸出各存儲單元的數據:')xn1=A*1/N; %結果乘以1/Nsubplot(211)stem(xn1); title('旋轉因子指數變極性法實現IFFT計算xn')%disp('調用fft函數運算的結果:'),xn2=ifft(Xk,N);subplot(212)stem(xn2);title('對結果進行IFFT驗證xn')b) 圖像結果：圖2

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

46、存儲單元的數據:'),disp(A); J=N/2; %給倒序數賦初值 for I=1:N/2-1; %按序交換數據和算倒序數 if I<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('倒序后各存儲單元的數據:'),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('本級運算后各存儲單元的數據:'),disp(A); end %disp('輸出各存儲單元的數據:')xn1=conj(A)*1/N; %共軛并乘以1/Nsubplot(211)stem(xn1); title('直接調用FFT子程序法實現IFFT計算xn')%disp('調用fft函數運算的結果:'

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

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

50、數：y,Fs=audioread('E曉霞m文件實驗3大王叫我來巡山.wav');%讀取音頻，獲得采樣頻率N=length(y); %采樣點數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)%橢圓型低通濾波器設計Wp=2*1000/Fs; %通帶截止頻率Ws=2*1110/Fs; %阻帶截止頻率Rp=3; %通帶內最大衰減系數Rs=60; %阻帶內最大衰減系數No,Wc=ellipord(Wp,Ws,Rp,Rs); %設計濾波器的階數以及3dB截止頻率Bz,Az=ellip(No,Rp,Rs,Wc,'low'); %計算濾波器系統函數分子分母多項式 ylow=filter(Bz,Az,y); % 進行低通濾波ylowf1=fftshift(fft(ylow(:,1); % 對信號做FFT變換 figure(2)plot(furn,abs(ylowf1)xlabel('f'),ylabel('s')title('信號低