第一篇：MATLAB實現(xiàn)數(shù)字信號處理

數(shù)字信號處理

說 明 書

目錄

一．摘要…………………………………3 二．課程設計目的………………………3 三．設計內(nèi)容……………………………3 四．設計原理……………………………4 4.1．語音信號的采集…………………………….4 4.2．濾波器……………………………………….4 4．21．IIR濾波器原理…………………………………….4 4.22．FIR濾波器原理………………………………………5 五．設計步驟……………………………6 5.1錄制女音………………………………………6 5.2采樣語音信號并畫出時域波形和頻譜圖……7 5.3采用雙線性變換法設計IIR濾波器…………10 5.4窗函數(shù)法設計FFR濾波器………………......12 5.5用IIR濾波器對信號進行濾波………………14 5.6用FIR濾波器對信號進行濾波………………16 5.7男女聲語音信號頻譜特點分析………………19 5.8有背景噪聲的信號分析………………………20 六．心得體會…………………………….22 七．參考文獻…………………………….23

一．摘要：

這次課程設計的主要目的是綜合運用本課程的理論知識進行頻譜分析以及濾波器設計，通過理論推導得出相應結(jié)論，并利用MATLAB或者ＤＳＰ開發(fā)系統(tǒng)作為工具進行實現(xiàn)，從而復習鞏固課堂所學的理論知識，提高對所學知識的綜合應用能力，并從實踐上初步實現(xiàn)對數(shù)字信號的處理。通過對聲音的采樣，將聲音采樣后的頻譜與濾波。

MATLAB全稱是Matrix Laboratory，是一種功能強大、效率高、交互性好的數(shù)值和可視化計算機高級語言，它將數(shù)值分析、矩陣運算、信號處理和圖形顯示有機地融合為一體，形成了一個極其方便、用戶界面友好的操作環(huán)境。經(jīng)過多年的發(fā)展，已經(jīng)發(fā)展成為一種功能全面的軟件，幾乎可以解決科學計算中所有問題。MATLAB軟件還提供了非常廣泛和靈活的用于處理數(shù)據(jù)集的數(shù)組運算功能。

在本次課程設計中，主要通過MATLAB來編程對語音信號處理與濾波，設計濾波器來處理數(shù)字信號并對其進行分析。

二．課程設計目的：

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

三．設計內(nèi)容：

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

四．設計原理：

4.1．語音信號的采集

熟悉并掌握MATLAB中有關(guān)聲音（wave）錄制、播放、存儲和讀取的函數(shù)，在MATLAB環(huán)境中，有關(guān)聲音的函數(shù)有：

a：y=wavrecord(N,fs,Dtype);利用系統(tǒng)音頻輸入設備錄音，以fs為采樣頻率，默認值為11025，即以11025HZ進行采樣。Dtype為采樣數(shù)據(jù)的存儲格式，用字符串指定，可以是：‘double’、‘single’、’int16’、‘int8’其中只有int8是采用8位精度進行采樣，其它三種都是16位采樣結(jié)果轉(zhuǎn)換為指定的MATLAB數(shù)據(jù)；

b：wavplay(y,fs)；利用系統(tǒng)音頻輸出設備播放，以fs為播放頻率，播放語音信號y；

c：wavwrite((y,fs,wavfile);創(chuàng)建音頻文件； d：y=wavread(file);讀取音頻文件；

關(guān)于聲音的函數(shù)還有sound()；soundsc()；等。4.2濾波器： 4．21.IIR濾波器原理

沖激響應不變法是使數(shù)字濾波器在時域上模擬濾波器，但是它們的缺點是產(chǎn)生頻率響應的混疊失真，這是由于從s平面到z平面是多值的映射關(guān)系所造成的。

雙線性變換法是使數(shù)字濾波器的頻率響應與模擬濾波器的頻率響應相似的一種變換方法。為了克服多值映射這一缺點，我們首先把整個s平面壓縮變換到某一中介的s1平面的一條橫帶里，再通過變換關(guān)系將此橫帶變換到整個z平面上去，這樣就使得s平面與z平面是一一對應的關(guān)系，消除了多值變換性，也 就消除了頻譜混疊現(xiàn)象。

雙線性法設計IIR數(shù)字濾波器的步驟：

1）將數(shù)字濾波器的頻率指標{ ?k}由Wk＝（2/T）*tan(wk),轉(zhuǎn)換為模擬濾波器的頻率指標{?k}.2)由模擬濾波器的指標設計H(s).3)由H(s)轉(zhuǎn)換為H(z)21?z?1H(z)?H(s)s?T1?z?1

4.22.FIR濾波器原理

FIR濾波器與IIR濾波器特點不同，IIR濾波器的相位是非線性的，若需線性相位則要采用全通網(wǎng)絡進行相位校正。而有限長單位沖激響應（FIR）數(shù)字濾波器就可以做成具有嚴格的線性相位，同時又可以具有任意的幅度特性。

由于FIR系統(tǒng)的沖激響應就是其系統(tǒng)函數(shù)各次項的系數(shù)，所以設計FIR濾波器的方法之一可以從時域出發(fā)，截取有限長的一段沖激響應作為H（z）的系數(shù)，沖激響應長度N就是系統(tǒng)函數(shù)H（z）的階數(shù)。只要N足夠長，截取的方法合理，總能滿足頻域的要求。這種時域設計、頻域檢驗的方法一般要反復幾個回合，不像IIR DF設計靠解析公式一次計算成功。給出的理想濾波器頻率響應是，它是w的周期函數(shù)，周期

由傅立葉反變換導出，即

hd(n)?1Hd(ejw)ejwndw?2?，再將hd(n)與窗函數(shù)，因此可展開成傅氏級數(shù)w(n)相乘就可以得到h(n)。、的計算可采用傅氏變換的現(xiàn)成公式和程序，窗函數(shù)也是現(xiàn)成的。但整個設計過程不能一次完成，因為窗口類型和大小的選擇沒有解析公式可一次算，整個設計可用計算機編程來做。

窗函數(shù)的傅式變換W(ejω)的主瓣決定了H(ejω)過渡帶寬。W(ejω)的旁瓣大小和多少決定了H(ejω)在通帶和阻帶范圍內(nèi)波動幅度，常用的幾種窗函數(shù)有：

矩形窗

w(n)=RN(n)；

Hanning窗

；

Hamming窗

；

Blackmen窗

；

Kaiser窗。

式中Io(x)為零階貝塞爾函數(shù)。

五．設計步驟：

5.1錄制女音：

利用MATLAB中的函數(shù)錄制聲音。function nvyin()fs=11025;

%采樣頻率

str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

開始錄音');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);y=wavrecord(3*fs,fs,'double');

%錄制聲音3秒

str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

錄音結(jié)束');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

播放錄音');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];5 disp(str);wavplay(y,fs);

%播放錄音

str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

播放錄音結(jié)束');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);wavwrite(y,fs,'原女音');

%聲音的存儲

5.2采樣語音信號并畫出時域波形和頻譜圖

讀取語音信號，畫出其時域波形和頻譜圖，與截取后的語音信號的時域波形和頻譜圖比較，觀察其變化。程序如下：

[x,fs,bits]=wavread('女音.wav');

%讀取聲音

N=length(x);

%計數(shù)讀取信號的點數(shù) t=(1:N)/fs;

%信號的時域采樣點 f0=fs/N;

%采樣間隔 n=1:N/2;

%取信號的一半 figure(1);subplot(2,2,1);

%把畫圖區(qū)域劃分為2行2列，指定第一個圖 plot(t, x);

%畫出聲音采樣后的時域波形 title('原女音信號的時域波形');

%給圖形加注標簽說明 xlabel('時間/t');ylabel('振幅/A');grid;

%添加網(wǎng)格

y=fft(x);

%對信號做N點FFT變換 k=(n-1)*f0;

%頻域采樣點

subplot(2,2,3);

%把畫圖區(qū)域劃分為2行2列，指定第三個圖 plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應圖 title('FFT變換后聲音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;

%添加網(wǎng)格

subplot(2,2,4);

%把畫圖區(qū)域劃分為2行2列，指定第四個圖 if y~=0

%判斷指數(shù)是否為0

plot(k,20*log10(abs(y(n))));

%畫信號頻譜的分貝圖 end xlabel('Hz');ylabel('振幅/分貝');title('FFT變換后聲音的頻譜特性');grid;

%添加網(wǎng)格

%實際發(fā)出聲音落后錄制動作半拍的現(xiàn)象的解決 siz=wavread('女音.wav','size');x1=wavread('女音.wav',[3500 32076]);

%截取語音信號 N=length(x1);

%計數(shù)讀取信號的點數(shù) t=(1:N)/fs;

%信號的時域采樣點 f0=fs/N;

%采樣間隔 n=1:N/2;

%取信號的一半

figure(2);subplot(2,2,1);

%把畫圖區(qū)域劃分為2行2列，指定第一個圖 plot(t,x1);

%畫出聲音采樣后的時域波形 title('截取后女音信號的時域波形');

%給圖形加注標簽說明 xlabel('時間/t');ylabel('振幅/A');grid;

%添加網(wǎng)格

y1=fft(x1);

%對信號做N點FFT變換

subplot(2,2,3);

%把畫圖區(qū)域劃分為2行2列，指定第三個圖 k=(n-1)*f0;

%頻域采樣點

plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應圖 title('FFT變換后聲音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;

%添加網(wǎng)格

subplot(2,2,4);

%把畫圖區(qū)域劃分為1行2列，指定第二個圖 if y1~=0

%判斷指數(shù)是否為0

plot(k,20*log10(abs(y1(n))));

%畫信號頻譜的分貝圖 end xlabel('Hz');ylabel('振幅/分貝');title('FFT變換后聲音的頻譜特性');grid;

%添加網(wǎng)格

原女音信號的時域波形10.5A/幅0振-0.5-10123時間/tFFT變換后聲音的頻譜特性FFT變換后聲音的頻譜特性30050A200貝/值分/幅0幅100振00200040006000-***頻率/HzHz 截取后女音信號的時域波形10.5振幅/A0-0.5-10123FFT變換后聲音的頻譜特性50時間/tFFT變換后聲音的頻譜特性300200振幅/分貝幅值/A01000020004000頻率/Hz6000-5002000Hz40006000

結(jié)果分析：

由原女音信號的時域波形可知錄取開始時實際發(fā)出聲音大概落后3500個采樣點，我們把前3500點去除即可解決實際發(fā)出聲音落后錄制動作半拍的現(xiàn)象。由原女音的的頻譜圖和截取后聲音的頻譜圖可看出，對聲音的截取并不會影響它們頻譜分布。

5.3采用雙線性變換法設計IIR濾波器：

人的聲音頻率一般在（1~~4）kHZ之間，則我們只需要設計一個帶通濾波器即可濾去聲音頻帶以外的無用噪聲，得到比較清晰的聲音。根據(jù)聲音的頻譜圖分析，設計一個帶通濾波器性能指標如下：

fp1＝1000 Hz，fp2＝3000 Hz，fsc1＝500 Hz，fsc2＝3500Hz，As＝100dB，Ap＝1dB，fs=10000 程序如下：

%iir帶通的代碼： %w=2*pi*f/fs Ap=1;

%通帶波紋系數(shù)

Az=100;

%最小阻帶衰減

wp=[0.2 0.6];

%歸一化通帶數(shù)字截止頻率 wz=[0.1 0.7];

%歸一化阻帶數(shù)字截止頻率 [N,wn]=cheb1ord(wp,wz,Ap,Az);

%估計契比雪夫I型濾波器階數(shù) [b,a]=cheby1(N,Ap,wn);

%N指定濾波器階數(shù)，wn歸一化

截 %止頻率，Ap通帶波動

[h,w]=freqz(b,a);

%求數(shù)字濾波器的復頻率響應 figure(1);subplot(2,1,1);plot(w/pi,abs(h));

%繪制數(shù)字濾波器的頻譜圖 grid;xlabel('omega/pi');ylabel('振幅（幅值）');title('契比雪夫Ⅰ型帶通濾波器的幅頻響應');subplot(2,1,2);if abs(h)~=0

%判斷指數(shù)是否為0

plot(w/pi,20*log10(abs(h)));

%繪制數(shù)字濾波器頻譜的分貝圖 end grid;xlabel('omega/pi');ylabel('振幅(分貝)');title('契比雪夫Ⅰ型帶通濾波器的幅頻響應');契比雪夫Ⅰ型帶通濾波器的幅頻響應1振幅（幅值）0.5000.10.20.50.60.70.8?/?契比雪夫Ⅰ型帶通濾波器的幅頻響應0.30.40.910振幅(分貝)-200-400-60000.10.20.30.40.5?/?0.60.70.80.91

5.4窗函數(shù)法設計FFR濾波器

線性相位FIR濾波器通常采用窗函數(shù)法設計。窗函數(shù)法設計FIR濾波器的基本思想是：根據(jù)給定的濾波器技術(shù)指標，選擇濾波器長度N和窗函數(shù)ω（n），使其具有最窄寬度的主瓣和最小的旁瓣。其核心是從給定的頻率特性，通過加窗確定有限長單位脈沖響應序列h(n)。工程中常用的窗函數(shù)共有6種，即矩形窗、巴特利特（Bartlett）窗、漢寧（Hanning）窗、漢明（Hamming）窗、布萊克曼（Blackman）窗和凱澤（Kaiser）窗。

這次設計我采用的是布萊克曼來設計給定數(shù)字帶通濾波器的參數(shù)如下： wp1=0.3pi, wp2=0.6pi, wz1=0.2pi, wz2=0.7pi, Ap=1dB, Az=70dB 程序如下：

Ap=1;

%通帶波紋系數(shù) Az=100;

%最小阻帶衰減 fs=10000;

%采樣頻率 wp1=0.3*pi;wp2=0.6*pi;wz1=0.2*pi;wz2=0.7*pi;wc1=(wz1+wp1)/2;wc2=(wz2+wp2)/2;deltaW=min((wp1-wz1),(wz2-wp2));

%---取兩個過渡帶中的小者 N0=ceil(2*5.5*pi/deltaW);

%---查表7-3（P342）布拉克曼窗 N=N0+mod(N0+1,2);

%---確保N為奇數(shù) hdWindow=ideallp(wc2,N)-ideallp(wc1,N);%理想帶通濾波器 wdWindow=blackman(N);

%布拉克曼窗 hr=wdWindow.*hdWindow';

%點乘

n=0:N-1;

%階數(shù) subplot(2,2,1);stem(n,wdWindow);

%繪制布拉克曼窗時域波形 xlabel('時間');ylabel('振幅');title('布拉克曼窗');[H,W]=freqz(hr,1);

%求濾波器頻率響應 subplot(2,2,3);plot(W/pi,abs(H))

%繪制濾波器頻域波形 xlabel('omega/pi');ylabel('振幅');title('FIR帶通濾波器幅頻特性');subplot(2,2,4);

if abs(H)~=0

%判斷指數(shù)是否為0

plot(W/pi,20*log10(abs(H)));

%畫濾波器頻譜的分貝圖 end xlabel('omega/');ylabel('振幅/分貝');title('FIR帶通濾波器幅頻特性');grid;

%添加網(wǎng)格 %---ideallp（）函數(shù)（非系統(tǒng)自有函數(shù)）在系統(tǒng)安裝目錄的WORK子目錄ideallp.m function hd = ideallp(wc,N);% 理想低通濾波器的脈沖響應子程序 % hd = 點0 到 N-1之間的理想脈沖響應 % wc = 截止頻率(弧度)% N = 理想濾波器的長度

tao =(N-1)/2;

% 理想脈沖響應的對稱中心位置 n = [0:(N-1)];

% 設定脈沖響應長度 m = n-tao + eps;

% 加一個小數(shù)以避免零作除數(shù)

hd = sin(wc*m)./(pi*m);

% 理想脈沖響應

布拉克曼窗1振幅0.500406080時間FIR帶通濾波器幅頻特性500振幅/分貝20FIR帶通濾波器幅頻特性1.51振幅-50-100-15000.5?/10.5000.5?/?1

5.5用IIR濾波器對信號進行濾波

用自己設計的IIR濾波器分別對采集的信號進行濾波，在Matlab中，IIR濾波器利用函數(shù)filter對信號進行濾波。程序如下： [x,fs,bits]=wavread('女音.wav');N=length(x);

%計數(shù)讀取信號的點數(shù) t=(1:N)/fs;

%信號的時域采樣點 f0=fs/N;

%采樣間隔 n=1:N/2;

%取信號的一半 y=fft(x);

%對信號做N點FFT變換 k=(n-1)*f0;

%頻域采樣點

subplot(2,1,1);

%把畫圖區(qū)域劃分為2行1列，指定第一個圖 plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應圖 title('濾波前女音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;

%iir帶通的代碼：

Ap=1;

%通帶波紋系數(shù)

Az=100;

%最小阻帶衰減

wp=[0.2 0.6];

%歸一化通帶數(shù)字截止頻率 wz=[0.1 0.7];

%歸一化阻帶數(shù)字截止頻率 [N,wn]=cheb1ord(wp,wz,Ap,Az);

%估計契比雪夫I型濾波器階數(shù)

[b,a]=cheby1(N,Ap,wn);

%N指定濾波器階數(shù)，wn歸一化截止頻率，Ap通帶波動 x1=filter(b,a,x);

%對聲音濾波 wavplay(x1)wavwrite(x1,'IIR濾波后女音.wav');N=length(x1);

%計數(shù)讀取信號的點數(shù) t=(1:N)/fs;

%信號的時域采樣點 f0=fs/N;

%采樣間隔 n=1:N/2;

%取信號的一半

y=fft(x1);

%對信號做N點FFT變換 k=(n-1)*f0;

%頻域采樣點

subplot(2,1,2);

%把畫圖區(qū)域劃分為2行1列，指定第一個圖 plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應圖 title('l濾波后女音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;

濾波前女音的頻譜特性300幅值/A***030004000頻率/Hz濾波后女音的頻譜特性500060006040幅值/A***0頻率/Hz400050006000

結(jié)果分析：

由上面濾波前后的頻譜圖可看出，濾波器濾除了小于1000Hz和大于3400Hz的頻譜成分。回放語音信號，由于低頻和高頻成分被濾除，聲音變得較低沉。

5.6用FIR濾波器對信號進行濾波

用自己設計的FIR濾波器分別對采集的信號進行濾波，在Matlab中,FIR濾波器利用函數(shù)fftfilt對信號進行濾波 程序如下：

[x,fs,bits]=wavread('女音.wav');N=length(x);

%計數(shù)讀取信號的點數(shù)

t=(1:N)/fs;

%信號的時域采樣點 f0=fs/N;

%采樣間隔 n=1:N/2;

%取信號的一半

y=fft(x);

%對信號做N點FFT變換 k=(n-1)*f0;

%頻域采樣點

subplot(2,1,1);

%把畫圖區(qū)域劃分為2行1列，指定第一個圖 plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應圖 title('濾波前女音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/omega');ylabel('幅值/A');grid;

%FIR帶通濾波器代碼 fs=10000;wp1=0.3*pi;wp2=0.6*pi;wz1=0.2*pi;wz2=0.7*pi;wc1=(wz1+wp1)/2;wc2=(wz2+wp2)/2;deltaW=min((wp1-wz1),(wz2-wp2));

%---取兩個過渡帶中的小者 N0=ceil(2*5.5*pi/deltaW);

%---查表7-3（P342）布拉克曼窗 N=N0+mod(N0+1,2);

%---確保N為奇數(shù) hdWindow=ideallp(wc2,N)-ideallp(wc1,N);wdWindow=blackman(N);hr=wdWindow.*hdWindow';x1=fftfilt(hr,x);

%對聲音濾波 wavplay(x1)wavwrite(x1,'FIR濾波后女音.wav');N=length(x1);

%計數(shù)讀取信號的點數(shù) t=(1:N)/fs;

%信號的時域采樣點 f0=fs/N;

%采樣間隔 n=1:N/2;

%取信號的一半

y=fft(x1);

%對信號做N點FFT變換 k=(n-1)*f0;

%頻域采樣點

subplot(2,1,2);

%把畫圖區(qū)域劃分為2行1列，指定第一個圖 plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應圖 title('l濾波后女音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;

濾波前女音的頻譜特性300200幅值/A***004000頻率/?l濾波后女音的頻譜特性500060006040幅值/A20005001000***03000頻率/Hz***0

結(jié)果分析：

由上面濾波前后的頻譜圖可看出，濾波器濾除了小于1000Hz和大于3500Hz的頻譜成分。和用IIR濾波器濾波一樣，回放語音信號，由于低頻和高頻成分被濾除，聲音變得較低沉。5.7男女聲語音信號頻譜特點分析

換一個男音錄制同樣一段語音內(nèi)容，分析兩段內(nèi)容相同的語音信號頻譜之間有什么特點。程序如下：

[x,fs,bits]=wavread('女音.wav');

%讀取聲音

N=length(x);

%計數(shù)讀取信號的點數(shù) t=(1:N)/fs;

f0=fs/N;

n=1:N/2;

y=fft(x);

k=(n-1)*f0;

subplot(2,1,1);

plot(k,abs(y(n)));

title('FFT變換后女音的頻譜特性');xlabel('頻率/omega');ylabel('幅值/A');grid;

[x,fs,bits]=wavread('明明.wav');

N=length(x);

t=(1:N)/fs;

f0=fs/N;

n=1:N/2;

y=fft(x);

k=(n-1)*f0;

subplot(2,1,2);

plot(k,abs(y(n)));

title('FFT變換后男音的頻譜特性');xlabel('頻率/omega');ylabel('幅值/A');grid;

%信號的時域采樣點

%采樣間隔

%取信號的一半

%對信號做N點FFT變換

%頻域采樣點

%把畫圖區(qū)域劃分為2行1列，指定第一個圖%繪制原始語音信號的幅頻響應圖

%給圖形加注標簽說明

%添加網(wǎng)格

%讀取聲音

%計數(shù)讀取信號的點數(shù)

%信號的時域采樣點

%采樣間隔

%取信號的一半

%對信號做N點FFT變換

%頻域采樣點

%把畫圖區(qū)域劃分為2行1列，指定第二個圖%繪制原始語音信號的幅頻響應圖

%給圖形加注標簽說明

%添加網(wǎng)格

axis([0 6000 0 300])；

%改變橫縱坐標便于比較頻譜圖

FFT變換后女音的頻譜特性300200幅值/A***00頻率/?FFT變換后男音的頻譜特性***200幅值/A***00頻率/?400050006000

結(jié)果分析：

通過比較上面女音頻譜圖和男音頻譜圖可知，男音的頻譜集中在低頻部分，高頻成分底，譜線較平滑，聲音聽起來低沉。5.8有背景噪聲的信號分析

從硬盤中把一段噪聲（頻譜能量集中在某個小范圍內(nèi)）疊加到語音信號中，分析疊加前后信號頻譜的變化，設計一個合適的濾波器，能夠把該噪聲濾除； 程序如下：

z=wavread('女音.wav',[1 24000]);

%讀取聲音在1-24000之間 f=wavread('noise.wav',[1 24000]);x=z+f;wavplay(x);fs=11025;N=length(x);f0=fs/N;

%采樣間隔

n=1:N;

%取信號的一半 y=fft(x,N);%對信號做N點FFT變換

k=(n-1)*f0;

%頻域采樣點

subplot(2,1,1);

%把畫圖區(qū)域劃分為1行2列，指定第二個圖 plot(k,abs(y(n)));

%繪制原始語音信號的幅頻響應圖 title('加噪聲后聲音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;%添加網(wǎng)格

%iir帶通濾波器的代碼：

Ap=1;

%通帶波紋系數(shù)

Az=70;

%最小阻帶衰減

wp=[0.2 0.7];

%歸一化通帶數(shù)字截止頻率 wz=[0.1 0.8];

%歸一化阻帶數(shù)字截止頻率 [N,wn]=cheb1ord(wp,wz,Ap,Az);

%估計契比雪夫I型濾波器階數(shù)

[b,a]=cheby1(N,Ap,wn);

%N指定濾波器階數(shù)，wn歸一化截止頻率，Ap通帶波動 x1=filter(b,a,x);

%對聲音濾波 wavplay(x1)；

wavwrite(x1,'濾除噪音后女音.wav');N=length(x1);f0=fs/N;

%采樣間隔 n=1:N;

%取信號的一半

y1=fft(x1,N);

%對信號做fs點FFT變換

subplot(2,1,2);

%把畫圖區(qū)域劃分為1行2列，指定第二個圖 k=(n-1)*f0;

%頻域采樣點

plot(k,abs(y1(n)));

%繪制原始語音信號的幅頻響應圖 title('濾除噪聲后聲音的頻譜特性');

%給圖形加注標簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;%添加網(wǎng)格

加噪聲后聲音的頻譜特性3000幅值/A***0008000頻率/Hz濾除噪聲后聲音的頻譜特性***030幅值/A***000頻率/Hz80001000012000

結(jié)果分析

觀察加噪聲后聲音的頻譜圖可知，噪音頻率主要在4000Hz處，只要我們設計一個，濾波器濾除大概在4000Hz的頻譜即可，回放濾波后的語音信號，可證噪音基本濾除。

六．心得體會：

通過這次課程設計，讓我對MATLAB的基本應用有了更深的了解，還有數(shù)字信號處理在MATLAB中的一些函數(shù)的用法。通過理論推導得出相應結(jié)論，并利用MATLAB作為工具進行實現(xiàn)，從而復習鞏固課堂所學的理論知識，提高對所學知識的綜合應用能力，并從實踐上初步實現(xiàn)對數(shù)字信號的處理。

在這次實驗中，也遇到了很多問題，比如畫信號頻譜的分貝圖時（20*log10(abs(y))）指數(shù)為零時的處理。濾波器的設計也花了好大的功夫，剛開始不會設計參數(shù)，一頭霧水，通過同學的指導和討論，得知通過觀察信號的頻譜圖，看噪音頻率集中在那一部分，設計濾波器把其濾除即可。可反復設置參數(shù)直到濾波后語音信號的效果好為止。

七．參考文獻：

（1）《MATLAB LabVIEW SystemView》翁劍楓 葉志前 編著, 機械工業(yè)出版社；

（2）《MATLAB及在電子信息課程中的應用》陳懷琛 吳大正 高西全編著，電子工業(yè)出版社；

（3）《MATLAB在數(shù)字信號處理中的應用》（弟2版）薛年喜 編著，清華大學出版社；

（4）《MATLAB擴展編程》何強 何英

編著，清華大學出版社；（5）《MATLAB7簡明教程》吳清 曹輝林 編著，清華大學出版社；（6）MATLAB5.3精要.編程及高級應用》程衛(wèi)國 馮峰 王雪梅 劉藝 編著，機械工程出版社。

第二篇：數(shù)字信號處理課后習題Matlab作業(yè)

數(shù)字信號處理MATLAB

第1頁

習題數(shù)字信號處理MATLAB習題

M1-1 已知g1(t)?cos(6?t)，g2(t)?cos(14?t)，g3(t)?cos(26?t)，以抽樣頻率fsam?10Hz對上述三個信號進行抽樣。在同一張圖上畫出g1(t)，g2(t)和g3(t)及抽樣點，對所得結(jié)果進行討論。

解：

第2頁

從以上兩幅圖中均可看出，三個余弦函數(shù)的周期雖然不同，但它們抽樣后相應抽樣點所對應的值都相同。那么這樣還原回原先的函數(shù)就變成相同的，實際上是不一樣的。這是抽樣頻率太小的原因，我們應該增大抽樣頻率才能真實還原。如下圖：f=50Hz

第3頁

程序代碼

f=10;

t=-0.2:0.001:0.2;g1=cos(6.*pi.*t);g2=cos(14.*pi.*t);g3=cos(26.*pi.*t);k=-0.2:1/f:0.2;h1=cos(6.*pi.*k);h2=cos(14.*pi.*k);h3=cos(26.*pi.*k);% subplot(3,1,1);

% plot(k,h1,'r.',t,g1,'r');% xlabel('t');% ylabel('g1(t)');% subplot(3,1,2);

% plot(k,h2,'g.',t,g2,'g');% xlabel('t');% ylabel('g2(t)');% subplot(3,1,3);

% plot(k,h3,'b.',t,g3,'b');% xlabel('t');% ylabel('g3(t)');

plot(t,g1,'r',t,g2,'g',t,g3,'b',k,h1,'r.',k,h2,'g.',k,h3,'b.')

第4頁

xlabel('t');ylabel('g(t)');

legend('g1(t)','g2(t)','g3(t)');

M2-1 利用DFT的性質(zhì)，編寫一MATLAB程序，計算下列序列的循環(huán)卷積。

(1)g[k]={1,-3,4,2,0,-2,},h[k]={3,0,1,-1,2,1}；(2)x[k]=cos(?k/2),y[k]=3k,k=0,1,2,3,4,5。解：(1)循環(huán)卷積結(jié)果

6.0000-3.0000 17.0000-2.0000 7.0000-13.0000

程序代碼

第5頁

g=[1-3 4 2 0-2];h=[3 0 1-1 2 1];l=length(g);L=2*l-1;GE=fft(g,L);HE=fft(h,L);y1=ifft(GE.*HE);for n=1:l

if n+l<=L

y2(n)=y1(n)+y1(n+l);else

y2(n)=y1(n);

end end y2

stem(0:l-1,y2)xlabel('k')ylabel('y(k)')title('循環(huán)卷積')

(2)循環(huán)卷積結(jié)果

-71.0000-213.0000 89.0000 267.0000 73.0000 219.0000

第6頁

程序代碼

k=0:5;

x=cos(pi.*k./2);y=3.^k;l=length(x);L=2*l-1;GE=fft(x,L);HE=fft(y,L);y1=ifft(GE.*HE);for n=1:l

if n+l<=L

y2(n)=y1(n)+y1(n+l);

else

y2(n)=y1(n);

end end y2

stem(0:l-1,y2)xlabel('k')ylabel('y’(k)')title('循環(huán)卷積')

第7頁

M2-2 已知序列x[k]???cos(k?/2N),|k|?N

0,其他?(1)計算序列DTFT的表達式X(ej?),并畫出N=10時，X(ej?)的曲線。

(2)編寫一MATLAB程序，利用fft函數(shù)，計算N=10時，序列x[k]的DTFT在?m?2?m/N的抽樣值。利用hold函數(shù)，將抽樣點畫在X(ej?)的曲線上。

解：

(1)X(e)?DTFT{x[k]}?j?k????x[k]e??j?k?k??N?cos(k?/2N)eN?j?k

程序代碼

N=10;k=-N:N;

x=cos(k.*pi./(2*N));W=linspace(-pi,pi,512);

第8頁

X=zeros(1,length(W));for k=-N:N

X1=x(k+N+1).*exp(-j.*W.*k);X=X+X1;end

plot(W,abs(X))xlabel('W');ylabel('abs(X)');

(2)

程序代碼

N=10;k=-N:N;

x=cos(k.*pi./(2*N));X_21=fft(x,21);L=-10:10;

W=linspace(-pi,pi,1024);X=zeros(1,length(W));for k=-N:N

X1=x(k+N+1).*exp(-j.*W.*k);X=X+X1;end

第9頁

plot(W,abs(X));hold on;

plot(2*pi*L/21,fftshift(abs(X_21)),'o');xlabel('W');ylabel('abs(X)');

M2-3 已知一離散序列為x[k]?Acos?0k?Bcos[(?0???)k]。用長度N=64的Hamming窗對信號截短后近似計算其頻譜。試用不同的A和B的取值，確定用Hamming窗能分辨的最小的譜峰間隔??w?c的值。

解：f1=100Hz f2=120Hz時

2?中cN

f2=140Hz時

第10頁

f2=160Hz時

第11頁

由以上三幅圖可見

f2=140Hz時，各譜峰可分辨。則?f又

??w?c2?N

?40Hz

且

??w???T?2??fT?2??40?1 800所以c=3.2（近似值）

程序代碼

N=64;L=1024;

f1=100;f2=160;;fs=800;

A=1;B1=1;B2=0.5;B3=0.25;B4=0.05;T=1/fs;ws=2*pi*fs;k=0:N-1;

x1=A*cos(2*pi*f1*T*k)+B1*cos(2*pi*f2*T*k);x2=A*cos(2*pi*f1*T*k)+B2*cos(2*pi*f2*T*k);x3=A*cos(2*pi*f1*T*k)+B3*cos(2*pi*f2*T*k);x4=A*cos(2*pi*f1*T*k)+B4*cos(2*pi*f2*T*k);hf=(hamming(N))';x1=x1.*hf;x2=x2.*hf;x3=x3.*hf;x4=x4.*hf;

X1=fftshift(fft(x1,L));X2=fftshift(fft(x2,L));X3=fftshift(fft(x3,L));X4=fftshift(fft(x4,L));

W=T*(-ws/2+(0:L-1)*ws/L)/(2*pi);subplot(2,2,1);plot(W,abs(X1));title('A=1,B=1');xlabel('W');ylabel('X1');subplot(2,2,2);

第12頁

plot(W,abs(X2));title('A=1,B=0.5');xlabel('W');ylabel('X2');subplot(2,2,3);plot(W,abs(X3));title('A=1,B=0.25');xlabel('W');ylabel('X3');subplot(2,2,4);plot(W,abs(X4));title('A=1,B=0.05');xlabel('W');ylabel('X4');

M2-4 已知一離散序列為x[k]?cos?0k?0.75cos?1k,0?k?63。其中, ?0?2?/15，?1?2.3?/15。

(1)對x[k]做64點FFT, 畫出此時信號的譜。

(2)如果(1)中顯示的譜不能分辨兩個譜峰，是否可對(1)中的64點信號補0而分辨出兩個譜峰。通過編程進行證實，并解釋其原因。

解：(1)

第13頁

程序代碼

W0=2*pi/15;W1=2.3*pi/15;N=64;k=0:N-1;

x=cos(W0*k)+0.75*cos(W1*k);X=fft(x);

plot(k/N,abs(X));grid on;

title('64點FFT');

(2)

第14頁

第15頁

由以上三幅圖看出：不能對(1)中的64點信號補零而分辨出兩個譜峰，這樣的方法只能改變屏幕分辨率，但可以通過加hamming窗來實現(xiàn)對譜峰的分辨。程序代碼

W0=2*pi/15;W1=2.3*pi/15;N=64;L=1024;k=0:N-1;

x=cos(W0*k)+0.75*cos(W1*k);X=fft(x,L);

plot((0:L-1)/N,abs(X));grid on;

title('1024點FFT');

M2-5 已知一連續(xù)信號為x(t)=exp(-3t)u(t)，試利用DFT近似分析

第16頁

其頻譜。若要求頻率分辨率為1Hz，試確定抽樣頻率fsam、抽樣點數(shù)N以及持續(xù)時間Tp。

解：

本題使用矩形窗，則N?fsamfsam1??fsam，Tp??1 ?f1?f

第17頁

由以上三幅圖可以看出當fsam越來越大時，近似值越來越接近

第18頁

于實際值。即fsam越大擬合效果越好，造成的混疊也是在可以允許的范圍內(nèi)。程序代碼

fs=100;ws=2*pi*fs;Ts=1/fs;N=fs;

x=exp(-3*Ts*(0:N-1));y=fft(x,N);l=length(y);

k=linspace(-ws/2,ws/2,l);

plot(k,Ts*fftshift(abs(y)),'b:');hold on;

w=linspace(-ws/2,ws/2,1024);y1=sqrt(1./(9+w.^2));plot(w,y1,'r')

title('fs=100Hz時的頻譜')legend('近似值','實際值);

M2-6 試用DFT近似計算高斯信號g(t)?exp(?dt2)的頻譜抽樣值。

π?2通過和頻譜的理論值G(j?)?exp(?)比較，討論如何根據(jù)時域的信

d4d號來恰當?shù)剡x取截短長度和抽樣頻率使計算誤差能滿足精度要求。

解：

第19頁

第20頁

由以上三幅圖可以看出：

當時域截取長度相同時，抽樣間隔越小時誤差越小，當抽樣間隔一定時，時域截取長度越長，誤差越小。當取抽樣間隔為1S，時域截取長度為2S時，誤差較大，絕對誤差在0.5左右；當抽樣間隔為0,5S，時域截取長度為2S時，誤差比間隔為1S時小，絕對誤差不大于0.2；當抽樣間隔為0.5S時域截取長度為4S時，誤差更小，絕對誤差不大于0.04。因為時域截取長度越長，保留下來的原信號中的信息越多，抽樣間隔越小，頻譜越不容易發(fā)生混疊，所以所得頻譜與理論值相比，誤差更小。

程序代碼

Ts=0.5;N=4;N0=64;

k=(-N/2:(N/2))*Ts;

第21頁

x=exp(-pi*(k).^2);X=Ts*fftshift(fft(x,N0));

w=-pi/Ts:2*pi/N0/Ts:(pi-2*pi/N0)/Ts;XT=(pi/pi)^0.5*exp(-w.^2/4/pi);subplot(2,1,1)

plot(w/pi,abs(X),'-o',w/pi,XT);xlabel('omega/pi');ylabel('X(jomega)');

legend('試驗值','理論值');

title(['Ts=',num2str(Ts)subplot(2,1,2)plot(w/pi,abs(X)-XT)ylabel('實驗誤差')

xlabel('omega/pi');

'N=',num2str(N)]);第22頁

' '

第三篇：Matlab的數(shù)字信號處理課程實驗設計的論文

摘要:本文設計了一個基于Matlab的“數(shù)字信號處理”課程綜合性實驗。該實驗把“數(shù)字信號處理”課程中的許多離散的知識點串接了起來，包括采樣、量化、濾波器設計、濾波器實現(xiàn)、DFT/FFT和濾波器的有限字長效應等。教學實踐表明該實驗有利于鞏固學生課堂上學到的理論知識，提高學生的理論聯(lián)系實際的能力和動手解決問題的能力。

關(guān)鍵詞:數(shù)字信號處理;綜合性實驗;Matlab

0引言

“數(shù)字信號處理”課程的主要內(nèi)容包括z變換、離散傅里葉變換(DFT)、快速傅里葉變換(FFT)、數(shù)字濾波器設計和實現(xiàn)以及數(shù)字信號處理中的有限字長效應等等［1］。在學習理論知識的同時或之后，引入實驗將有助于學生更好地理解和掌握課程內(nèi)容［2－3］。筆者在教學過程中，設計了Matlab綜合性實驗。該實驗在不失趣味性的同時，能把該課程中許多分散的知識點串接起來。教學實踐表明，該實驗可以幫助學生更深入地理解本門課程，取得了較好的教學效果。

1綜合實驗內(nèi)容設計

筆者所設計的Matlab實驗如下:對下式所示的輸入信號進行濾波。x=sin(100πt)+sin(480πt)(1)具體步驟為(1)將輸入的模擬信號x進行采樣和量化，得到12位精度的數(shù)字信號;(2)設計一個低通無限沖激響應(IIＲ)濾波器，將輸入信號中的240Hz的干擾濾除，要求濾波器的輸出信號中240Hz處的噪聲功率比50Hz處的信號功率低60dB。(3)設計一個高通有限沖激響應(FIＲ)濾波器，將輸入信號中的50Hz的干擾濾除，要求濾波器的輸出信號中50Hz處的噪聲功率比240Hz處的信號功率低60dB。(4)對于上述兩個濾波器，要求:給出理想濾波器的傳輸函數(shù)及頻率響應;給出系數(shù)量化后所得的新的濾波器的傳輸函數(shù)及頻率響應;確定濾波器實現(xiàn)所采用的結(jié)構(gòu)，并給出該結(jié)構(gòu)中所用加法器和乘法器的位數(shù);將輸入的數(shù)字信號通過前一步實現(xiàn)的濾波器，畫出輸出信號的頻譜，確保濾波器性能滿足設計要求。順利完成上述Matlab實驗，需要解決以下問題:(1)采樣頻率和FFT點數(shù)的選取:根據(jù)采樣定理，采樣頻率只要不低于信號中所包含的最高頻率的兩倍，就可以從采樣后的離散時間信號中恢復出原始的模擬信號。根據(jù)式(1)，采樣頻率只要不小于480Hz即可。但是當需要使用FFT對信號進行頻譜分析時，在確定采樣頻率時，除了要滿足采樣定理外，還需要考慮其他條件。例如:在做FFT時，信號頻率應為頻率分辨率的整數(shù)倍，這樣才能準確地從頻譜中看到該頻率信號的功率，避免譜泄漏，即下式中的k應為整數(shù):k=ffs=N(2)其中f，fs和N分別為信號頻率、采樣頻率和FFT的點數(shù)。fs/N為頻率分辨率，N一般為2的冪次方。在k不為整數(shù)時，為了減小譜泄漏的影響，可以在做FFT之前對采樣所得的信號進行加窗處理［1］。(2)模數(shù)轉(zhuǎn)換器的實現(xiàn):實驗中要求對輸入信號進行量化，得到12位精度的數(shù)字信號。在將輸入信號進行量化時，涉及到如何確定模數(shù)轉(zhuǎn)換器的滿量程范圍、結(jié)構(gòu)、量化方式(舍入還是截斷)以及如何進行有符號數(shù)的量化等。(3)IIＲ濾波器類型的選擇和設計:雙線性變換是設計數(shù)字IIＲ濾波器的常用方法。它首先要將所要設計的數(shù)字濾波器的歸一化邊界角頻率進行預畸變，然后再設計出滿足性能要求的模擬濾波器。模擬濾波器有四種類型，分別為巴特沃斯濾波器，切比雪夫I型濾波器、切比雪夫II型濾波器以及橢圓濾波器。只有了解了這四種濾波器的特性，才能根據(jù)實際需求來選擇合適的濾波器類型。在選擇好濾波器類型后，將濾波器的性能指標輸入相應的Matlab函數(shù)，就可以得到濾波器的傳輸函數(shù)，完成濾波器的設計。以橢圓濾波器為例，可以依次調(diào)用函數(shù)elli-pord()，函數(shù)ellipap()和函數(shù)zp2tf()來獲得濾波器的階數(shù)、零極點、增益和s域傳輸函數(shù);也可以直接調(diào)用函數(shù)ellip()來得到濾波器的s域傳輸函數(shù)。最后再通過調(diào)用函數(shù)bilinear()得到相應數(shù)字濾波器的傳輸函數(shù)。(4)FIＲ濾波器的設計:在用窗函數(shù)法來設計FIＲ濾波器時，首先要根據(jù)濾波器的性能參數(shù)(如過渡帶寬度、阻帶衰減等)選取合適的窗函數(shù)以及確定窗函數(shù)的長度，之后將得到的窗函數(shù)與理想濾波器的單位脈沖響應序列相乘得到FIＲ濾波器的單位脈沖響應序列。以Kaiser窗為例，在Matlab中，函數(shù)kaiserord()用于預估FIＲ濾波器的階數(shù)，函數(shù)kaiser()用于產(chǎn)生相應長度的Kaiser窗函數(shù)，函數(shù)fir1()用于實現(xiàn)采用該Kaiser窗設計的FIＲ濾波器，輸出為濾波器的單位脈沖響應序列。(5)濾波器的實現(xiàn):在用硬件實現(xiàn)濾波器時，必須考慮濾波器的有限字長效應，即濾波器系數(shù)的量化、濾波器中加法器和乘法器的有限字長效應以及運算結(jié)果的有限字長等等。濾波器的實現(xiàn)結(jié)構(gòu)有直接型、級聯(lián)型和并聯(lián)型等。由于IIＲ濾波器存在量化噪聲的積累，所以在選擇結(jié)構(gòu)時，需要考慮各種結(jié)構(gòu)對有限字長效應的靈敏度。高階IIＲ濾波器通常采用級聯(lián)型或并聯(lián)型結(jié)構(gòu)來實現(xiàn)。Matlab中的函數(shù)residuez(B，A)用于計算傳輸函數(shù)B(z)/A(z)的留數(shù)、極點和直接項，從而得到有理式的部分分式展開;利用傳輸函數(shù)的部分分式展開，并通過適當?shù)暮喜ⅲ梢缘玫綖V波器的并聯(lián)型結(jié)構(gòu)。函數(shù)tf2sos()則可用于將傳輸函數(shù)轉(zhuǎn)換成二階節(jié)，得到濾波器的級聯(lián)型結(jié)構(gòu)。圖3給出了系數(shù)量化前后高通濾波器的頻率響應。為了能夠判斷所設計和實現(xiàn)的濾波器的性能是否達到設計指標，需要對濾波器的輸出序列做N點的FFT。這時需要注意兩點:一要能正確地區(qū)分輸出序列中的暫態(tài)響應部分和穩(wěn)態(tài)響應部分;二要從穩(wěn)態(tài)響應部分選取連續(xù)的N個輸出值做N點的FFT。

2教學反饋

根據(jù)學生上交的實驗報告，從他們所寫的實驗收獲和實驗心得可以看出這個實驗對他們學好這門功課所起的作用。總結(jié)如下:(1)本次實驗是FIＲ濾波器與IIＲ濾波器的設計，綜合使用了大量數(shù)字濾波器的設計方法，比如雙線性變換法，窗函數(shù)法等，加深了對課堂學習的理論知識的理解，如IIＲ和FIＲ濾波器的優(yōu)缺點、濾波器的暫態(tài)響應和穩(wěn)態(tài)響應、各種模擬濾波器的性能比較以及各種窗函數(shù)之間的差異等。(2)學生對采樣定理和FFT有了更深的認識，明白了采樣頻率、FFT點數(shù)等對頻譜分析結(jié)果的影響，并通過不斷的摸索與嘗試，總結(jié)出了使用FFT時的一些注意事項。(3)對數(shù)字信號處理中的有限字長效應有了更加直觀的體會，認識到在設計濾波器的傳輸函數(shù)時，需要考慮量化對濾波器性能的影響，設計指標需要留出一定的裕量。(4)提高了用Matlab實現(xiàn)數(shù)字信號處理功能的能力，包括:熟悉了使用Matlab設計FIＲ和IIＲ濾波器的流程;學會使用Matlab中的一些函數(shù)，如fft，cheb1ord，cheby，bilinear，fir1等;學會了用Matlab編寫程序來實現(xiàn)指定結(jié)構(gòu)的濾波器;學會了從時域和頻域觀察濾波器的輸出是否正確以及是否達到性能要求等。總而言之，通過這次實驗，使學生真正了解了如何利用Matlab來進行濾波器的設計，感覺受益匪淺，對他們學好“數(shù)字信號處理”課程很有幫助。

3結(jié)語

筆者所設計的基于Matlab的綜合性實驗涵蓋了“數(shù)字信號處理”課程中的主要知識點。從學生反饋的意見可以看出，本實驗取得了良好的教學效果，這有利于提高學生學習興趣以及增強他們解決實際問題的能力。

參考文獻:

［1］程佩青，數(shù)字信號處理教程［M］，北京:清華大學出版社，2007．

［2］曹建玲，劉煥淋，雷宏江．基于MATLAB的“數(shù)字信號處理”仿真實驗［J］．北京:中國電力教育，2012(32):88－89．

［3］易婷．“數(shù)字信號處理”課程課內(nèi)配套實驗的設計［J］．南京:電氣電子教學學報，2013，35(4):89－90．

第四篇：《數(shù)字信號處理原理及實現(xiàn)》課程小結(jié)

時間過得好快，轉(zhuǎn)眼半學期結(jié)束了。這半學期數(shù)字信號的學習讓我受益匪淺。前兩章和信號與線性系統(tǒng)相關(guān)，介紹了離散時間信號與系統(tǒng)的時域分析方法最深刻的是采樣，時域采樣定理與采樣恢復和采樣內(nèi)插公式。第二章介紹了離散時間信號與系統(tǒng)的頻域分析，DTFT與Z變換，系統(tǒng)函數(shù)的零極點分布。第三章主要講了離散傅里葉變換DFT及其性質(zhì)，和頻域采樣定理。第四章介紹了傅里葉快速變換FFT，熟悉了其原理特點及方法。五六兩章分別介紹了IIR和FIR濾波器，知道了IIR的脈沖響應不變法與雙線性變換法及其優(yōu)缺點，并學會其MATLAB應用設計濾波器，F(xiàn)IR的窗函數(shù)法與頻域采樣法設計濾波器及其MATLAB實現(xiàn)。第七章主要介紹了IIR和FIR濾波器的基本網(wǎng)絡結(jié)構(gòu)，通過老師上課習題的練習基本掌握了其結(jié)構(gòu)圖的畫法。

先說說對課程的建議吧，張曉光老師是個很負責講課思路也很清晰的老師，知道從學生的角度來講問題，根據(jù)學生的反應來調(diào)整課程進度。我們都很喜歡這樣的老師，老師開新課之前總是列提綱復習上節(jié)課講的知識，每章結(jié)束都根據(jù)章節(jié)的重要性開一節(jié)總結(jié)課，這種方式個人覺得很好，希望老師堅持。但是，感覺老師講題講的不是很多，或許是課時原因，但我覺得每章結(jié)束后開一節(jié)例題課，把知識點融進去，畢竟大學生現(xiàn)在做題比較少，這樣強制一下效果會更好。這次考試的試題覺得有不少都見過，有的是課后題，但做起來還是有點吃力，應該就是習題練的少，計算跟不上去。至于教材，我覺得編的很好，每章都有相關(guān)的MATLAB編程方法，在原理講清之后就來實踐，免去了學生盲目做實驗，提高了效率。還有就是老師也很重視實驗，總是把相關(guān)的MATLAB語句語義講解清楚，這樣我們在編程序時也就相對容易點。但好像老師講程序時都注重程序的意思了，希望老師以后再講程序時把它先部分后整體，就是在講完程序意思后把程序設計思路或框架結(jié)構(gòu)，及各部分要實現(xiàn)什么再講講，這樣有助于學生設計時設計思路更清晰。再說說考試，老師分卷面成績和實驗成績及平時成績，將實驗單獨考試，可見對實驗的重視，也說明MATLAB的重要性，這樣確實提高了學生的重視心理，雖然實驗做完了，但做完50道題并看完相關(guān)講解，我又收獲了不少，理清了設計方法與思路，所以我覺的考試方式還是挺不錯的，鍛煉了我們各方面的知識。

數(shù)字信號課程結(jié)束了，真希望您還能教我們別的課。

小組成員：陳文斌、李亞偉、王猛、汪子雄、吳官寶

第五篇：數(shù)字信號處理實驗-FFT的實現(xiàn)

實

驗

報

告

學生姓名：

學 號：

指導教師：

一、實驗室名稱：數(shù)字信號處理實驗室

二、實驗項目名稱：FFT的實現(xiàn)

三、實驗原理：

一．FFT算法思想：

1．DFT的定義：

對于有限長離散數(shù)字信號{x[n]}，0 ? n ? N-1，其離散譜{x[k]}可以由離散付氏變換（DFT）求得。DFT的定義為：

N?1X[k]?通常令e?j2?N?x[n]en?0?j2?Nnk，k=0,1,…N-1 ?WN，稱為旋轉(zhuǎn)因子。

2．直接計算DFT的問題及FFT的基本思想：

由DFT的定義可以看出，在x[n]為復數(shù)序列的情況下，完全直接運算N點DFT需要（N-1）2次復數(shù)乘法和N（N-1）次加法。因此，對于一些相當大的N值（如1024）來說，直接計算它的DFT所作的計算量是很大的。

FFT的基本思想在于，將原有的N點序列分成兩個較短的序列，這些序列的DFT可以很簡單的組合起來得到原序列的DFT。例如，若N為偶數(shù)，將原有的N

22點序列分成兩個（N/2）點序列，那么計算N點DFT將只需要約[(N/2)·2]=N/2次復數(shù)乘法。即比直接計算少作一半乘法。因子（N/2）2表示直接計算（N/2）點DFT所需要的乘法次數(shù)，而乘數(shù)2代表必須完成兩個DFT。上述處理方法可以反復使用，即（N/2）點的DFT計算也可以化成兩個（N/4）點的DFT（假定N/2為偶數(shù)），從而又少作一半的乘法。這樣一級一級的劃分下去一直到最后就劃分成兩點的FFT運算的情況。

3．基2按時間抽取（DIT）的FFT算法思想：

設序列長度為N?2L，L為整數(shù)（如果序列長度不滿足此條件，通過在后面補零讓其滿足）。

將長度為N?2L的序列x[n](n?0,1,...,N?1)，先按n的奇偶分成兩組：

x[2r]?x1[r]x[2r?1]?x2[r]，r=0,1,…,N/2-1 DFT化為：

N?1N/2?1N/2?1X[k]?DFT{x[n]}?N/2?1?n?0x[n]WnkN?2rk?r?0x[2r]W2rkN??r?0x[2r?1]WN(2r?1)kN/2?1???r?0N/2?1x1[r]Wx1[r]W2rkN?W?WkN?r?0N/2?1x2[r]WN

?r?0rkN/2kN?r?0x2[r]WN/22rkrk上式中利用了旋轉(zhuǎn)因子的可約性，即：WNN/?21NrkN?/21rk?WN/2。又令

rkX1[k]??r?0x[1r]W,/X2[k]?2?r?0x[r]WN2，則上式可以寫成： /2X[k]?X1[k]?WNX2[k](k=0,1,…,N/2-1)

k可以看出，X1[k],X2[k]分別為從X[k]中取出的N/2點偶數(shù)點和奇數(shù)點序列的N/2點DFT值，所以，一個N點序列的DFT可以用兩個N/2點序列的DFT組合而成。但是，從上式可以看出，這樣的組合僅表示出了X[k]前N/2點的DFT值，還需要繼續(xù)利用X1[k],X2[k]表示X[k]的后半段本算法推導才完整。利用旋轉(zhuǎn)因子的周期性，有：WN/2?WN/2X1[N2N/2?1rkr(k?N/2)，則后半段的DFT值表達式：

rk?k]??r?0x1[r]W2N/2r(N?k)N/2?1??r?0x1[r]WN/2?X1[k]，同樣，X2[N2?k]?X2[k]

（k=0,1,…,N/2-1），所以后半段（k=N/2,…,N-1）的DFT值可以用前半段k值表達式獲得，中間還利用到WN(N2?k)N?WN2Wk得到后半段的X[k]值表達式??W，k為：X[k]?X1[k]?WNkX2[k]（k=0,1,…,N/2-1）。

這樣，通過計算兩個N/2點序列x1[n],x2[n]的N/2點DFTX1[k],X2[k]，可以組合得到N點序列的DFT值X[k]，其組合過程如下圖所示：

X1[k] X1[k]?WNkX2[k]

X2[k] WNnk-1 X1[k]?WNkX2[k]

比如，一個N = 8點的FFT運算按照這種方法來計算FFT可以用下面的流程圖來表示：

x(0)W0x(1)W0x(2)W0x(3)W2W0W1W0x(5)W0x(6)W0x(7)W2X(7)W3X(6)W2X(5)X(3)X(2)X(1)X(0)x(4)X(4)

4．基2按頻率抽取（DIF）的FFT算法思想：

設序列長度為N?2L，L為整數(shù)（如果序列長度不滿足此條件，通過在后面補零讓其滿足）。

在把X[k]按k的奇偶分組之前，把輸入按n的順序分成前后兩半：

N?1N/2?1nkNN?1X[k]?DFT{x[n]}?N/2?1N/2?1?x[n]Wn?0?(n??n?0N2)kx[n]WnkN??n?N/2x[n]WNnk??n?0N/2?1x[n]WnkN??n?0x[n?NkN2]WNnk

?N?n?0[x[n]?x[n?N2NkN2]W2N]?WN,k?0,1,...,N?1因為W2N??1，則有WX[k]???(?1)，所以：

kkN/2?1?n?0[x[n]?(?1)x[n?N2]]?WN,k?0,1,...,N?1

nk按k的奇偶來討論，k為偶數(shù)時：

N/2?1X[2r]??n?0[x[n]?x[n?N2]]?WN,k?0,1,...,N?1 N22rnN/2?1k為奇數(shù)時：X[2r?1]?前面已經(jīng)推導過WNN/2?1?n?0[x[n]?x[n?]]?WN(2r?1)n,k?0,1,...,N?1

2rk?WN/2，所以上面的兩個等式可以寫為：

N2]]?WN/2,r?0,1,...,N/2?1 N2rnrkX[2r]??n?0[x[n]?x[n?N/2?1X[2r?1]??n?0{[x[n]?x[n?]]?WN}WN/2,r?0,1,...,N/2?1

nnr通過上面的推導，X[k]的偶數(shù)點值X[2r]和奇數(shù)點值X[2r?1]分別可以由組合而成的N/2點的序列來求得，其中偶數(shù)點值X[2r]為輸入x[n]的前半段和后半段之和序列的N/2點DFT值，奇數(shù)點值X[2r?1]為輸入x[n]的前半段和后半段之差再與WN相乘序列的N/2點DFT值。

令x1[n]?x[n]?x[n?N/2?1nN2]，x2[n]?[x[n]?x[n?N/2?1N2]]?WN，則有：

nX[2r]??n?0x1[n]?WrnN/2,X[2r?1]??n?0x2[n]?WrnN/2,r?0,1,...,N2?1

這樣，也可以用兩個N/2點DFT來組合成一個N點DFT，組合過程如下圖所示：

x[n] x[n]?x[n?N2]

x[n?N2]-1 WNn [x[n]?x[n?N2]]WNn

二．在FFT計算中使用到的MATLAB命令：

函數(shù)fft(x)可以計算R點序列的R點DFT值；而fft(x,N)則計算R點序列的N點DFT，若R>N，則直接截取R點DFT的前N點，若R

四、實驗目的：

離散傅氏變換（DFT）的目的是把信號由時域變換到頻域，從而可以在頻域分析處理信息，得到的結(jié)果再由逆DFT變換到時域。FFT是DFT的一種快速算法。在數(shù)字信號處理系統(tǒng)中，F(xiàn)FT作為一個非常重要的工具經(jīng)常使用，甚至成為DSP運算能力的一個考核因素。

本實驗通過直接計算DFT，利用FFT算法思想計算DFT，以及使用MATLAB函數(shù)中的FFT命令計算離散時間信號的頻譜，以加深對離散信號的DFT變換及FFT算法的理解。

五、實驗內(nèi)容：

a)計算實數(shù)序列x(n)?cos5?16n,0?n?256的256點DFT。

b)計算周期為1kHz的方波序列（占空比為50％，幅度取為＋/-512，采樣頻率為25kHz，取256點長度）256點DFT。

六、實驗器材（設備、元器件）：

安裝MATLAB軟件的PC機一臺，DSP實驗演示系統(tǒng)一套。

七、實驗步驟：

（1）先利用DFT定義式，編程直接計算2個要求序列的DFT值。

（2）利用MATLAB中提供的FFT函數(shù)，計算2個要求序列的DFT值。（3）(拓展要求)不改變序列的點數(shù)，僅改變DFT計算點數(shù)（如變?yōu)橛嬎?024點DFT值），觀察畫出來的頻譜與前面頻譜的差別，并解釋這種差別。通過這一步驟的分析，理解頻譜分辨力的概念，解釋如何提高頻譜分辨力。

（4）利用FFT的基本思想（基2－DIT或基2－DIF），自己編寫FFT計算函數(shù)，并用該函數(shù)計算要求序列的DFT值。并對前面3個結(jié)果進行對比。

（5）(拓展要求)嘗試對其他快速傅立葉變換算法（如Goertzel算法）進行MATLAB編程實現(xiàn)，并用它來計算要求的序列的DFT值。并與前面的結(jié)果進行對比。

（6）（拓展要求）在提供的DSP實驗板上演示要求的2種序列的FFT算法（基2－DIT），用示波器觀察實際計算出來的頻譜結(jié)果，并與理論結(jié)果對比。

八、實驗數(shù)據(jù)及結(jié)果分析：

程序：（1）對要求的2種序列直接進行DFT計算的程序

（2）對要求的2種序列進行基2－DIT和基2－DIF FFT算法程序（3）對要求的2種序列用MATLAB中提供的FFT函數(shù)進行計算的程序

結(jié)果：（1）對2種要求的序列直接進行DFT計算的頻域波形

（2）對2種要求的序列進行基2－DIT和基2－DIF FFT算法頻域波形（3）對2種要求的序列用MATLAB中提供的FFT函數(shù)計算的頻域波形。（4）（拓展要求）分析利用上面的方法畫出的信號頻譜與理論計算出來的頻譜之間的差異，并解釋這種差異。

（5）（拓展要求）保持序列點數(shù)不變，改變DFT計算點數(shù)（變?yōu)?024點），觀察頻譜的變化，并分析這種變化，由此討論如何提高頻譜分辨力的問題。

九、實驗結(jié)論：

十、總結(jié)及心得體會：

十一、對本實驗過程及方法、手段的改進建議：