第一篇：數字信號處理課設

信息科學與工程學院

數字信號處理課程設計實驗報告

課題名稱： 簡單信號濾波演示系統 學生姓名： 學 號： 專業班級： 指導老師： 實驗時間： 2014.10.8

目 錄

第一章

概述.................................3

第二章

第三章

第四章第五章第六章

1.1 FIR、IIR概述.................................3 1.2題目要求......................................3 設計分析.............................5 2.1算法分析......................................5 2.2 在matlab中實現的分析........................6

程序實現.......................8 3.1 程序主體介紹..................................8 3.2 子程序........................................9 3.3 程序調試及運行結果............................9 3.4 結果分析及問題分析...........................16 心得體會............................17 參考文獻............................18 源代碼..............................19

MATLAB

第一章 概述

1.1 FIR、IIR概述

數字濾波器是指輸入輸出均為數字信號，通過一定的運算關系改變輸入信號所含頻率成分的相對比例或者濾除某些頻率成分的器件。數字濾波器與模擬濾波器相比數字濾波器具有精度高、穩定、體積小、重量小、靈活等特點。主要分為兩種：有限脈沖響應FIR和無限脈沖響應IIR。設計濾波器的主要要求有兩種，一是幅頻特性，一是相頻特性。一般的濾波器主要是對幅頻特性作出要求，如果對輸出相頻特性也有要求，就需要用到線性相位濾波器。IIR濾波器的設計主要有兩類，一是借助于模擬濾波器設計進行，二是直接在頻域或時域中進行設計。FIR濾波器的設計不能借助于模擬濾波器，也有兩類設計方法，一是窗函數法，二是頻率采樣法。還有一種比較有效的方法是切比雪夫等波紋逼近法，需通過計算機輔助進行。

1.2題目要求

設計一個工作流程如圖所示的信號濾波演示系統：

圖2 濾波演示圖

⑴ 信號發生器—根據信號選擇分為兩大類：

① 靜態型：直接輸入(或從文件讀取)測試信號序列；

② 動態型：輸入由多個不同頻率正弦信號疊加組合而成的模擬信號公式 100sin(2?f1t)?100sin(2?f2t)?...?100sin(2?fnt)采樣頻率（Hz）以及采樣點數，動態生成該信號的采樣序列，作為測試信號。

⑵ 頻譜分析—使用FFT 對產生的測試信號進行頻譜分析并展示其幅頻、相頻特

性，指定需要濾除或保留的頻帶，通過選擇濾波器類型（IIR/FIR），確定對

應的濾波器（低通、高通、帶通、帶阻）技術指標。

⑶ 濾波器設計—根據IIR/FIR 數字濾波器技術指標設計濾波器，生成相應的濾

波器系數，并展示對應的濾波器幅頻（衰減）、相頻特性。

① IIR DF 設計：使用雙線性變換法，可選擇濾波器基型（巴特沃斯或切比雪夫型）；

② FIR DF 設計：使用窗口法，可選擇窗口類型。

⑷ 數字濾波—根據設計的濾波器系數，對測試信號進行濾波。

① IIR DF：要求通過差分方程迭代實現濾波，未知初值置零處理； ② FIR DF：要求通過快速卷積實現濾波。可以選擇使用重疊相加或重疊保留法進行卷積運算，并動態展示卷積運算的詳細過程。

⑸ 輸出信號分析—展示濾波后信號的幅頻與相頻特性，分析是否滿足濾波要求。

對同一濾波要求，根據輸出信號頻譜，對比分析各類濾波器的差異。選做部分

將一段語音作為測試信號，通過頻譜展示和語音播放，對比分析濾波前后語音信號的變化，進一步加深對數字信號處理的理解。

第二章 設計分析

2.1算法分析

此題目的實現可分為三個某塊的設計實現：輸入信號模塊，設計濾波器模塊，濾波模塊。

首先明確各模塊間的數據依賴關系：在輸入信號模塊得到信號后，對信號

進行頻域分析，從而確定濾波器的相關技術指標，根據濾波器的技術指標與類型，在濾波器設計模塊完成濾波器的設計，然后將濾波器的設計結果傳遞給濾波模塊，濾波模塊同時接收輸入信號，從而通過運算，實現信號的濾波處理。從數據傳遞關系上分析，濾波模塊的實現依賴于其他模塊的數據輸出，因此放在最后設計。先設計輸入模塊。

因為此設計相對復雜，分模塊設計，通過參數傳遞和接口實現分模塊設計即檢驗，提高程序的穩定性與健壯性。

輸入的實現可以有兩種方式：靜態輸入和動態輸入。靜態輸入選擇從文本輸入數據，將信號取樣值以矩陣的形式存放在文本中。采用文件的讀取就可以實現，比較容易。動態輸入，即輸入由一系列頻率的正弦信號相加的組成的信號，需要經過采樣的，注意在設置采樣頻率時一定要符合奈奎斯特準則，提高采樣點數，增加頻譜分辨率。最后輸出一采樣信號向量，傳遞給其余兩模塊。

濾波器的設計，通過輸入信號的頻譜分析，設置濾波器的參數，然后才可以設計濾波器。第一步需要總結設計濾波器需要哪些參數，通過學習可以總結，所有濾波器的參數有四個：通帶截止頻率、阻帶截止頻率、通帶最大衰減、阻帶最小衰減。

對濾波器的設計分兩類：FIR和IIR，二則所需的函數及設計方法不同。IIR采用借助于模擬濾波器的方式，包括巴特沃斯濾波器和切比雪夫濾波器兩種類型。FIR采用窗函數方式，有矩形窗、三角窗、漢明窗、漢寧窗、布萊克曼和凱塞窗。通過調用不同的函數來實現濾波器的設計。特別在實現窗函數濾波器時，各個函數的主要區別是不同的頻率采樣，可以通過選擇結構實現，簡化程序。通過濾波器的設計最后可以得到濾波器的系統函數H(z)的系數。分析濾波器的幅頻特性和相頻特性，如果不符合要求重新設定濾波器參數或者換成其他濾波器類型。如果性能符合要求，則將系數傳遞給濾波模塊。濾波模塊不調用濾波函數，實現濾波功能根據濾波器類型的不同，有兩種方式可以選擇，一種是通過差分方程運算，一種是通過線性卷積運算。前者適合對IIR濾波器進行濾波，后者適合對FIR濾波器進行濾波。且線性卷積為實現對長序列的卷積運算，采用重疊相加法，采用動態變化展示運算過程。

2.2 在matlab中實現的分析

輸入模塊通過讀取文件和直接輸入數據運算可以很容易實現。在輸入信號確定后，要對其進行頻譜分析，從而確定濾波器的參數和類型（低通、高通、帶通、帶阻），此模塊作用也就完成，將數據分別用全局變量傳遞給下面的模塊。

設置模擬信號：

我采用的測試信號是兩個正弦信號疊加而成的信號，信號為

y=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts)其中頻率f1=30;頻率f2=50;頻率f3=200;采樣頻率fs=3000;采樣間隔 Ts=1/fs;采樣點數

N=1024;n=1:N 采集模擬信號的程序代碼： f1=30;% 頻率1

f2=50;% 頻率2 f3=200;% 頻率3

fs=3000;% 采樣頻率

Ts=1/fs;% 采樣間隔

N=1024;% 采樣點數

n=1:N;

y=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts)+sin(2*pi*f1*n*Ts）;% 正弦波混合

頻譜分析：

使用FFT對產生的測試信號進行頻譜分析并展示其幅頻特性與相頻特性，指定需要濾除的頻帶，通過選擇濾波器類型IIR，確定對應的濾波器（低通、高通）技術指標Fp、Fc、Rp、Rs。

濾波器的設計：

根據以上技術指標（通帶截止頻率、通帶最大衰減、阻帶截止頻率、阻帶最

小衰減），設計數字濾波器，生成相應的濾波器系數，并畫出對應的濾波器

幅頻特性與相頻特性。分別設計巴特沃斯濾波器、切比雪夫I型濾波器、切

比雪夫II型濾波器、橢圓濾波器、yulewalk濾波器。

巴特沃斯和切比雪夫的濾波函數調用為：

[N,Wc]=buttord(wp,ws,rp,rs);[N,Wc]=cheb1ord(wp,ws,rp,rs);[B,A]=butter(N,Wc,’property’);[B,A]=cheby1(N,rp,Wc,’property’);

property對于低通和高通為’’，帶通’high’,帶阻’stop’;窗函數濾波器設計的調用函數：

求窗函數的階數：

N=ceil((h*pi)/wdel);%wdel為窗函數的過渡帶寬，h對應不同窗函數的值 wn=boxcar(N+1);wn=triang(N+1);wn=hanning(N+1);wn=hamming(N+1);wn=blackman(N+1);wn=kaiser(N+1,beta);%bata為kaiser的a參數 B=fir1(N,ws,'property',wn);property對于低通和高通為’’，帶通’high’,帶阻’stop’;

數字濾波：

根據設計的濾波器系數，對測試信號通過差分方程迭代實現濾波數字濾波，展示濾波后信號的幅頻特性與相頻特性，分析是否滿足濾波要求（對同一濾波要求，對比分析各類濾波器的差異）。

需要注意的是，窗函數對濾波參數的使用，需要通過運算得到一窗函數階數，對通帶最大衰減無使用。需要注意FIR和IIR對參數的不同利用。

濾波模塊：運用差分方程運算的濾波函數，可以用循環調用簡單實現。動態展示線性卷積的重疊相加法可以用流程圖來說明： 第三章 MATLAB程序實現

3.1 程序主體介紹

此程序界面也是采用GUI，不過是采用圖形用戶界面開發環境。和第一個對比，布局更加容易，很容易調整界面。但是在一些函數操作行會受一些限制，比如函數參數的傳遞，目前來說還是采用全局變量才能完成。

pushbutton1_Callback(hObject, eventdata, handles)%獲得輸入信號

uipanel11_SelectionChangeFcn(hObject, eventdata, handles)%獲取濾波器的設計參數

uipanel15_SelectionChangeFcn(hObject, eventdata, handles)%設計濾波器

pushbutton4_Callback(hObject, eventdata, handles)%進行濾波

以上是一些功能控件，需要注意的是他們的view ballcak屬性，對于groupbutton一組按鈕，需要選用SelectionChangeFcn，其余選擇callback屬性即可。另外，handles函數作用為傳遞的句柄，是一結構體，調用時注意此處。其余和直接用函數創建GUI界面無太大區別。

3.2 子程序

3.2.1

IIR濾波器

FIR濾波器

巴特沃斯濾波器

3.2.2

矩形窗

三角窗

漢寧窗

哈明窗 3.3 程序調試及運行結果

原始信號幅頻特性及相頻特性

巴特沃斯的幅頻特性及相頻特性

信號濾波后的幅頻特性及相頻特性

矩形窗的幅頻特性及相頻特性

信號濾波后的幅頻特性及相頻特性

三角窗的幅頻特性及相頻特性

信號濾波后的幅頻特性及相頻特性

漢寧窗的幅頻特性及相頻特性

信號濾波后的幅頻特性及相頻特性

哈明窗的幅頻特性及相頻特性

信號濾波后的幅頻特性及相頻特性

3.4 結果分析及問題分析

用雙線性變換法設計無限脈沖響應數字濾波器（IIF DF）時，先把數字濾波器指標轉換成模擬濾波器的指標，然后根據模擬濾波器的指標設計模擬濾波器，再經過線性變換把模擬濾波器轉換成數字濾波器。該系統要能夠設計巴特沃茲型低通、帶通、高通濾波器，并能夠輸入數字濾波器的性能指標，顯示出濾波器的階數和系數。該系統的關鍵部分是濾波器的設計部分，按照雙線性變換法設計濾波器的步驟進行設計即可。

第四章 心得體會

本次課程設計可以稱為matlab和數字信號處理的綜合設計，因為這次有一半的時間在研究matlab，有一半的時間在思考數字信號問題的解決。幾天下來收獲是很多的，對數字信號處理有了一次全面的回顧，而且也看到了matlab的神秘面容。只要是自己真正的努力去做了，就絕對不會后悔在課程設計上花費的時間。現在發現自己確實會的太少，而這次又學了太多，對自己以后的學習還是有鞭策意義的。

這次課程設計也使我認識到了matlab的強大，或者進一步說工具軟件有著你意想不到的功效，能節省你的時間，同時又能夠讓你從實踐上更深層次的認識到所學知識的奇妙，而且可以同時明白了學習與實踐的相輔相成。對課程設計也是一直保持很高興趣的，雖然有時為它焦頭爛額，但是也會因此有柳暗花明的喜悅，任何事情都折射著一個普通的道理——付出才有回報。淺顯的道理，卻是需要我們用毅力來堅持見證的。

具體來說，通過本次課程設計，我掌握了MATLAB的基本運用，尤其是其中GUI可視化圖形用戶界面的設計，包括函數調用與在圖形用戶界面開發環境下的調用。函數的調用與參數的傳遞是兩個簡單卻很容易出錯的地方，自由細心才可以保證程序的健壯性。

對數字信號本身內容的理解來說，全書的內容確實是可以融合在這兩個課程設計題目中的，一個是DFT運用，另一個是濾波器的設計和利用。對全書的內容可以做出最好的概括。其中濾波器的設計中，調用了許多濾波函數，這是本次實驗有點欠缺的地方，但是仍然從整體上把握了整個濾波器的設計過程。

對于課程設計中出現的問題，解決的方式有兩種，一是自己解決?？梢陨暇W，查閱圖書和matlab的本身的幫助文件，不斷嘗試，自己修改錯誤，可以更好的反省。二是與同學相商。一加一不是等于二那么簡單的，相互交流才是進步的最好方式。其中解決問題最重要的應該是堅持不懈，在遇到問題時，不放棄才有可能稱為最后的勝利者。

每次課程設計都收獲頗多，而且每次都對自己的學習態度做一次調整，這個的作用確實是很大的，希望以后更加珍惜的課程設計的時間。

第五章 參考文獻

[1] 丁玉美等.數字信號處理 [M].西安:西安電子科技大學出版社,2002

[2] 程佩青.數字信號處理教程,第二版[M].北京:清華大學出版社,2001

[3] 趙樹杰等.數字信號處理[M].西安:西安電子科技大學出版社,1997

[4] 陳懷琛等.MATLAB 及在電子信息課程中的應用[M],北京:電子工業出版社出

版,2002

[5] 余成波.數字信號處理及MATLAB實現，第一版[M].北京：清華大學出版社，2008

[6] S.K.Mitra.Digital Signal Processing: A Computer Based Approach, 3rdEdition[M],New York, USA:McGraw-Hill,2000

[7] R.G.Lyons.Understanding Digital Signal Processing,2nd Edition[M].New Jersey, USA:Prentice Hall,2005

第六章 源代碼

function varargout = df(varargin)

gui_Singleton = 1;

gui_State = struct('gui_Name', mfilename,...'gui_Singleton', gui_Singleton,...'gui_OpeningFcn', @df_OpeningFcn,...'gui_OutputFcn', @df_OutputFcn,...'gui_LayoutFcn', [] ,...'gui_Callback', []);if nargin && ischar(varargin{1})

gui_State.gui_Callback = str2func(varargin{1});end

if nargout

[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});else

gui_mainfcn(gui_State, varargin{:});end

function df_OpeningFcn(hObject, eventdata, handles, varargin)

handles.output = hObject;

guidata(hObject, handles);

function varargout = df_OutputFcn(hObject, eventdata, handles)

varargout{1} = handles.output;

function edit1_Callback(hObject, eventdata, handles)

function edit1_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function edit2_Callback(hObject, eventdata, handles)

function edit2_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function edit3_Callback(hObject, eventdata, handles)

function edit3_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function edit4_Callback(hObject, eventdata, handles)

function edit4_CreateFcn(hObject, eventdata, handles)if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function edit5_Callback(hObject, eventdata, handles)

function edit5_CreateFcn(hObject, eventdata, handles)if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function pushbutton1_Callback(hObject, eventdata, handles)f=zeros(1,4);

f(1)=str2num(get(handles.edit2,'String'));f(2)=str2num(get(handles.edit1,'String'));f(3)=str2num(get(handles.edit3,'String'));f(4)=str2num(get(handles.edit4,'String'));

T=str2num(get(handles.edit12,'String'));%采樣周期% fre=str2num(get(handles.edit5,'String'));%采樣頻率% t=0:1/fre:T;

xn=10*sin(2*pi*t*f(1))+10*sin(2*pi*t*f(2))+10*sin(2*pi*t*f(3))+10*sin(2*pi*t*f(4));

set(handles.axes15,'userdata',xn);%將Xn放在用戶數據userdata% yn3=abs(fft(xn));%快速傅立葉變換（符頻特性）% n1=[0:length(yn3)-1]/length(yn3)*fre;%橫坐標% axes(handles.axes12);%坐標系編號% stem(n1,yn3,'.');

axis([0,fre/2,0,max(yn3)]);%坐標軸單位控制% title('信號的幅頻特性');xlabel('頻率（Hz）');ylabel('|X(ejw)|');

yn4=angle(fft(xn));%相頻特性%

n4=[0:length(yn4)-1]/length(yn4)*fre;axes(handles.axes16);stem(n4,yn4,'.');

axis([0,fre/2,min(yn4),max(yn4)]);title('信號的相頻特性');xlabel('頻率（Hz）');ylabel('相位');

function edit6_Callback(hObject, eventdata, handles)

function edit6_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function edit7_Callback(hObject, eventdata, handles)

function edit7_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function edit8_Callback(hObject, eventdata, handles)

function edit8_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function edit9_Callback(hObject, eventdata, handles)

function edit9_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function pushbutton2_Callback(hObject, eventdata, handles)xn=get(handles.axes15,'userdata');fre=str2num(get(handles.edit5,'String'));

fp=str2num(get(handles.edit6,'String'));%通帶最大頻率 % fs=str2num(get(handles.edit7,'String'));%阻帶最小頻率% rp=str2num(get(handles.edit8,'String'));rs=str2num(get(handles.edit9,'String'));wp=2*fp/fre;ws=2*fs/fre;

[N,wc]=buttord(wp,ws,rp,rs,'s');%求階數，截止頻率%

[B,A]=butter(N,wc,'s');%巴特沃茲模擬低通濾波器系數%

[Bz,Az]=bilinear(B,A,1);

[H,w]=freqz(Bz,Az);%分析數字濾波器% axes(handles.axes14);plot(w/pi,abs(H));

axis([0,1,0,max(abs(H))]);title('巴特沃茲的幅頻特性');xlabel('w/π');

ylabel('|X(ejw)|');axes(handles.axes17);plot(w/pi,angle(H));title('巴特沃茲的相頻特性');xlabel('頻率（Hz）');ylabel('相位');

yn2=filter(Bz,Az,xn);%迭代法求解濾波信號% yn=abs(fft(yn2));

n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes15);stem(n1,yn,'.');

title('信號濾波后的幅頻特性');xlabel('頻率（Hz）');ylabel('|X(ejw)|');axis([0,fre/2,0,max(yn)]);yn=angle(fft(yn2));

n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes18);stem(n1,yn,'.');

title('信號濾波后的相頻特性');xlabel('頻率（Hz）');

ylabel('相位');

axis([0,fre/2,min(yn),max(yn)]);

pushbutton1_Callback(hObject, eventdata, handles);

function edit10_Callback(hObject, eventdata, handles)function edit10_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function edit11_Callback(hObject, eventdata, handles)

function edit11_CreateFcn(hObject, eventdata, handles)if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function pushbutton5_Callback(hObject, eventdata, handles)

xn=get(handles.axes15,'userdata');fp=str2num(get(handles.edit10,'String'));fs=str2num(get(handles.edit11,'String'));fre=str2num(get(handles.edit5,'String'));wp=2*fp/fre*pi;ws=2*fs/fre*pi;

Bt=ws-wp;

N0=ceil(1.8*pi/Bt);N=N0+mod(N0+1,2);wc=(wp+ws)/2/pi;

hn=fir1(N-1,wc,boxcar(N));%單位序列相應%

yn=abs(fft(hn));yn=20*log10(yn);

n1=[0:length(yn)-1]/length(yn)*2;axes(handles.axes14);plot(n1,yn);

title('矩形窗的損耗函數');xlabel('w/π');

ylabel('20log|Hg(w)|');axis([0,1,min(yn),max(yn)]);axes(handles.axes17);yn=angle(hn);

n1=[0:length(yn)-1]/length(yn)*2;plot(n1,yn);

title('矩形窗的相頻特性');xlabel('w/π');

ylabel('相位');

axis([0,1,min(yn),max(yn)]);

yn2=fftfilt(hn,xn,512);%重疊相加法求濾波后序列% yn=abs(fft(yn2));n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes15);stem(n1,yn,'.');

title('信號濾波后的幅頻特性');xlabel('頻率（Hz）');ylabel('|X(ejw)|');axis([0,fre/2,0,max(yn)]);yn=angle(fft(yn2));

n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes18);stem(n1,yn,'.');

title('信號濾波后的相頻特性');xlabel('頻率（Hz）');

ylabel('相位');

axis([0,fre/2,min(yn),max(yn)]);

pushbutton1_Callback(hObject, eventdata, handles);

function pushbutton6_Callback(hObject, eventdata, handles)

xn=get(handles.axes15,'userdata');fp=str2num(get(handles.edit10,'String'));fs=str2num(get(handles.edit11,'String'));fre=str2num(get(handles.edit5,'String'));wp=2*fp/fre*pi;ws=2*fs/fre*pi;

Bt=ws-wp;

N0=ceil(6.1*pi/Bt);N=N0+mod(N0+1,2);wc=(wp+ws)/2/pi;

hn=fir1(N-1,wc,bartlett(N));

yn=abs(fft(hn));yn=20*log10(yn);

n1=[0:length(yn)-1]/length(yn)*2;axes(handles.axes14);plot(n1,yn);

title('三角窗的損耗函數');xlabel('w/π');

ylabel('20log|Hg(w)|');axis([0,1,min(yn),max(yn)]);axes(handles.axes17);yn=angle(hn);n1=[0:length(yn)-1]/length(yn)*2;plot(n1,yn);

title('三角窗的相頻特性');xlabel('w/π');

ylabel('相位');

axis([0,1,min(yn),max(yn)]);

yn2=fftfilt(hn,xn,512);yn=abs(fft(yn2));

n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes15);stem(n1,yn,'.');

title('信號濾波后的幅頻特性');xlabel('頻率（Hz）');ylabel('|X(ejw)|');axis([0,fre/2,0,max(yn)]);yn=angle(fft(yn2));

n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes18);stem(n1,yn,'.');

title('信號濾波后的相頻特性');xlabel('頻率（Hz）');

ylabel('相位');

axis([0,fre/2,min(yn),max(yn)]);

pushbutton1_Callback(hObject, eventdata, handles);

function pushbutton7_Callback(hObject, eventdata, handles)

xn=get(handles.axes15,'userdata');fp=str2num(get(handles.edit10,'String'));fs=str2num(get(handles.edit11,'String'));fre=str2num(get(handles.edit5,'String'));wp=2*fp/fre*pi;ws=2*fs/fre*pi;if(fp

N0=ceil(6.2*pi/Bt);N=N0+mod(N0+1,2);wc=(wp+ws)/2/pi;

hn=fir1(N-1,wc,hanning(N));else

Bt=wp-ws;

N0=ceil(6.2*pi/Bt);N=N0+mod(N0+1,2);wc=(wp+ws)/2/pi;

hn=fir1(N-1,wc,'high',hanning(N));end

yn=abs(fft(hn));yn=20*log10(yn);

n1=[0:length(yn)-1]/length(yn)*2;axes(handles.axes14);plot(n1,yn);

title('漢寧窗的損耗函數');xlabel('w/π');

ylabel('20log|Hg(w)|');axis([0,1,min(yn),max(yn)]);axes(handles.axes17);yn=angle(hn);

n1=[0:length(yn)-1]/length(yn)*2;plot(n1,yn);

title('漢寧窗的相頻特性');xlabel('w/π');

ylabel('相位');

axis([0,1,min(yn),max(yn)]);

yn2=fftfilt(hn,xn,512);yn=abs(fft(yn2));

n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes15);stem(n1,yn,'.');

title('信號濾波后的幅頻特性');xlabel('頻率（Hz）');ylabel('|X(ejw)|');axis([0,fre/2,0,max(yn)]);yn=angle(fft(yn2));

n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes18);stem(n1,yn,'.');

title('信號濾波后的相頻特性');xlabel('頻率（Hz）');

ylabel('相位');

axis([0,fre/2,min(yn),max(yn)]);

pushbutton1_Callback(hObject, eventdata, handles);

function pushbutton8_Callback(hObject, eventdata, handles)

xn=get(handles.axes15,'userdata');fp=str2num(get(handles.edit10,'String'));fs=str2num(get(handles.edit11,'String'));fre=str2num(get(handles.edit5,'String'));wp=2*fp/fre*pi;ws=2*fs/fre*pi;if(fp

N0=ceil(6.2*pi/Bt);N=N0+mod(N0+1,2);wc=(wp+ws)/2/pi;

hn=fir1(N-1,wc,hamming(N));else

Bt=wp-ws;

N0=ceil(6.2*pi/Bt);N=N0+mod(N0+1,2);wc=(wp+ws)/2/pi;

hn=fir1(N-1,wc,'high',hamming(N));end

yn=abs(fft(hn));yn=20*log10(yn);

n1=[0:length(yn)-1]/length(yn)*2;axes(handles.axes14);plot(n1,yn);

title('哈明窗的損耗函數');xlabel('w/π');

ylabel('20log|Hg(w)|');axis([0,1,min(yn),max(yn)]);axes(handles.axes17);yn=angle(hn);

n1=[0:length(yn)-1]/length(yn)*2;plot(n1,yn);

title('哈明窗的相頻特性');xlabel('w/π');

ylabel('相位');

axis([0,1,min(yn),max(yn)]);

yn2=fftfilt(hn,xn,512);yn=abs(fft(yn2));

n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes15);stem(n1,yn,'.');

title('信號濾波后的幅頻特性');xlabel('頻率（Hz）');ylabel('|X(ejw)|');axis([0,fre/2,0,max(yn)]);yn=angle(fft(yn2));

n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes18);stem(n1,yn,'.');

title('信號濾波后的相頻特性');xlabel('頻率（Hz）');

ylabel('相位');

axis([0,fre/2,min(yn),max(yn)]);

pushbutton1_Callback(hObject, eventdata, handles);

%---Executes during object creation, after setting all properties.function axes12_CreateFcn(hObject, eventdata, handles)% hObject handle to axes12(see GCBO)

% eventdata reservedhandles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate axes12

%---Executes on mouse press over axes background.function axes12_ButtonDownFcn(hObject, eventdata, handles)% hObject handle to axes12(see GCBO)

% eventdata reservedto be defined in a future version of MATLAB % handles structure with handles and user data(see GUIDATA)

% Hints: get(hObject,'String')returns contents of edit12 as text % str2double(get(hObject,'String'))returns contents of edit12 as a double

%---Executes during object creation, after setting all properties.function edit12_CreateFcn(hObject, eventdata, handles)% hObject handle to edit12(see GCBO)

% eventdata reservedhandles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.% See ISPC and COMPUTER.if ispc

&&

isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))set(hObject,'BackgroundColor','white');end

%---If Enable == 'on', executes on mouse press in 5 pixel border.%---Otherwise, executes on mouse press in 5 pixel border or over edit5.function edit5_ButtonDownFcn(hObject, eventdata, handles)% hObject handle to edit5(see GCBO)

% eventdata reservedto be defined in a future version of MATLAB % handles structure with handles and user data(see GUIDATA)

%---Executes during object deletion, before destroying properties.function edit12_DeleteFcn(hObject, eventdata, handles)% hObject handle to edit12(see GCBO)

% eventdata reserved-to be defined in a future version of MATLAB % handles structure with handles and user data(see GUIDATA)

第二篇：數字信號處理知識總結課案（模版）

1.傅里葉變換

x(n)的一個周期； 有限長序列

x(n)

可看成周期序列

~

x(n)看成x把 ~(n)的以N為周期的周期延拓。有限長序列的離散傅里葉變換（DFT）：

N?1?nk ?X(k)?DFT[x(n)]??x(n)WNx(n)?x((n))N?~?n?0 ??~1N?1?nk?x(n)?x(n)RN(n)?x((n))NRN(n)x(n)?IDFT[X(k)]?X(k)W ??N?Nk?0?

① 長度為N的有限長序列 x(n)，其離散傅里葉變換 X(k)仍是一個長度為N 的有限長序列；

② x(n)與X(k)是一個有限長序列離散傅里葉變換對，已知x(n)就能唯一地確定 X(k)；同樣已知X(k)也就唯一地確定x(n)。實際上x(n)與 X(k)都是長度為 N 的序列（復序列）都有N個獨立值，因而具有等量的信息； ③ 有限長序列隱含著周期性。

2.循環卷積（有可能會讓畫出卷積過程或結果）

循環卷積過程為：

最后結果為：

3.（見課本）

課本

3、線性卷積（有可能會讓畫出卷積過程或結果）

以下為PPT上的相關題目：

4.計算分段卷積：重疊相加法和重疊保留法（一定會考一種）

重疊相加法解題基本步驟:

將長序列均勻分段，每段長度為M；

基于DFT快速卷積法，通過循環卷積求每一段的線性卷積；

依次將相鄰兩段的卷積的N-1個重疊點相加，得到最終的卷積結果。

4.級聯、并聯、直接形（畫圖）以下為課后作業相關題目：

1.已知系統用下面差分方程描述: 311y(n)＝y(n?1)－y(n?2)＋x(n)?x(n?1)483

試分別畫出系統的直接型、級聯型和并聯型結構。式中x(n)和y(n)分別表示系統的輸入和輸出信號。

解： 將原式移項得

311y(n)?y(n?1)?y(n?2)?x(n)?x(n?1)483

將上式進行Z變換，得到

311?1?2?1Y(z)?Y(z)z?Y(z)z?X(z)?X(z)z483

11?z?13H(z)?311?z?1?z?248(1)按照系統函數

H(z)，根據Masson公式，畫出直接型結構如題1解圖

（一）所示。1?z?11(2)將H(z)的分母進行因式分解：

H(z)? ? 3111?z?11?z?1

241?11?z按照上式可以有兩種級聯型結構: 13H(z)? ? 111?z?11?z?1 241畫出級聯型結構如題1解圖

（二）（a)所示 1?z?11 H(z)? ? 3111?z?11?z?1

24畫出級聯型結構如題1解圖

（二）（b)所示

111?z?11?z?133H(z)??3?11?21?11?11?z?z(1?z)(1?z)4824

(3)

將H(z)進行部分分式展開：

1z?H(z)AB3???H(z)?11111?11?1z(z?)(z?)z?z?(1?z)(1?z)2424 2411?z?***zz?H(z)33?3???3?3H(z)?3111?11?111z(z?)z?1?z1?zz?z?24242410

根據上式畫出并聯型結構如題1解圖

（三）所示。

3.設系統的差分方程為

y(n)=(a+b)y(n－1)－aby(n－2)+x(n－2)+(a+b)x(n－1)+ab 式中, |a|<1，|b|<1，x(n)和y(n)分別表示系統的輸入和輸出信號, 試畫出系統的直接型和級聯型結構。

解：

(1)直接型結構。將差分方程進行Z變換，得到

Y(z)=(a+b)Y(z)z－1－abY(z)z－2+X(z)z－2－(a+b)X(z)z－1+ab Y(z)ab?(a?b)z?zH(z)??X(z)1?(a?b)z?1?abz?2按照Masson公式畫出直接型結構如題3解圖

（一）所示。

?1?2

(2)級聯型結構。將H(z)的分子和分母進行因式分解，得到

(a?z?1)(b?z?1)H(z)??H1(z)H2(z)?1?1(1?az)(1?bz)

按照上式可以有兩種級聯型結構:①

z?1?aH1(z)?1?az?1畫出級聯型結構如題3解圖

（二）(a)所示

z?1?bH2(z)?1?bz?1?1

z?aH1(z)??11?bz畫出級聯型結構如題3解圖

（二）(b)所示 ?1z?bH2(z)??11?az

四．設計模擬濾波器（考試時不能編代碼）一般步驟：

根據Ap、As、Ωs、Ωp，確定濾波器階次N和截止頻率Ωc。

P161 【例6.2.2】

設計一個模擬低通巴特沃斯濾波器，指標如下：

(1)通帶截止頻率：Ωp=0.2π；通帶最大衰減：Ap=7 dB。

(2)阻帶截止頻率：Ωs=0.3π；阻帶最小衰減：As=16dB。解： ?lg[(100.7?1)/(101.6?1)]? N?????2.79??32lg(0.2?/0.3?)??

由Ωp，得:

0.2?

Qc??0.49850.76 10?1由Ωs，得：

0.3?Q??0.5122 c1.6610?1

在上面兩個Ωc之間選Ωc=0.5。

最后可得（級聯型）：

0.125Ha(s)?

(s?0.5)(s2?0.5s?0.25)

五、脈沖響應不變法（P177 第6.3節）156-158頁

脈沖響應不變法的優點：

? 時域逼近。使數字濾波器的單位脈沖響應完全模仿模擬濾波器的單位沖激響應，即時域逼近良好。? 線性頻率關系。

模擬頻率Ω和數字頻率ω之間呈線性關系ω=ΩT。脈沖響應不變法的缺點： ? 混疊失真效應

因此，只適用于限帶的模擬濾波器（例如衰減特性很好的低通或帶通濾波器），而且高頻衰減越快，混疊效應越小；而對于高通和帶阻濾波器，由于它們在高頻部分不衰減，因此會產生混疊現象。

六、雙線性變換法

七，與實驗相關

本題中老師會給出類似于下列表達式的信號：

S(t)?2?3cos(2?50?t?30?/180)?1.5cos(2?75?t?90?/180)

要求用脈沖相應不變法或雙線性法編寫主要的代碼（如下面代碼）來達到濾除其中的部分信號，并畫出你所設計的濾波器的頻響曲線，并標明Ωs、Ωp，以及濾波后信號的時域波形（波形中要體現相位特征）。1）脈沖響應不變法濾除第三個信號： Fs=256;

% 采樣頻率 fp=60;

% 通帶截止頻率 fs=70;

% 阻帶截止頻率 Rp=1;

Rs=25;Wp=(fp/Fs)*2*pi;

%臨界頻率采用角頻率表示 Ws=(fs/Fs)*2*pi;

%臨界頻率采用角頻率表示 OmegaP=Wp*Fs;OmegaS=Ws*Fs;[n,Wc]=buttord(OmegaP,OmegaS,Rp,Rs,'s');[b,a]=butter(n,Wc,'s');[Bz,Az]=impinvar(b,a,Fs);

2）雙線性法濾除第三個信號： Fs=256;

% 采樣頻率 fp=60;

% 通帶截止頻率 fs=70;

% 阻帶截止頻率 Rp=1;

Rs=25;Wp=(fp/Fs)*2*pi;

% 臨界頻率采用角頻率表示 Ws=(fs/Fs)*2*pi;

% 臨界頻率采用角頻率表示 OmegaP=2*Fs*tan(Wp/2);

% 頻率預畸 OmegaS=2*Fs*tan(Ws/2);[n,Wc]=buttord(OmegaP,OmegaS,Rp,Rs,'s');[b,a]=butter(n,Wc,'s');

[Bz,Az]=bilinear(b,a,Fs);注：要好好看實驗中關于低通，高通，帶通，帶阻的設計代碼。

帶通：

fp1=40;

% 通帶截止頻率 fs1=30;

% 阻帶截止頻率 fp2=60;

% 通帶截止頻率 fs2=70;

% 阻帶截止頻率 Rp=1;

Rs=25;Wp1=(fp1/Fs)*2*pi;Ws1=(fs1/Fs)*2*pi;Wp2=(fp2/Fs)*2*pi;Ws2=(fs2/Fs)*2*pi;Wp=[Wp1,Wp2];

% 向量 Ws=[Ws1,Ws2];

% 向量

帶阻：

fp1=30;

% 通帶截止頻率 fs1=40;

% 阻帶截止頻率 fp2=70;

% 通帶截止頻率 fs2=60;

% 阻帶截止頻率 Rp=1;

Rs=25;Wp1=(fp1/Fs)*2*pi;Ws1=(fs1/Fs)*2*pi;Wp2=(fp2/Fs)*2*pi;Ws2=(fs2/Fs)*2*pi;Wp=[Wp1,Wp2];Ws=[Ws1,Ws2];若S(??5t)?信號表達

i(n?2

式為

t2??70???3c?to?st?s?則相關代碼為：

1)低通濾波器代碼

fp=110;% 通帶截止頻率 fs=130;% 阻帶截止頻率

Rp=1;Rs=25;Wp=(fp/Fs)*2*pi;Ws=(fs/Fs)*2*pi;%臨界頻率采用角頻率表示

(1):脈沖響應不變法

OmegaP=Wp*Fs;OmegaS=Ws*Fs;[n,Wc]=buttord(OmegaP,OmegaS,Rp,Rs,'s');[b,a]=butter(n,Wc,'s');% 指明為高通濾波器 [Bz,Az]=impinvar(b,a,Fs);

（2）雙線性變換法

OmegaP=2*Fs*tan(Wp/2);OmegaS=2*Fs*tan(Ws/2);% 頻率預畸 [n,Wc]=buttord(OmegaP,OmegaS,Rp,Rs,'s');[b,a]=butter(n,Wc,'s');[Bz,Az]=bilinear(b,a,Fs);

2)高通濾波器

fp=280;% 通帶截止頻率 fs=260;% 阻帶截止頻率 Rp=1;Rs=25;Wp=(fp/Fs)*2*pi;%臨界頻率采用角頻率表示 Ws=(fs/Fs)*2*pi;%臨界頻率采用角頻率表示(2):雙線性變換法

OmegaP=2*Fs*tan(Wp/2);% 頻率預畸 OmegaS=2*Fs*tan(Ws/2);[n,Wc]=buttord(OmegaP,OmegaS,Rp,Rs,'s');[b,a]=butter(2*n,Wc,'high','s');[Bz, Az]=bilinear(b,a,Fs);3)帶通濾波器

fp1=130;% 通帶截止頻率 fs1=110;% 阻帶截止頻率 fp2=255;% 通帶截止頻率 fs2=265;% 阻帶截止頻率 Rp=1;Rs=25;Wp1=(fp1/Fs)*2*pi;Ws1=(fs1/Fs)*2*pi;Wp2=(fp2/Fs)*2*pi;Ws2=(fs2/Fs)*2*pi;Wp=[Wp1,Wp2];Ws=[Ws1,Ws2];(2):雙線性變換法

OmegaP=2*Fs*tan(Wp/2);% 頻率預畸 OmegaS=2*Fs*tan(Ws/2);[n,Wc]=buttord(OmegaP,OmegaS,Rp,Rs,'s');[b,a]=butter(2*n,Wc,'s');[Bz,Az]=bilinear(b,a,Fs);4)帶阻濾波器的代碼如下：

fp1=110;% 通帶截止頻率 fs1=240;% 阻帶截止頻率 fp2=265;% 通帶截止頻率 fs2=255;% 阻帶截止頻率 Rp=1;Rs=25;Wp1=(fp1/Fs)*2*pi;Ws1=(fs1/Fs)*2*pi;Wp2=(fp2/Fs)*2*pi;Ws2=(fs2/Fs)*2*pi;

Wp=[Wp1,Wp2];Ws=[Ws1,Ws2];(2):雙線性變換法

OmegaP=2*Fs*tan(Wp/2);% 頻率預畸 OmegaS=2*Fs*tan(Ws/2);[n,Wc]=buttord(OmegaP,OmegaS,Rp,Rs,'s');[b,a]=butter(2*n,Wc,'stop','s');[Bz,Az]=bilinear(b,a,Fs);

第三篇：基于微課的數字信號處理教學方法的探討

基于微課的數字信號處理教學方法的探討

【摘 要】基于目前數字信號處理課程在教學方法及教學手段上的不足，本文提出了基于“微課”的教學方法，通過對該課程重點難點內容做成微課視頻的形式，供學生課前或課后學習，以提高學生的學習興趣及自主學習能力，再結合雙向互動式研究型教學方法及教學實踐力度的加大，從而改善教學效果，使獨立院校學生能更好的掌握這門課程。

【關鍵詞】數字信號處理；微課；教學改革；實踐教學

0 引言

數字信號處理作為通信、電子信息類專業的一門專業基礎課程，要求學生有微積分、信號與系統分析等課程的基礎，另外數字信號處理是一門理論性非常強的課程，抽象概念較多，涉及的數學知識也較多，而獨立學院學生數學基礎、信號與系統分析基礎都相對比較薄弱，所以學習本課程有一定的難度，易出現學生學習興趣低、怕學以及厭學等現象。針對該問題，本文提出了一種基于微課的教學方式，將課程內容制作成微課，分享給學生，以提高學生的學習興趣。

“微課”采用了精細化管理的教學模式，把“慕課”平臺上的“高級軟件工程”課程加以優化及改編之后得到“微課”平臺，加州大學伯克利分校的 Armando Fox 教授將“微課”平臺引入了校園，并取得了意想不到的好效果。2013 年，哈佛大學也開始嘗試“微課”平臺的教學，效果很好，且學術委員會主席Robert Lue 教授對“微課”給出了極高的評價，認為相對于“慕課”教學方式而言，“微課”體現了“幾乎不可避免的進化”。科羅拉多大學全球分校最初采用“慕課”大班教學，經過調研，對學生的需求進行分析，最終也將所有課程都改成采用“微課”的在線授課模式。

從人的認知規律出發，人在學習的時候前十分鐘注意力最集中，所以“微課”長度一般設置在5-15分鐘的?！拔⒄n”在闡述問題的方式上除了教師講解外，還結合動畫、視頻等各種信息技術，使得對知識的講解過程更為直觀。且由于一個微課視頻時間較短，可以突破傳統課堂對時間和地點的限制，學生可以隨時隨地學習，也可以根據個人學習情況調整學習進度，使得學生真正成為學習的主角。數字信號處理教學中出現的問題

近幾年來，數字信號處理課程多媒體與黑板板書相結合的教學方式，再結合Matlab軟件實驗，教學效果有一定的改善，但仍存在不足之處，主要體現為一下幾個方面：

1.1 教學內容方面

公式推導比較多，概念比較抽象，內容枯燥，學生無法將所學的內容與工程實踐相結合，再者獨立學院學生基礎較差，所以易出現“太難，聽不懂”等不學、厭學現象。

1.2 教學方法方面

傳統教學多采用“灌輸式”的方法進行授課，教師作為課堂主體，單向地向學生灌輸知識，教學的方式過于單調，不足以吸引學生的注意力，和學習興趣的提高。

1.3 實踐教學方面

實驗目前僅局限于理論仿真，而沒有硬件的實踐操作，學生很難把所學的理論知識和理論仿真同實際應用聯系起來，而獨立學院學生更應該把學習的重心往實踐方向傾斜。數字信號處理教學改革措施

針對數字信號處理教學過程中出現的3個問題，我們可以從以下3個方面進行改進。

2.1 教學內容的改革

數字信號處理課程的知識結構可分為變換域和數字濾波器兩部分。對于變換域部分的知識，如離散時間信號的時域頻域分析部分與信號與系統分析Z變換部分是重復的，所以在實際教學中，可以對這部分內容進行精簡和提煉，并且刪除一些復雜的公式推導。關于時域與變換域之間的關系、波形變化以及濾波器的原理等抽象的內容，設計原理的講解可以精簡，而應側重如何通過Matlab實現具體濾波器的設計，另外，在課堂中增加數字信號應用實例的MALTAB程序講解，演示仿真波形，通過多波形的分析，使學生理解時域信號好頻域信號的特性，從而進一步加深學生對理論知識的理解。

2.2 教學方法和教學手段的改革

2.2.1 引入微課形式的教學方法

在傳統“灌輸式”、PPT講解的基礎上，借助視頻學習、慕課、微課豐富教學手段，將課程重要內容制作成微課，在課堂上播放，吸引學生注意力，提高學生學習興趣；也可放到網絡平臺上，供學生學習，充分利用學生的碎片化時間提供在線學習的機會。甚至能夠翻轉課堂，學生可在課前通過微課自主完成課堂學習，而在課堂上對重點難點進行討論，使學生在討論的過程中深化對知識點的理解，同時可以活躍課堂氣氛，讓學生成為課堂的主體。

2.2.2 雙向互動式研究型學習

針對課程部分重點難點教學內容，提取出合適的討論課題，進行定期或不定期的專題討論。首先，通過課堂或是其他網絡平臺向學生發布題目，讓學生提前做好討論或者報告準備，采用多媒體上講臺進行專題討論或報告。學生通過對課題的準備及討論過程可對該課程內容有更深入的理解，且能培養學生的辯證思考的能力，另外教師對課題討論的主要參與人做相應的獎勵，比如加平時分，來誘導學生自主學習。

2.3 實踐教學的改革

部分獨立學院目前已經開展了數字信號處理的軟件實驗課程，即基于Matlab軟件，通過編程實現對信號的仿真，可一定程度上加深學生對信號時域頻域的形象理解，可提高學生的編程能力，但畢竟Matlab軟件實驗只是做信號建模及仿真，所以很多學生仍然存在“學了這門課，卻不知道怎么在實際中應用”的疑問，因此，增加DSP硬件實驗是必要的。然而增加DSP硬件實驗需要校方新建DSP實驗室，該問題不是朝夕之間能解決的問題。所以針對此問題，可先進行探索性實驗，即通過收集整理其他院校的DSP硬件實現演示視頻，使學生對理論的仿真到實際的硬件實現有一個形象全面的了解，比如較抽象的知識點像卷積運算、FIR/IIR數字濾波器以及快速傅里葉變換FFT的DSP實現，通過實驗課使學生把抽象的知識與實際應用相結合?？偨Y

實踐證明，基于微課及雙向互動式的教學方法能充分調動學生的學習積極性，提高學生的自學能力。另外，引入DSP實踐教學，使學生更能將理論知識與工程實踐相結合。

【參考文獻】

[1]曾偉梁，劉穎，鮑曼，徐宏艷，趙沛豐.用 MATLAB 實踐數字信號處理課程的輔助教學[J].黑龍江生態工程職業學院學報，2014（1）.[2]高靜，王鳳文，舒冬梅.《數字信號處理》課程教學實踐與探索[J].科技創新導報，2011（4）：153-154.[3]曾偉梁，劉穎，鮑曼，徐宏艷，趙沛豐.用MATLAB 實踐數字信號處理課程的輔助教學[J].黑龍江生態工程職業學院學報，2014（1）.[4]劉伯紅，閻英，方義秋.高校計算機實踐教學質量保障體系改革探索與實踐[J].實驗室研究與探索，2012，31（12）：121-123，150.[責任編輯：楊玉潔]

第四篇：數字信號處理課程設計..

課程設計報告

課程名稱： 數字信號處理 課題名稱： 語音信號的處理與濾波

姓 名： 學 號： 院 系： 專業班級： 指導教師： 完成日期： 2013年7月2日

目錄

第1部分 課程設計報告………………………………………3 一．設計目的……………………………………………3 二．設計內容……………………………………………3 三．設計原理……………………………………………3 四．具體實現……………………………………………5 1.錄制一段聲音…………………………………5 2.巴特沃斯濾波器的設計………………………8 3.將聲音信號送入濾波器濾波…………………13 4.語音信號的回放………………………………19 5.男女語音信號的頻譜分析……………………19 6.噪聲的疊加和濾除……………………………22 五． 結果分析……………………………………………27 第2部分 課程設計總結………………………………28 一． 參考文獻……………………………………………28

第1部分 課程設計報告

一．設計目的

綜合運用本課程的理論知識進行頻譜分析以及濾波器設計，通過理論推導得出相應結論，并利用MATLAB作為工具進行實現，從而復習鞏固課堂所學的理論知識，提高對所學知識的綜合應用能力，并從實踐上初步實現對數字信號的處理。

二．設計內容

錄制一段個人自己的語音信號，并對錄制的信號進行采樣；畫出采樣后語音信號的時域波形和頻譜圖；給定濾波器的性能指標，采用窗函數法和雙線性變換法設計濾波器，并畫出濾波器的頻率響應；然后用自己設計的濾波器對采集的信號進行濾波，畫出濾波后信號的時域波形和頻譜，并對濾波前后的信號進行對比，分析信號的變化；回放語音信號；換一個與你性別相異的人錄制同樣一段語音內容，分析兩段內容相同的語音信號頻譜之間有什么特點；再錄制一段同樣長時間的背景噪聲疊加到你的語音信號中，分析疊加前后信號頻譜的變化，設計一個合適的濾波器，能夠把該噪聲濾除；

三．設計原理

1.在Matlab軟件平臺下，利用函數wavrecord()，wavwrite(),wavread(),wavplay()對語音信號進行錄制，存儲，讀取，回放。

2.用y=fft(x)對采集的信號做快速傅立葉變換，并用[h1,w]=freqz(h)進行DTFT變換。

3.掌握FIR DF線性相位的概念，即線性相位對h(n)、H(?)及零點的約束，了解四種FIR DF的頻響特點。

4.在Matlab中，FIR濾波器利用函數fftfilt對信號進行濾波。

5.抽樣定理

連續信號經理想抽樣后時域、頻域發生的變化（理想抽樣信號與連續信號頻譜之間的關系）

理想抽樣信號能否代表原始信號、如何不失真地還原信號即由離散信號恢復連續信號的條件（抽樣定理）

理想采樣過程描述： 時域描述：

?a(t)?xa(t)?T(t)??xa(t)?(t?nT)??xa(nT)?(t?nT)xn???n??????T(t)?頻域描述：利用傅氏變換的性質，時域相乘頻域卷積，若

n?????(t?nT)??a(t)Xa(j?)?xXa(j?)?xa(t)?T(j?)??T(t)

則有

?(j?)?1X(j?)??(j?)XaaT2?1?2?1??Xa(j?)??Xa(j??jk)??Xa(j??jk?s)Tk???TTk????(j?)與X(j?)的關系：理想抽樣信號的頻譜是連續信號頻譜的Xaa

周期延拓,重復周期為?s(采樣角頻率)。如果：

?X(j?)?Xa(j?)??a??0???s/2???s/2即連續信號是帶限的，且信號最高頻率不超過抽樣頻率的二分之一，則可不失真恢復。

奈奎斯特采樣定理：要使實信號采樣后能夠不失真還原，采樣頻率必須大于信號最高頻率的兩倍：?s?2?h 或 fs?2fh

四．具體實現

1.錄制一段聲音

1.1錄制并分析

在MATLAB中用wavrecord、wavread、wavplay、wavwrite對聲音進行錄制、讀取、回放、存儲。

程序如下：

Fs=8000;%抽樣頻率 time=3;%錄音時間 fprintf('按Enter鍵錄音%ds',time);%文字提示 pause;%暫停命令 fprintf('錄音中......');x=wavrecord(time*Fs,Fs,'double');%錄制語音信號 fprintf('錄音結束');%文字提示 fprintf('按Enter鍵回放錄音');pause;%暫停命令

wavplay(x,Fs);%按任意鍵播放語音信號

wavwrite(x,Fs,'C:UsersacerDesktop數字信號sound.wav');%存儲語音信號

N=length(x);%返回采樣點數 df=fs/N;%采樣間隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%頻帶寬度 figure(2);subplot(2,1,1);plot(x);%錄制信號的時域波形 title('原始信號的時域波形');%加標題 ylabel('幅值/A');%顯示縱坐標的表示意義 grid;%加網格

y0=fft(x);%快速傅立葉變換 figure(2);subplot(2,1,2);plot(f,abs(y0(n1)));%原始信號的頻譜圖 title('原始信號的頻譜圖');%加標題 xlabel('頻率w/pi');%顯示橫坐標表示的意義 ylabel('幅值 ');%顯示縱坐標表示的意義 title('原始信號的頻譜圖');%加標題

grid;%加網格

圖1.1 原始信號的時域與頻譜圖

1.2濾除無效點

針對實際發出聲音落后錄制動作半拍的現象，如何拔除對無效點的采樣的問題: 出現這種現象的原因主要是錄音開始時，人的反應慢了半拍，導致出現了一些無效點，而后而出現的無效的點，主要是已經沒有聲音的動作，先讀取聲音出來，將原始語音信號時域波形圖畫出來，根據己得到的信號，可以在第二次讀取聲音的后面設定采樣點，取好有效點，畫出濾除無效點后的語音信號時域波形圖，對比可以看出。這樣就可以解決這個問題。

x=wavread('C:UsersacerDesktop數字信號sound.wav', 7

[4000,24000]);%從4000點截取到24000結束 plot(x);%畫出截取后的時域圖形 title('截取后的聲音時域圖形');%標題 xlabel('頻率');ylabel('振幅');grid;%畫網格

圖1.2 去除無效點

2.巴特沃斯濾波器的設計

2.1設計巴特沃思低通濾波器

MATLAB程序如下。濾波器圖如圖3.3所示。

%低通濾波

fp=1000;fs=1200;Fs=22050;rp=1;rs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'s');[Bz,Az]=bilinear(B,A,Fs1);figure(1);[h,w]=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h));title('巴特沃斯低通濾波器');xlabel('頻率（HZ）');ylabel('耗損（dB）');gridon;9

圖2.1 巴特沃思低通濾波器

2.2設計巴特沃思高通濾波器

MATLAB程序如下。濾波器圖如圖3.5所示。%高通濾波

fp=4800;fs=5000;Fs=22050;rp=1;rs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;T=1;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);10

[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'high','s');[Bz,Az]=bilinear(B,A,Fs1);figure(1);[h,w]=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h));title('巴特沃斯高通濾波器');xlabel('頻率（HZ）');ylabel('耗損（dB）');grid on;

圖2.2巴特沃思高通濾波器

2.3設計巴特沃思帶通濾波器

MATLAB程序如下。濾波器圖如圖3.7所示。%帶通濾波

fp=[1200,3000];fs=[1000,3200];Fs=8000;rp=1;rs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;T=1;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'s');[Bz,Az]=bilinear(B,A,Fs1);figure(4);[h,w]=freqz(Bz,Az,512,Fs1*1000);plot(w,abs(h));title('巴特沃斯帶通濾波器');xlabel('頻率（HZ）');ylabel('耗損（dB）');grid on;12

圖2.3巴特沃思帶通濾波器

3．將聲音信號送入濾波器濾波

x=wavread('C:UsersacerDesktop數字信號sound.wav');%播放原始信號

wavplay(x,fs);%播放原始信號 N=length(x);%返回采樣點數 df=fs/N;%采樣間隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%頻帶寬度 figure(4);subplot(4,2,1);plot(x);%錄制信號的時域波形

title('原始信號的時域波形');%加標題 ylabel('幅值/A');%顯示縱坐標的表示意義 grid;%加網格

y0=fft(x);%快速傅立葉變換 subplot(4,2,3);plot(f,abs(y0(n1)));%原始信號的頻譜圖 title('原始信號的頻譜圖');%加標題 xlabel('頻率w/pi');%顯示橫坐標表示的意義 ylabel('幅值 ');%顯示縱坐標表示的意義 title('原始信號的頻譜圖');%加標題 grid;%加網格

3．1低通濾波器濾波 fs=8000;beta=10.056;wc=2*pi*1000/fs;ws=2*pi*1200/fs;width=ws-wc;wn=(ws+wc)/2;n=ceil(12.8*pi /width);h=fir1(n,wn/pi,'band',kaiser(n+1,beta));[h1,w]=freqz(h);

ys=fftfilt(h,x);%信號送入濾波器濾波，ys為輸出 fftwave=fft(ys);%將濾波后的語音信號進行快速傅立葉變換 figure(4);subplot(4,2,2);%在四行兩列的第二個窗口顯示圖形 plot(ys);%信號的時域波形

title('低通濾波后信號的時域波形');%加標題 xlabel('頻率w/pi');ylabel('幅值/A');%顯示標表示的意義 grid;%網格

subplot(4,2,4);%在四行兩列的第四個窗口顯示圖形 plot(f, abs(fftwave(n1)));%繪制模值 xlabel('頻率w/pi');ylabel('幅值/A');%顯示標表示的意義

title('低通濾波器濾波后信號的頻譜圖');%標題 grid;%加網格

wavplay(ys,8000);%播放濾波后信號

3．2高通濾波器濾波 fs=8000;beta=10.056;ws=2*5000/fs;wc=2*4800/fs;

width=ws-wc;wn=(ws+wc)/2;n=ceil(12.8*pi/width);h=fir1(n,wn/pi, 'high',kaiser(n+2,beta));[h1,w]=freqz(h);ys=fftfilt(h,x);%將信號送入高通濾波器濾波 subplot(4,2,5);%在四行兩列的第五個窗口顯示圖形 plot(ys);%信號的時域波形 xlabel('頻率w/pi');ylabel('幅值/A');%顯示標表示的意義 title('高通濾波后信號的時域波形');%標題 ylabel('幅值/A');%顯示縱坐標的表示意義 grid;%網格

fftwave=fft(ys);%將濾波后的語音信號進行快速傅立葉變換 subplot(4,2,7);%在四行兩列的第七個窗口顯示圖形 plot(f,abs(fftwave(n1)));%繪制模值 axis([0 1 0 50]);xlabel('頻率w/pi');ylabel('幅值/A');%顯示標表示的意義

title('高通濾波器濾波后信號的頻譜圖');%標題 grid;%加網格

wavplay(ys,8000);%播放濾波后信號

3．3帶通濾波器 fs=8000;beta=10.056;wc1=2*pi*1000/fs;wc2=2*pi*3200/fs;ws1=2*pi*1200/fs;ws2=2*pi*3000/fs;width=ws1-wc1;wn1=(ws1+wc1)/2;wn2=(ws2+wc2)/2;wn=[wn1 wn2];n=ceil(12.8/width*pi);h=fir1(n,wn/pi,'band',kaiser(n+1,beta));[h1,w]=freqz(h);ys1= fftfilt(h,x);%將信號送入高通濾波器濾波 figure(4);subplot(4,2,6);%在四行兩列的第六個窗口顯示圖形 plot(ys1);%繪制后信號的時域的圖形 title('帶通濾波后信號的時域波形');%加標題 xlabel('頻率w/pi');ylabel('幅值/A');%顯示縱坐標表示的意義 grid;%網格

fftwave=fft(ys1);%對濾波后的信號進行快速傅立葉變換 subplot(4,2,8);%在四行兩列的第八個窗口顯示圖形

plot(f, abs(fftwave(n1)));%繪制模值 axis([0 1 0 50]);xlabel('頻率w/pi');ylabel('幅值/A');%顯示標表示的意義 title('帶通濾波器濾波后信號的頻譜圖');%加標題 grid;%網格

wavplay(ys1,8000);%播放濾波后信號 圖形如下：

原始信號的時域波形幅值/A0-1012x 10原始信號的頻譜圖34幅值/A1低通濾波后信號的時域波形0.50-0.5012頻率w/pi3400.51頻率w/pi高通濾波后信號的時域波形幅值/A0幅值/A0幅值/Ax 10高通濾波器濾波后信號的頻譜圖5012頻率w/pi34幅值/A0.20-0.2幅值/A2001000x 10低通濾波器濾波后信號的頻譜圖200100000.51頻率w/pi帶通濾波后信號的時域波形0.50-0.501234頻率w/pix 10帶通濾波器濾波后信號的頻譜圖50幅值 00.5頻率w/pi1000.5頻率w/pi1

分析：三個濾波器濾波后的聲音與原來的聲音都發生了變化。其中低

通的濾波后與原來聲音沒有很大的變化，其它兩個都又明顯的變化

4.語音信號的回放

sound(xlow,Fs,bits);%在Matlab中，函數sound可以對聲音進行回放，其調用格式: sound(xhigh, Fs,bits);%sound(x, Fs, bits);sound(xdaitong, Fs,bits);5.男女語音信號的頻譜分析

5.1 錄制一段異性的聲音進行頻譜分析

Fs=8000;%抽樣頻率 time=3;%錄音時間 fprintf('按Enter鍵錄音%ds',time);%文字提示 pause;%暫停命令 fprintf('錄音中......');x=wavrecord(time*Fs,Fs,'double');%錄制語音信號 fprintf('錄音結束');%文字提示 fprintf('按Enter鍵回放錄音');pause;%暫停命令 wavplay(x,Fs);%按任意鍵播放語音信號

wavwrite(x,Fs,'C:UsersacerDesktop數字信號sound2.wav');%存儲語音信號

5.2 分析男女聲音的頻譜

x=wavread(' C:UsersacerDesktop數字信號sound2.wav ');%播放原始信號,解決落后半拍

wavplay(x,fs);%播放原始信號 N=length(x);%返回采樣點數 df=fs/N;%采樣間隔 n1=1:N/2;

f=[(n1-1)*(2*pi/N)]/pi;%頻帶寬度 figure(1);subplot(2,2,1);plot(x);%錄制信號的時域波形

title('原始女生信號的時域波形');%加標題 ylabel('幅值/A');%顯示縱坐標的表示意義 grid;%加網格

y0=fft(x);%快速傅立葉變換 subplot(2,2,2);plot(f,abs(y0(n1)));%原始信號的頻譜圖 title('原始女生信號的頻譜圖');%加標題 xlabel('頻率w/pi');%顯示橫坐標表示的意義 ylabel('幅值 ');%顯示縱坐標表示的意義 grid;%加網格

[y,fs,bits]=wavread(' C:UsersacerDesktop數字信號sound.wav ');% 對語音信號進行采樣

wavplay(y,fs);%播放原始信號 N=length(y);%返回采樣點數 df=fs/N;%采樣間隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%頻帶寬度 subplot(2,2,3);plot(y);%錄制信號的時域波形

title('原始男生信號的時域波形');%加標題 ylabel('幅值/A');%顯示縱坐標的表示意義 grid;%加網格

y0=fft(y);%快速傅立葉變換

subplot(2,2,4);%在四行兩列的第三個窗口顯示圖形 plot(f,abs(y0(n1)));%原始信號的頻譜圖 title('原始男生信號的頻譜圖');%加標題 xlabel('頻率w/pi');%顯示橫坐標表示的意義 ylabel('幅值 ');%顯示縱坐標表示的意義 grid;%加網格

5.3男女聲音的頻譜圖

原始女生信號的時域波形0.50-0.5-1150100原始女生信號的頻譜圖幅值/A幅值 012345000x 10原始男生信號的時域波形0.50.5頻率w/pi原始男生信號的頻譜圖1300200幅值/A0幅值 012x 1034100-0.5000.5頻率w/pi1

圖5.3男女聲音信號波形與頻譜對比

分析：就時域圖看，男生的時域圖中振幅比女生的高，對于頻譜圖女生的高頻成分比較多

6.噪聲的疊加和濾除

6.1錄制一段背景噪聲

Fs=8000;%抽樣頻率 time=3;%錄音時間 fprintf('按Enter鍵錄音%ds',time);%文字提示 pause;%暫停命令 fprintf('錄音中......');x=wavrecord(time*Fs,Fs,'double');%錄制語音信號

fprintf('錄音結束');%文字提示 fprintf('按Enter鍵回放錄音');pause;%暫停命令 wavplay(x,Fs);%按任意鍵播放語音信號 wavwrite(x,Fs,'C:UsersacerDesktop數字信號噪音.wav');%存儲語音信號

6.2 對噪聲進行頻譜的分析

[x1,fs,bits]=wavread(' C:UsersacerDesktop數字信號噪音.wav ');%對語音信號進行采樣

wavplay(x1,fs);%播放噪聲信號 N=length(x1);%返回采樣點數 df=fs/N;%采樣間隔

n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%頻帶寬度 figure(5);subplot(3,2,1);plot(x1);%信號的時域波形 title('噪聲信號的時域波形');grid;ylabel('幅值/A');y0=fft(x1);%快速傅立葉變換

subplot(3,2,2);plot(f,abs(y0(n1)));%噪聲信號的頻譜圖 ylabel('幅值');title('噪聲信號的頻譜圖');

6.3原始信號與噪音的疊加

fs=8000;[x,fs,bits]=wavread(' C:UsersacerDesktop數字信號sound.wav ');%對錄入信號進行采樣

[x1,fs,bits]=wavread(' C:UsersacerDesktop數字信號噪音.wav ');%對噪聲信號進行采樣

yy=x+x1;%將兩個聲音疊加

6.4疊加信號的頻譜分析：

wavplay(yy,fs);%播放疊加后信號 N=length(yy);%返回采樣點數 df=fs/N;%采樣間隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%頻帶寬度 figure(5);subplot(3,2,3);plot(yy,'LineWidth',2);%信號的時域波形

title('疊加信號的時域波形');xlabel('時間/t');ylabel('幅值/A');grid;y0=fft(yy);%快速傅立葉變換 subplot(3,2,4);plot(f,abs(y0(n1)));%疊加信號的頻譜圖 title('疊加信號的頻譜圖');xlabel('頻率w/pi');ylabel('幅值/db');grid;

6.5 設計一個合適的濾波器將噪聲濾除 fs=18000;%采樣頻率 Wp=2*1000/fs;%通帶截至頻率 Ws=2*2000/fs;%阻帶截至頻率 Rp=1;%最大衰減 Rs=100;%最小衰減

[N,Wn]=buttord(Wp,Ws,Rp,Rs);%buttord函數（n為階數，Wn為截至頻率）

[num,den]=butter(N,Wn);%butter函數（num為分子系數den為分母系數）

[h,w]=freqz(num,den);%DTFT變換

ys=filter(num,den,yy);%信號送入濾波器濾波，ys為輸出 fftwave=fft(ys);%將濾波后的語音信號進行快速傅立葉變換 figure(5);subplot(3,2,5);plot(ys);%信號的時域波形

title('低通濾波后信號的時域波形');%加標題 ylabel('幅值/A');%顯示標表示的意義 grid;%網格 subplot(3,2,6);plot(f, abs(fftwave(n1)));%繪制模值 title('低通濾波器濾波后信號的頻譜圖');%標題 xlabel('頻率w/pi');ylabel('幅值/A');%顯示標表示的意義 grid;%加網格

wavplay(ys,8000);%播放濾波后信號 grid;圖形如下：

噪聲信號的時域波形1100噪聲信號的頻譜圖幅值/A0-1幅值0123450000.5疊加信號的頻譜圖1x 10疊加信號的時域波形10-101時間/t2200幅值/db34幅值/A100000.5頻率w/pi1x 10低通濾波后信號的時域波形0.5低通濾波器濾波后信號的頻譜圖200幅值/A0-0.5幅值/A012x 1034100000.5頻率w/pi1

圖6.1噪音的疊加與濾除前后頻譜對比

7.結果分析

1.錄制剛開始時，常會出現實際發出聲音落后錄制動作半拍，可在[x,fs,bits]=wavread('d:matlavworkwomamaaiwo.wav')加 窗[x,fs,bits]=wavread('d:matlavworkwomamaaiwo.wav'，[100 10000])，窗的長度可根據需要定義。

2.語音信號通過低通濾波器后，把高頻濾除，聲音變得比較低沉。當通過高通濾波器后，把低頻濾除，聲音變得比較就尖銳。通過帶通濾波器后，聲音比較適中。

3.通過觀察男生和女生圖像知：時域圖的振幅大小與性別無關，只與說話人音量大小有關，音量越大，振幅越大。頻率圖中，女生高 27

頻成分較多。

4.疊加噪聲后，噪聲與原信號明顯區分,但通過低通濾波器后，噪聲沒有濾除，信號產生失真。原因可能為噪聲與信號頻率相近無法濾除。

第2部分 課程設計總結

通過本次課程設計，使我們對數字信號處理相關知識有了更深刻的理解，尤其是對各種濾波器的設計。在設計的過程中遇到了很多問題，剛剛開始時曾天真的認為只要把以前的程序改了參數就可以用了，可是問題沒有我想象中的那么簡單，單純的搬程序是不能解決問題的。通過查閱資料和請教同學收獲了很多以前不懂的理論知識。再利用所學的操作，發現所寫的程序還是沒有能夠運行，通過不斷地調試，運行，最終得出了需要的結果。整個過程中學到了很多新的知識，特別是對Matlab的使用終于有些了解。在以后的學習中還需要深入了解這方面的內容。在這次的課程設計中讓我體會最深的是：知識來不得半點的馬虎。也認識到自己的不足，以后要進一步學習。

八．參考文獻

[1]數字信號處理教程(第三版)程佩青 清華大學出版社 [2]MATLAB信號處理 劉波 文忠 電子工業出版社 [3]MATLAB7.1及其在信號處理中的應用 王宏 清華大學出版社

[4]MATLAB基礎與編程入門 張威 西安電子科技大學出版社

[5] 數字信號處理及其MATLAB實驗 趙紅怡 張常 化學工業出版社

[6]MATLAB信號處理詳解 陳亞勇等 人民郵電出版社 [7] 數字信號處理

錢同惠 機械工業出版社 29

第五篇：數字信號處理實驗報告

JIANGSU

UNIVERSITY OF TECHNOLOGY

數字信號處理實驗報告

學院名稱： 電氣信息工程學院

專 業：

班 級： 姓 名： 學 號： 指導老師： 張維璽（教授）

2013年12月20日

實驗一 離散時間信號的產生

一、實驗目的

數字信號處理系統中的信號都是以離散時間形態存在的，所以對離散時間信號的研究是數字信號的基本所在。而要研究離散時間信號，首先需要產生出各種離散時間信號。使用MATLAB軟件可以很方便地產生各種常見的離散時間信號，而且它還具有強大繪圖功能，便于用戶直觀地處理輸出結果。

通過本實驗，學生將學習如何用MATLAB產生一些常見的離散時間信號，實現信號的卷積運算，并通過MATLAB中的繪圖工具對產生的信號進行觀察，加深對常用離散信號和信號卷積和運算的理解。

二、實驗原理

離散時間信號是指在離散時刻才有定義的信號，簡稱離散信號，或者序列。離散序列通常用x（n）來表示，自變量必須是整數。常見的離散信號如下：（1）單位沖激序列δ（n）

如果δ（n）在時間軸上延遲了k個單位，得到δ（n-k），即長度為N的單位沖激序列δ（n）可以通過下面的MATLAB命令獲得。

n=-（N-1）：N-1 x=[zeros(1,N-1)1 zeros(1,N-1)]； stem(n,x)延遲K個采樣點的長度為N的單位沖激序列δ（n-k）（k

n=0：N-1 y=[zeros(1,M)1 zeros(1,N-M-1)]； stem（n，y）

（2）單位階躍序列u（n）

如果u（n）在時間軸上延遲了k個單位，得到u（n-k），即長度為N的單位階躍序列u（n）可以通過下面的MATLAB命令獲得。

n=-（N-1）：N-1 x=[zeros(1,N-1)ones(1,N)]； stem(n,x)延遲的單位階躍序列可以使用類似于單位沖激序列的方法獲得。（3）矩形序列

矩形序列有一個重要的參數，就是序列的寬度N。矩形序列與u（n）之間的關系為矩形序列等= u（n）— u（n-N）。

因此，用MATLAB表示矩形序列可利用上面的單位階躍序列組合而成。（4）正弦序列x（n）

這里，正弦序列的參數都是實數。與連續的正弦信號不同，正弦序列的自變量n必須為整數。可以證明，只有當2π/w為有理數時，正弦序列具有周期性。

長度為N的正弦序列x（n）可以通過下面的MATLAB命令獲得。n=0：N-1 x=A*cos（2*pi*f*n/Fs+phase）（5）單邊實指數序列x（n）

長度為N的實指數序列x（n）可以通過下面的MATLAB命令實現。n=0：N-1 x=a.^n stem(n,x)單邊指數序列n的取值范圍為n>=0。當|a|>1時，單邊指數序列發散；當|a|<1時，單邊指數序列收斂。當a>0時，該序列均取正值；當a<0時，序列在正負擺動。

（6）負指數序列x（n）

當a=0時，得到虛指數序列x（n）。

與連續負指數信號一樣，我們將負指數序列實部和虛部的波形分開討論，得到如下結論：

1）當a>0時，負指數序列x（n）的實部和虛部分別是按指數規律增長的正弦振蕩序列；

2)當a<0時，負指數序列x（n）的實部和虛部分別是按指數規律衰減的正弦振蕩序列；

3)當a=0時，負指數序列x（n）即為虛指數序列，其實部和虛部分別是等幅的正弦振蕩序列；

長度為N的實指數序列x（n）可以通過下面的MATLAB命令實現。n=0：N-1 x=exp（（a.+j*w）*n）stem(n,real(x))或

stem(n,imag(x))

三、實驗內容及分析

?1n?01、編制程序產生單位沖激序列??n???“?并繪出其圖及??n?”學號后兩位0n?0?形。程序:(1)N=4;

n=-(N-1):N-1;

x=[zeros(1,N-1)1 zeros(1,N-1)];stem(n,x);

title('單位沖激序列');

grid on;

(2)N=6;

M=1;%學號01 n=-(N-1):N-1;

y=[zeros(1,N-M+1)1 zeros(1,N-M-1)];stem(n,y);

title('單位沖激序列');grid on;

分析:在上圖的基礎上向右平移了1個單位。

?1n?02、編制程序產生單位階躍序列u?n???、u?n?“學號后兩位”?及

0n?0?u?n??u?n?“學號后兩位”?，并繪出其圖形。程序: 4

(1)N=5;

n=-(N-1):N-1;

x=[zeros(1,N-1)ones(1,N)];stem(n,x);

title('單位階躍序列');grid on;

(2)N=6;

M=1;%學號01 n=-(N-1):N-1;

x=[zeros(1,N-M+1)ones(1,N-M)];stem(n,x);

title('單位階躍序列');grid on;

分析:在上圖的基礎上平移了1個單位.(3)N=6;

M=1;%學號01 n=-(N-1):N-1;

x=[zeros(1,N-1)ones(1,N)];y=[zeros(1,N-M+1)ones(1,N-M)];z=x-y;stem(n,z);

title('單位階躍序列');grid on;

2??

3、編制程序產生正弦序列x?n??cos?2?n?、x?n??cos??n?及

?學號后兩位?x?n??sin?2n?并繪出其圖形。

程序:(1)N=5;

A=1;

w=2*pi;phi=0;n=0:0.05:N-1;x=A*cos(w*n+phi);stem(n,x);title('余弦信號');grid on;

分析:該序列具有周期性,且輸出為余弦信號.(2)N=5;

A=1;

w=2*pi/1;%學號01 phi=0;n=0:0.05:N-1;x=A*cos(w*n+phi);stem(n,x);title('余弦信號');grid on;

;

分析:該序列具有周期性,且輸出為余弦信號.(3)N=5;

A=1;

w=2*pi;phi=0;

n=0:0.05:N-1;x=A*sin(w*n+phi);stem(n,x);title('正弦信號');grid on;

分析:該序列具有周期性,且輸出為正弦信號.4、編制程序產生復正弦序列x?n??e(2?j學號后兩位）n，并繪出其圖形。N=3;

n=0:0.2:N-1;

w=1;%學號01 x=exp((2+j*w)*n);subplot(2,1,1)

stem(n,real(x)),title('實部');grid on;subplot(2,1,2)

stem(n,imag(x)),title('虛部');grid on;

5、編制程序產生指數序列x?n??an，并繪出其圖形。其中a=學號后兩位、a=1/“學號后兩位”。

(1)N=10;

n=0:N-1;

a=1;%學號01 x=a.^n;stem(n,x);title('指數序列');grid on;

(2)N=10;

n=0:N-1;

a=1;%學號01 x=a.^(-n);stem(n,x);title('指數序列');grid on;

實驗三 離散時間信號的頻域分析

一、實驗目的

信號的頻域分析是信號處理中一種有效的工具。在離散信號的頻域分析中，通常將信號表示成單位采樣序列的線性組合，而在頻域中，將信號表示成復變量或的線性組合。通過這樣的表示，可以將時域的離散序列映射到頻域以便于進一步的處理。

在本實驗中，將學習利用MATLAB計算離散時間信號的DTFT和DFT，并加深對其相互關系的理解。

二、實驗原理

（１）DTFT和DFT的定義及其相互關系。

（２）使用到的MATLAB命令有基于DTFT離散時間信號分析函數以及求解序列的DFT函數。

三、實驗內容及分析

（１）編程計算并畫出下面DTFT的實部、虛部、幅度和相位譜。

X(e)?jw0.0518?0.1553e1?1.2828ex(n)?cos?jw?jw?0.1553e?j2w?1.0388e?j2w?0.0518e?j3w?0.3418e?j3w

（２）計算32點序列

5?n16，0≦n≦31的32點和64點DFT，分別繪出幅度譜圖形，并繪出該序列的DTFT圖形。

３－１

clear;

x=[0.0518,-0.1553,0.1553,0.0518];y=[1,1.2828,1.0388,0.3418];w=[0:500]*pi/500 H=freqz(x,y,w);

magX=abs(H);angX=angle(H);realX=real(H);imagX=imag(H);subplot(221);plot(w/pi,magX);grid;

xlabel('frequency in pi unit');ylabel('magnitude');title('幅度 part');axis([0 0.9 0 1.1]);

subplot(223);plot(w/pi,angX);grid;

xlabel('frequency in pi unit');ylabel('radians');title('相位 part');axis([0 1-3.2 3.2]);

subplot(222);plot(w/pi,realX);grid;

xlabel('frequency in pi unit');ylabel('real part');title('實部 part');axis([0 1-1 1]);

subplot(224);plot(w/pi,imagX);grid;

xlabel('frequency in pi unit');ylabel('imaginary');title('虛部 part');axis([0 1-1 1.1]);

３－２

N=32;n=0:N-1;

xn=cos(5*pi*n/16);k=0:1:N-1;Xk=fft(xn,N);subplot(2,1,1);stem(n,xn);subplot(2,1,2);stem(k,abs(Xk));title('32點');figure N=64;n=0:N-1;

xn=cos(5*pi*n/16);k=0:1:N-1;Xk=fft(xn,N);subplot(2,1,1);stem(n,xn);subplot(2,1,2);stem(k,abs(Xk));title('64點');

（１）

（２）

實驗四 離散時間LTI系統的Z域分析

一、實驗目的

本實驗通過使用MATLAB函數對離散時間系統的一些特性進行仿真分析，以加深對離散時間系統的零極點、穩定性，頻率響應等概念的理解。學會運用MATLAB分析離散時間系統的系統函數的零極點；學會運用MATLAB分析系統函數的零極點分布與其時域特性的關系；學會運用MATLAB進行離散時間系統的頻率特性分析。

二、實驗原理

離散時間系統的系統函數定義為系統零狀態響應的Z變化與激勵的Z變化之比。

在MATLAB中系統函數的零極點可通過函數roots得到，也可借助函數tf2zp得到，tf2zp的語句格式為

[Z,P,K]=tf2zp(B,A)其中，B與A分別表示H（z）的分子與分母多項式的系數向量。它的作用是將H（z）的有理分式表示式轉換為零極點增益形式。

若要獲得系統函數H（z）的零極點分布圖，可直接應用zplane函數，其語句格式為

Zplane（B，A）

其中，B與A分別表示H（z）的分子和分母多項式的系數向量。它的作用是在z平面上畫出單位圓、零點與極點。

離散系統中z變化建立了時域函數h（n）與z域函數H（z）之間的對應關系。因此，z變化的函數H（z）從形式可以反映h（n）的部分內在性質?？筛鶕到y的傳遞函數H（z）求單位沖激響應h（n）的函數impz、filter等。

利用系統的頻率響應，可以分析系統對各種頻率成分的響應特性，并推出系統的特性（高通、低通、帶通、帶阻等）。

MATLAB提供了求離散時間系統頻響特性的函數freqz，調用freqz的格式主要有兩種。一種形式為

[H,w]= reqz(B,A,N)其中，B與A分別表示H（z）分子和分母多項式的系數向量；N為正整數，默認值為512；返回值w包含[0，π]范圍內的N個頻率等分點；返回值H則是離散時間系統頻率響應在0~π范圍內N個頻率處的值。另一種形式為

[H，w]= freqz（B，A，N，‘whole’）

與第一種方式不同之處在于角頻率的范圍由[0，π]擴展到[0，2π]。

三、實驗內容與結果分析

已知LTI離散時間系統，要求由鍵盤實現系統參數輸入，并繪出幅頻和相頻響應曲線和零極點分布圖，進而分析系統的濾波特性和穩定性。

（一）程序

b=[0.0528,0.797,0.1295,0.1295,0.797,0.0528];

a=[1,-1.8107,2.4947,-1.8801,0.9537,-0.2336];w=[0:20:500]*pi/500;

x1=0.0528+0.797*exp(-1*j*w)+0.1295*exp(-2*j*w)+0.1295*exp(-3*j*w)+0.797*exp(-4*j*w)+0.0528*exp(-5*j*w);

x2=1-1.8107*exp(-1*j*w)+2.4947*exp(-2*j*w)+1.8801*exp(-3*j*w)+0.9537*exp(-4*j*w)+0.2336*exp(-5*j*w);x22=x2+(x2==0)*eps;x=x1./x22;magx=abs(x);

angx=angle(x).*180/pi;

subplot(2,2,3);zplane(b,a);title('零極點圖');subplot(2,2,2);stem(w/pi,magx);title('幅度部分');ylabel('振幅');subplot(2,2,4);stem(w/pi,angx);

xlabel('以pi為單位的頻率');title('相位部分');ylabel('相位');

（二）波形圖

圖4-1 幅頻、相頻響應曲線、零極點分布圖

實驗六 IIR數字濾波器的設計

一、實驗目的

從理論上講，任何的線性是不變（LTI）離散時間系統都可以看做一個數字濾波器，因此設計數字濾波器實際就是設計離散時間系統。數字濾波器你包括IIR（無限沖激響應）和FIR（有限沖激響應）型，在設計時通常采用不同的方法。

本實驗通過使用MATLAB函數對數字濾波器進行設計和和實現，要求掌握IIR數字巴特沃斯濾波器、數字切比雪夫濾波器的設計原理、設計方法和設計步驟；能根據給定的濾波器指標進行濾波器設計；同時也加深學生對數字濾波器的常用指標和設計過程的理解。

二、實驗原理

在IIR濾波器的設計中，常用的方法是：先根據設計要求尋找一個合適的模擬原型濾波器，然后根據一定的準則將此模擬原型濾波器轉換為數字濾波器。

IIR濾波器的階數就等于所選的模擬原型濾波器的階數，所以其階數確定主要是在模擬原型濾波器中進行的。

IIR數字濾波器的設計方法如下：（1）沖激響應不變法。（2）雙線性變化法。

一般來說，在要求時域沖激響應能模仿模擬濾波器的場合，一般使用沖激響應不變法。沖激響應不變法一個重要特點是頻率坐標的變化是線性的，因此如果模擬濾波器的頻率響應帶限于折疊頻率的話，則通過變換后濾波器的頻率響應可不失真地反映原響應與頻率的關系。

與沖激響應不變法比較，雙線性變化的主要優點是靠頻率的非線性關系得到s平面與z平面的單值一一對應關系，整個值對應于單位圓一周。所以從模擬傳遞函數可直接通過代數置換得到數字濾波器的傳遞函數。

MATLAB提供了一組標準的數字濾波器設計函數，大大簡化了濾波器的設計工程。

（1）butter。

（2）cheby1、cheby2。

三、實驗內容及分析

利用MATLAB編程方法或利用MATLAB中fdatool工具設計不同功能的IIR數字濾波器。

1、基于chebyshev I型模擬濾波器原型使用沖激不變轉換方法設計數字濾波器，要求參數為通帶截止頻率?p?0.4?；通帶最大衰減Ap?1dB；阻帶截止頻率?s?0.4?；阻帶最小衰減As?35dB。

程序:

wp=0.2*pi;

%通帶邊界頻率

ws=0.4*pi;

%阻帶截止頻率 rp=1;

%通帶最大衰減 rs=35;

%阻帶最小衰減

Fs=1000;

%?ùéè3é?ù??3?1000hz

[N,Wn]=cheb1ord(wp,ws,rp,rs,'s');

[Z,P,K]=cheby1(N,rp,Wn,'s');[H,W]=zp2tf(Z,P,K);

figure(1);freqs(H,W);[P,Q]=freqs(H,W);figure(2);plot(Q*Fs/(2*pi),abs(P));grid on;

xlabel('頻率/Hz');ylabel('幅度');

2、基于Butterworth型模擬濾波器原型使用雙線性變換方法設計數字濾波器的，要求參數為截止頻率?p?0.4?；通帶最大衰減Ap?1dB；阻帶截止頻率?s?0.25?；阻帶最小衰減AS?40dB。程序: wp=0.4*pi;ws=0.25*pi;rp=1;rs=40;fs=500;ts=1/fs;wp1=wp*ts;ws1=ws*ts;

wp2=2*fs*tan(wp1/2);ws2=2*fs*tan(ws1/2);

[N,Wn]=buttord(wp2,ws2,rp,rs,'s');[Z,P,K]=buttap(N);[Bap,Aap]=zp2tf(Z,P,K);[b,a]=lp2lp(Bap,Aap,Wn);[bz,az]=bilinear(b,a,fs);[H,W]=freqz(bz,az);subplot(2,1,1);plot(W/pi,abs(H));grid on;xlabel('頻率')ylabel('幅度')subplot(2,1,2);

plot(W/pi,20*log10(abs(H)));grid on;xlabel('頻率');ylabel('幅度(dB)');

實驗七 FIR數字濾波器的設計

一、實驗目的

掌握用窗函數設計FIR數字濾波的原理及其設計步驟；熟悉線性相位數字濾波器的特性。學習編寫數字濾波器的設計程序的方法，并能進行正確編程；根據給定的濾波器指標，給出設計步驟。

二、實驗原理

如果系統的沖激響應h（n）為已知，則系統的輸入輸出關系為

y（n）=x（n）*h（n）

對于低通濾波器，只要設計出低通濾波器的沖激響應函數，就可以由式得到系統的輸出了。

但是將h（n）作為濾波器的脈沖響應有兩個問題：一是它是無限長的；二是它是非因果的。對此，采取兩項措施：一是將h（n）截短；二是將其右移。

設計時，要根據阻帶的最小衰減和過渡帶寬度來選擇恰當的窗函數類型和窗口長度N。常用的窗函數有矩形窗、海明窗和布萊克曼窗等。

窗函數設計FIR濾波器步驟如下：

（1）給定理想頻率響應的幅頻特性和相頻特性；

（2）求理想單位脈沖響應，在實際計算中，可對理想頻率響應采樣。（3）根據過渡帶寬度和阻帶最小衰減，確定窗函數類型和窗口長度N；（4）求FIR濾波器單位脈沖響應；

（5）分析幅頻特性，若不滿足要求，可適當改變窗函數形式或長度N，重復上述設計過程，以得到滿意的結果。

三、實驗內容及分析

1、分別用海明窗和布萊克曼窗設計一個48階的FIR帶通濾波器，通帶為Wn??0.450.55?。程序1:海明窗設計

N=48;

Window=hamming(N+1);w1=0.45;w2=0.55;ws=[w1,w2];

b=fir1(N,ws/pi,Window);freqz(b,1,512);title('海明窗');grid on;

程序2:萊克曼窗設計

N=48;

Window=blackman(N+1);w1=0.45;w2=0.55;ws=[w1,w2];

b=fir1(N,ws/pi,Window);freqz(b,1,512);title('布萊克曼窗');grid on;

2、用矩形窗設計一個線性相位高通濾波器。其中He??jw?e?j????????00.3?????

0???0.3?程序: N=9;

alpha=(N-1)/2;Wc=0.7*pi;n=(0:8);i=n-alpha;i=i+(i==0)*eps;

h=(-1).^n.*sin((i).*Wc)./((i).*pi);%矩形窗函數設計的系統脈沖響應 w=(0:1:500)*2*pi/500;

H=h*exp(-j*n'*w);%矩形窗函數設計的頻響 magH=abs(H);% 矩形窗函數設計的振幅 subplot(211);stem(n,h);

axis([0,8,-0.4,0.4]);title('矩形窗設計h(n)');line([0,10],[0,0]);xlabel('n');ylabel('h');subplot(212);plot(w/pi,magH);

xlabel('以pi為單位的頻率');ylabel('H振幅');axis([0,2,0,1.7]);title('矩形窗設計振幅譜');

實驗心得體會：

這次實驗使我進一步加深了對MATLAB軟件的使用。從上次的信號系統實驗的初步使用到這一次的深入了解，有了更深刻的認識。對這種語言環境也有了新的了解。

在實驗的過程中，我對數字濾波器的整個過程有了很好的理解和掌握。IIR數字濾波器的設計讓我知道了巴特沃思濾波器和切比雪夫濾波器的頻率特性，還有雙線性變換及脈沖響應不變法設計的濾波器的頻率特性。做這兩個實驗的時候程序有點困難，但經過細心的改寫圖形最終出來了。FIR數字濾波器的設計出來的是兩種窗的圖形，通過兩種窗的比較，我了解了他們各自的特點，幅頻和相頻特性。

最后，感謝張老師對我的諄諄教導！