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

數(shù)字信號(hào)處理

說(shuō) 明 書(shū)

目錄

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

一．摘要：

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

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

在本次課程設(shè)計(jì)中，主要通過(guò)MATLAB來(lái)編程對(duì)語(yǔ)音信號(hào)處理與濾波，設(shè)計(jì)濾波器來(lái)處理數(shù)字信號(hào)并對(duì)其進(jìn)行分析。

二．課程設(shè)計(jì)目的：

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

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

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

四．設(shè)計(jì)原理：

4.1．語(yǔ)音信號(hào)的采集

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

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

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

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

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

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

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

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

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

4.22.FIR濾波器原理

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

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

由傅立葉反變換導(dǎo)出，即

hd(n)?1Hd(ejw)ejwndw?2?，再將hd(n)與窗函數(shù)，因此可展開(kāi)成傅氏級(jí)數(shù)w(n)相乘就可以得到h(n)。、的計(jì)算可采用傅氏變換的現(xiàn)成公式和程序，窗函數(shù)也是現(xiàn)成的。但整個(gè)設(shè)計(jì)過(guò)程不能一次完成，因?yàn)榇翱陬愋秃痛笮〉倪x擇沒(méi)有解析公式可一次算，整個(gè)設(shè)計(jì)可用計(jì)算機(jī)編程來(lái)做。

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

矩形窗

w(n)=RN(n)；

Hanning窗

；

Hamming窗

；

Blackmen窗

；

Kaiser窗。

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

五．設(shè)計(jì)步驟：

5.1錄制女音：

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

%采樣頻率

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

開(kāi)始錄音');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,'原女音');

%聲音的存儲(chǔ)

5.2采樣語(yǔ)音信號(hào)并畫出時(shí)域波形和頻譜圖

讀取語(yǔ)音信號(hào)，畫出其時(shí)域波形和頻譜圖，與截取后的語(yǔ)音信號(hào)的時(shí)域波形和頻譜圖比較，觀察其變化。程序如下：

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

%讀取聲音

N=length(x);

%計(jì)數(shù)讀取信號(hào)的點(diǎn)數(shù) t=(1:N)/fs;

%信號(hào)的時(shí)域采樣點(diǎn) f0=fs/N;

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

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

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

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

%給圖形加注標(biāo)簽說(shuō)明 xlabel('時(shí)間/t');ylabel('振幅/A');grid;

%添加網(wǎng)格

y=fft(x);

%對(duì)信號(hào)做N點(diǎn)FFT變換 k=(n-1)*f0;

%頻域采樣點(diǎn)

subplot(2,2,3);

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

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

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

%添加網(wǎng)格

subplot(2,2,4);

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

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

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

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

%添加網(wǎng)格

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

%截取語(yǔ)音信號(hào) N=length(x1);

%計(jì)數(shù)讀取信號(hào)的點(diǎn)數(shù) t=(1:N)/fs;

%信號(hào)的時(shí)域采樣點(diǎn) f0=fs/N;

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

%取信號(hào)的一半

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

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

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

%給圖形加注標(biāo)簽說(shuō)明 xlabel('時(shí)間/t');ylabel('振幅/A');grid;

%添加網(wǎng)格

y1=fft(x1);

%對(duì)信號(hào)做N點(diǎn)FFT變換

subplot(2,2,3);

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

%頻域采樣點(diǎn)

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

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

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

%添加網(wǎng)格

subplot(2,2,4);

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

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

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

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

%添加網(wǎng)格

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

結(jié)果分析：

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

5.3采用雙線性變換法設(shè)計(jì)IIR濾波器：

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

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);

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

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

截 %止頻率，Ap通帶波動(dòng)

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

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

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

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

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

這次設(shè)計(jì)我采用的是布萊克曼來(lái)設(shè)計(jì)給定數(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));

%---取兩個(gè)過(guò)渡帶中的小者 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';

%點(diǎn)乘

n=0:N-1;

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

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

%求濾波器頻率響應(yīng) 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);% 理想低通濾波器的脈沖響應(yīng)子程序 % hd = 點(diǎn)0 到 N-1之間的理想脈沖響應(yīng) % wc = 截止頻率(弧度)% N = 理想濾波器的長(zhǎng)度

tao =(N-1)/2;

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

% 設(shè)定脈沖響應(yīng)長(zhǎng)度 m = n-tao + eps;

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

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

% 理想脈沖響應(yīng)

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

5.5用IIR濾波器對(duì)信號(hào)進(jìn)行濾波

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

%計(jì)數(shù)讀取信號(hào)的點(diǎn)數(shù) t=(1:N)/fs;

%信號(hào)的時(shí)域采樣點(diǎn) f0=fs/N;

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

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

%對(duì)信號(hào)做N點(diǎn)FFT變換 k=(n-1)*f0;

%頻域采樣點(diǎn)

subplot(2,1,1);

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

%繪制原始語(yǔ)音信號(hào)的幅頻響應(yīng)圖 title('濾波前女音的頻譜特性');

%給圖形加注標(biāo)簽說(shuō)明 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);

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

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

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

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

%計(jì)數(shù)讀取信號(hào)的點(diǎn)數(shù) t=(1:N)/fs;

%信號(hào)的時(shí)域采樣點(diǎn) f0=fs/N;

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

%取信號(hào)的一半

y=fft(x1);

%對(duì)信號(hào)做N點(diǎn)FFT變換 k=(n-1)*f0;

%頻域采樣點(diǎn)

subplot(2,1,2);

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

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

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

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

結(jié)果分析：

由上面濾波前后的頻譜圖可看出，濾波器濾除了小于1000Hz和大于3400Hz的頻譜成分?；胤耪Z(yǔ)音信號(hào)，由于低頻和高頻成分被濾除，聲音變得較低沉。

5.6用FIR濾波器對(duì)信號(hào)進(jìn)行濾波

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

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

%計(jì)數(shù)讀取信號(hào)的點(diǎn)數(shù)

t=(1:N)/fs;

%信號(hào)的時(shí)域采樣點(diǎn) f0=fs/N;

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

%取信號(hào)的一半

y=fft(x);

%對(duì)信號(hào)做N點(diǎn)FFT變換 k=(n-1)*f0;

%頻域采樣點(diǎn)

subplot(2,1,1);

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

%繪制原始語(yǔ)音信號(hào)的幅頻響應(yīng)圖 title('濾波前女音的頻譜特性');

%給圖形加注標(biāo)簽說(shuō)明 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));

%---取兩個(gè)過(guò)渡帶中的小者 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);

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

%計(jì)數(shù)讀取信號(hào)的點(diǎn)數(shù) t=(1:N)/fs;

%信號(hào)的時(shí)域采樣點(diǎn) f0=fs/N;

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

%取信號(hào)的一半

y=fft(x1);

%對(duì)信號(hào)做N點(diǎn)FFT變換 k=(n-1)*f0;

%頻域采樣點(diǎn)

subplot(2,1,2);

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

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

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

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

結(jié)果分析：

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

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

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

%讀取聲音

N=length(x);

%計(jì)數(shù)讀取信號(hào)的點(diǎn)數(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;

%信號(hào)的時(shí)域采樣點(diǎn)

%采樣間隔

%取信號(hào)的一半

%對(duì)信號(hào)做N點(diǎn)FFT變換

%頻域采樣點(diǎn)

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

%給圖形加注標(biāo)簽說(shuō)明

%添加網(wǎng)格

%讀取聲音

%計(jì)數(shù)讀取信號(hào)的點(diǎn)數(shù)

%信號(hào)的時(shí)域采樣點(diǎn)

%采樣間隔

%取信號(hào)的一半

%對(duì)信號(hào)做N點(diǎn)FFT變換

%頻域采樣點(diǎn)

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

%給圖形加注標(biāo)簽說(shuō)明

%添加網(wǎng)格

axis([0 6000 0 300])；

%改變橫縱坐標(biāo)便于比較頻譜圖

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

結(jié)果分析：

通過(guò)比較上面女音頻譜圖和男音頻譜圖可知，男音的頻譜集中在低頻部分，高頻成分底，譜線較平滑，聲音聽(tīng)起來(lái)低沉。5.8有背景噪聲的信號(hào)分析

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

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;

%取信號(hào)的一半 y=fft(x,N);%對(duì)信號(hào)做N點(diǎn)FFT變換

k=(n-1)*f0;

%頻域采樣點(diǎn)

subplot(2,1,1);

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

%繪制原始語(yǔ)音信號(hào)的幅頻響應(yīng)圖 title('加噪聲后聲音的頻譜特性');

%給圖形加注標(biāo)簽說(shuō)明 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);

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

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

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

%對(duì)聲音濾波 wavplay(x1)；

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

%采樣間隔 n=1:N;

%取信號(hào)的一半

y1=fft(x1,N);

%對(duì)信號(hào)做fs點(diǎn)FFT變換

subplot(2,1,2);

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

%頻域采樣點(diǎn)

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

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

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

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

結(jié)果分析

觀察加噪聲后聲音的頻譜圖可知，噪音頻率主要在4000Hz處，只要我們?cè)O(shè)計(jì)一個(gè)，濾波器濾除大概在4000Hz的頻譜即可，回放濾波后的語(yǔ)音信號(hào)，可證噪音基本濾除。

六．心得體會(huì)：

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

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

七．參考文獻(xiàn)：

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

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

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

（4）《MATLAB擴(kuò)展編程》何強(qiáng) 何英

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

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

數(shù)字信號(hào)處理MATLAB

第1頁(yè)

習(xí)題數(shù)字信號(hào)處理MATLAB習(xí)題

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

解：

第2頁(yè)

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

第3頁(yè)

程序代碼

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頁(yè)

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

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

M2-1 利用DFT的性質(zhì)，編寫一MATLAB程序，計(jì)算下列序列的循環(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頁(yè)

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頁(yè)

程序代碼

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頁(yè)

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

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

(2)編寫一MATLAB程序，利用fft函數(shù)，計(jì)算N=10時(shí)，序列x[k]的DTFT在?m?2?m/N的抽樣值。利用hold函數(shù)，將抽樣點(diǎn)畫在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頁(yè)

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頁(yè)

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]。用長(zhǎng)度N=64的Hamming窗對(duì)信號(hào)截短后近似計(jì)算其頻譜。試用不同的A和B的取值，確定用Hamming窗能分辨的最小的譜峰間隔??w?c的值。

解：f1=100Hz f2=120Hz時(shí)

2?中cN

f2=140Hz時(shí)

第10頁(yè)

f2=160Hz時(shí)

第11頁(yè)

由以上三幅圖可見(jiàn)

f2=140Hz時(shí)，各譜峰可分辨。則?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頁(yè)

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)對(duì)x[k]做64點(diǎn)FFT, 畫出此時(shí)信號(hào)的譜。

(2)如果(1)中顯示的譜不能分辨兩個(gè)譜峰，是否可對(duì)(1)中的64點(diǎn)信號(hào)補(bǔ)0而分辨出兩個(gè)譜峰。通過(guò)編程進(jìn)行證實(shí)，并解釋其原因。

解：(1)

第13頁(yè)

程序代碼

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點(diǎn)FFT');

(2)

第14頁(yè)

第15頁(yè)

由以上三幅圖看出：不能對(duì)(1)中的64點(diǎn)信號(hào)補(bǔ)零而分辨出兩個(gè)譜峰，這樣的方法只能改變屏幕分辨率，但可以通過(guò)加hamming窗來(lái)實(shí)現(xiàn)對(duì)譜峰的分辨。程序代碼

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點(diǎn)FFT');

M2-5 已知一連續(xù)信號(hào)為x(t)=exp(-3t)u(t)，試?yán)肈FT近似分析

第16頁(yè)

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

解：

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

第17頁(yè)

由以上三幅圖可以看出當(dāng)fsam越來(lái)越大時(shí)，近似值越來(lái)越接近

第18頁(yè)

于實(shí)際值。即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時(shí)的頻譜')legend('近似值','實(shí)際值);

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

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

d4d號(hào)來(lái)恰當(dāng)?shù)剡x取截短長(zhǎng)度和抽樣頻率使計(jì)算誤差能滿足精度要求。

解：

第19頁(yè)

第20頁(yè)

由以上三幅圖可以看出：

當(dāng)時(shí)域截取長(zhǎng)度相同時(shí)，抽樣間隔越小時(shí)誤差越小，當(dāng)抽樣間隔一定時(shí)，時(shí)域截取長(zhǎng)度越長(zhǎng)，誤差越小。當(dāng)取抽樣間隔為1S，時(shí)域截取長(zhǎng)度為2S時(shí)，誤差較大，絕對(duì)誤差在0.5左右；當(dāng)抽樣間隔為0,5S，時(shí)域截取長(zhǎng)度為2S時(shí)，誤差比間隔為1S時(shí)小，絕對(duì)誤差不大于0.2；當(dāng)抽樣間隔為0.5S時(shí)域截取長(zhǎng)度為4S時(shí)，誤差更小，絕對(duì)誤差不大于0.04。因?yàn)闀r(shí)域截取長(zhǎng)度越長(zhǎng)，保留下來(lái)的原信號(hào)中的信息越多，抽樣間隔越小，頻譜越不容易發(fā)生混疊，所以所得頻譜與理論值相比，誤差更小。

程序代碼

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

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

第21頁(yè)

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('試驗(yàn)值','理論值');

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

xlabel('omega/pi');

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

' '

第三篇：Matlab的數(shù)字信號(hào)處理課程實(shí)驗(yàn)設(shè)計(jì)的論文

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

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

0引言

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

1綜合實(shí)驗(yàn)內(nèi)容設(shè)計(jì)

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

2教學(xué)反饋

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

3結(jié)語(yǔ)

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

參考文獻(xiàn):

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

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

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

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

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

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

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

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

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

實(shí)

驗(yàn)

報(bào)

告

學(xué)生姓名：

學(xué) 號(hào)：

指導(dǎo)教師：

一、實(shí)驗(yàn)室名稱：數(shù)字信號(hào)處理實(shí)驗(yàn)室

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

三、實(shí)驗(yàn)原理：

一．FFT算法思想：

1．DFT的定義：

對(duì)于有限長(zhǎng)離散數(shù)字信號(hào){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．直接計(jì)算DFT的問(wèn)題及FFT的基本思想：

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

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

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

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

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

將長(zhǎng)度為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點(diǎn)偶數(shù)點(diǎn)和奇數(shù)點(diǎn)序列的N/2點(diǎn)DFT值，所以，一個(gè)N點(diǎn)序列的DFT可以用兩個(gè)N/2點(diǎn)序列的DFT組合而成。但是，從上式可以看出，這樣的組合僅表示出了X[k]前N/2點(diǎn)的DFT值，還需要繼續(xù)利用X1[k],X2[k]表示X[k]的后半段本算法推導(dǎo)才完整。利用旋轉(zhuǎn)因子的周期性，有：WN/2?WN/2X1[N2N/2?1rkr(k?N/2)，則后半段的DFT值表達(dá)式：

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值表達(dá)式獲得，中間還利用到WN(N2?k)N?WN2Wk得到后半段的X[k]值表達(dá)式??W，k為：X[k]?X1[k]?WNkX2[k]（k=0,1,…,N/2-1）。

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

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

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

比如，一個(gè)N = 8點(diǎn)的FFT運(yùn)算按照這種方法來(lái)計(jì)算FFT可以用下面的流程圖來(lái)表示：

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按頻率抽?。―IF）的FFT算法思想：

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

在把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因?yàn)閃2N??1，則有WX[k]???(?1)，所以：

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

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

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

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

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通過(guò)上面的推導(dǎo)，X[k]的偶數(shù)點(diǎn)值X[2r]和奇數(shù)點(diǎn)值X[2r?1]分別可以由組合而成的N/2點(diǎn)的序列來(lái)求得，其中偶數(shù)點(diǎn)值X[2r]為輸入x[n]的前半段和后半段之和序列的N/2點(diǎn)DFT值，奇數(shù)點(diǎn)值X[2r?1]為輸入x[n]的前半段和后半段之差再與WN相乘序列的N/2點(diǎn)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

這樣，也可以用兩個(gè)N/2點(diǎn)DFT來(lái)組合成一個(gè)N點(diǎn)DFT，組合過(guò)程如下圖所示：

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

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

二．在FFT計(jì)算中使用到的MATLAB命令：

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

四、實(shí)驗(yàn)?zāi)康模?/p>

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

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

五、實(shí)驗(yàn)內(nèi)容：

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

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

六、實(shí)驗(yàn)器材（設(shè)備、元器件）：

安裝MATLAB軟件的PC機(jī)一臺(tái)，DSP實(shí)驗(yàn)演示系統(tǒng)一套。

七、實(shí)驗(yàn)步驟：

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

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

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

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

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

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

程序：（1）對(duì)要求的2種序列直接進(jìn)行DFT計(jì)算的程序

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

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

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

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

九、實(shí)驗(yàn)結(jié)論：

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

十一、對(duì)本實(shí)驗(yàn)過(guò)程及方法、手段的改進(jìn)建議：