第一篇:《數(shù)字信號處理原理及實現(xiàn)》課程小結(jié)
時間過得好快,轉(zhuǎn)眼半學(xué)期結(jié)束了。這半學(xué)期數(shù)字信號的學(xué)習(xí)讓我受益匪淺。前兩章和信號與線性系統(tǒng)相關(guān),介紹了離散時間信號與系統(tǒng)的時域分析方法最深刻的是采樣,時域采樣定理與采樣恢復(fù)和采樣內(nèi)插公式。第二章介紹了離散時間信號與系統(tǒng)的頻域分析,DTFT與Z變換,系統(tǒng)函數(shù)的零極點分布。第三章主要講了離散傅里葉變換DFT及其性質(zhì),和頻域采樣定理。第四章介紹了傅里葉快速變換FFT,熟悉了其原理特點及方法。五六兩章分別介紹了IIR和FIR濾波器,知道了IIR的脈沖響應(yīng)不變法與雙線性變換法及其優(yōu)缺點,并學(xué)會其MATLAB應(yīng)用設(shè)計濾波器,F(xiàn)IR的窗函數(shù)法與頻域采樣法設(shè)計濾波器及其MATLAB實現(xiàn)。第七章主要介紹了IIR和FIR濾波器的基本網(wǎng)絡(luò)結(jié)構(gòu),通過老師上課習(xí)題的練習(xí)基本掌握了其結(jié)構(gòu)圖的畫法。
先說說對課程的建議吧,張曉光老師是個很負(fù)責(zé)講課思路也很清晰的老師,知道從學(xué)生的角度來講問題,根據(jù)學(xué)生的反應(yīng)來調(diào)整課程進(jìn)度。我們都很喜歡這樣的老師,老師開新課之前總是列提綱復(fù)習(xí)上節(jié)課講的知識,每章結(jié)束都根據(jù)章節(jié)的重要性開一節(jié)總結(jié)課,這種方式個人覺得很好,希望老師堅持。但是,感覺老師講題講的不是很多,或許是課時原因,但我覺得每章結(jié)束后開一節(jié)例題課,把知識點融進(jìn)去,畢竟大學(xué)生現(xiàn)在做題比較少,這樣強(qiáng)制一下效果會更好。這次考試的試題覺得有不少都見過,有的是課后題,但做起來還是有點吃力,應(yīng)該就是習(xí)題練的少,計算跟不上去。至于教材,我覺得編的很好,每章都有相關(guān)的MATLAB編程方法,在原理講清之后就來實踐,免去了學(xué)生盲目做實驗,提高了效率。還有就是老師也很重視實驗,總是把相關(guān)的MATLAB語句語義講解清楚,這樣我們在編程序時也就相對容易點。但好像老師講程序時都注重程序的意思了,希望老師以后再講程序時把它先部分后整體,就是在講完程序意思后把程序設(shè)計思路或框架結(jié)構(gòu),及各部分要實現(xiàn)什么再講講,這樣有助于學(xué)生設(shè)計時設(shè)計思路更清晰。再說說考試,老師分卷面成績和實驗成績及平時成績,將實驗單獨考試,可見對實驗的重視,也說明MATLAB的重要性,這樣確實提高了學(xué)生的重視心理,雖然實驗做完了,但做完50道題并看完相關(guān)講解,我又收獲了不少,理清了設(shè)計方法與思路,所以我覺的考試方式還是挺不錯的,鍛煉了我們各方面的知識。
數(shù)字信號課程結(jié)束了,真希望您還能教我們別的課。
小組成員:陳文斌、李亞偉、王猛、汪子雄、吳官寶
第二篇:MATLAB實現(xiàn)數(shù)字信號處理
數(shù)字信號處理
說 明 書
目錄
一.摘要…………………………………3 二.課程設(shè)計目的………………………3 三.設(shè)計內(nèi)容……………………………3 四.設(shè)計原理……………………………4 4.1.語音信號的采集…………………………….4 4.2.濾波器……………………………………….4 4.21.IIR濾波器原理…………………………………….4 4.22.FIR濾波器原理………………………………………5 五.設(shè)計步驟……………………………6 5.1錄制女音………………………………………6 5.2采樣語音信號并畫出時域波形和頻譜圖……7 5.3采用雙線性變換法設(shè)計IIR濾波器…………10 5.4窗函數(shù)法設(shè)計FFR濾波器………………......12 5.5用IIR濾波器對信號進(jìn)行濾波………………14 5.6用FIR濾波器對信號進(jìn)行濾波………………16 5.7男女聲語音信號頻譜特點分析………………19 5.8有背景噪聲的信號分析………………………20 六.心得體會…………………………….22 七.參考文獻(xiàn)…………………………….23
一.摘要:
這次課程設(shè)計的主要目的是綜合運用本課程的理論知識進(jìn)行頻譜分析以及濾波器設(shè)計,通過理論推導(dǎo)得出相應(yīng)結(jié)論,并利用MATLAB或者DSP開發(fā)系統(tǒng)作為工具進(jìn)行實現(xiàn),從而復(fù)習(xí)鞏固課堂所學(xué)的理論知識,提高對所學(xué)知識的綜合應(yīng)用能力,并從實踐上初步實現(xiàn)對數(shù)字信號的處理。通過對聲音的采樣,將聲音采樣后的頻譜與濾波。
MATLAB全稱是Matrix Laboratory,是一種功能強(qiáng)大、效率高、交互性好的數(shù)值和可視化計算機(jī)高級語言,它將數(shù)值分析、矩陣運算、信號處理和圖形顯示有機(jī)地融合為一體,形成了一個極其方便、用戶界面友好的操作環(huán)境。經(jīng)過多年的發(fā)展,已經(jīng)發(fā)展成為一種功能全面的軟件,幾乎可以解決科學(xué)計算中所有問題。MATLAB軟件還提供了非常廣泛和靈活的用于處理數(shù)據(jù)集的數(shù)組運算功能。
在本次課程設(shè)計中,主要通過MATLAB來編程對語音信號處理與濾波,設(shè)計濾波器來處理數(shù)字信號并對其進(jìn)行分析。
二.課程設(shè)計目的:
綜合運用本課程的理論知識進(jìn)行頻譜分析以及濾波器設(shè)計,通過理論推導(dǎo)得出相應(yīng)結(jié)論,并利用MATLAB作為工具進(jìn)行實現(xiàn),從而復(fù)習(xí)鞏固課堂所學(xué)的理論知識,提高對所學(xué)知識的綜合應(yīng)用能力,并從實踐上初步實現(xiàn)對數(shù)字信號的處理。
三.設(shè)計內(nèi)容:
內(nèi)容:錄制一段個人自己的語音信號,并對錄制的信號進(jìn)行采樣;畫出采樣后語音信號的時域波形和頻譜圖;給定濾波器的性能指標(biāo),采用窗函數(shù) 法和雙線性變換法設(shè)計濾波器,并畫出濾波器的頻率響應(yīng);然后用自己設(shè)計的濾波器對采集的信號進(jìn)行濾波,畫出濾波后信號的時域波形和頻譜,并對濾波前后的信號進(jìn)行對比,分析信號的變化;回放語音信號;換一個與你性別相異的人錄制同樣一段語音內(nèi)容,分析兩段內(nèi)容相同的語音信號頻譜之間有什么特點;再錄制一段同樣長時間的背景噪聲疊加到你的語音信號中,分析疊加前后信號頻譜的變化,設(shè)計一個合適的濾波器,能夠把該噪聲濾除。
四.設(shè)計原理:
4.1.語音信號的采集
熟悉并掌握MATLAB中有關(guān)聲音(wave)錄制、播放、存儲和讀取的函數(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ù)的存儲格式,用字符串指定,可以是:‘double’、‘single’、’int16’、‘int8’其中只有int8是采用8位精度進(jìn)行采樣,其它三種都是16位采樣結(jié)果轉(zhuǎn)換為指定的MATLAB數(shù)據(jù);
b:wavplay(y,fs);利用系統(tǒng)音頻輸出設(shè)備播放,以fs為播放頻率,播放語音信號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ù)字濾波器在時域上模擬濾波器,但是它們的缺點是產(chǎn)生頻率響應(yīng)的混疊失真,這是由于從s平面到z平面是多值的映射關(guān)系所造成的。
雙線性變換法是使數(shù)字濾波器的頻率響應(yīng)與模擬濾波器的頻率響應(yīng)相似的一種變換方法。為了克服多值映射這一缺點,我們首先把整個s平面壓縮變換到某一中介的s1平面的一條橫帶里,再通過變換關(guān)系將此橫帶變換到整個z平面上去,這樣就使得s平面與z平面是一一對應(yīng)的關(guān)系,消除了多值變換性,也 就消除了頻譜混疊現(xiàn)象。
雙線性法設(shè)計IIR數(shù)字濾波器的步驟:
1)將數(shù)字濾波器的頻率指標(biāo){ ?k}由Wk=(2/T)*tan(wk),轉(zhuǎn)換為模擬濾波器的頻率指標(biāo){?k}.2)由模擬濾波器的指標(biāo)設(shè)計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)絡(luò)進(jìn)行相位校正。而有限長單位沖激響應(yīng)(FIR)數(shù)字濾波器就可以做成具有嚴(yán)格的線性相位,同時又可以具有任意的幅度特性。
由于FIR系統(tǒng)的沖激響應(yīng)就是其系統(tǒng)函數(shù)各次項的系數(shù),所以設(shè)計FIR濾波器的方法之一可以從時域出發(fā),截取有限長的一段沖激響應(yīng)作為H(z)的系數(shù),沖激響應(yīng)長度N就是系統(tǒng)函數(shù)H(z)的階數(shù)。只要N足夠長,截取的方法合理,總能滿足頻域的要求。這種時域設(shè)計、頻域檢驗的方法一般要反復(fù)幾個回合,不像IIR DF設(shè)計靠解析公式一次計算成功。給出的理想濾波器頻率響應(yīng)是,它是w的周期函數(shù),周期
由傅立葉反變換導(dǎo)出,即
hd(n)?1Hd(ejw)ejwndw?2?,再將hd(n)與窗函數(shù),因此可展開成傅氏級數(shù)w(n)相乘就可以得到h(n)。、的計算可采用傅氏變換的現(xiàn)成公式和程序,窗函數(shù)也是現(xiàn)成的。但整個設(shè)計過程不能一次完成,因為窗口類型和大小的選擇沒有解析公式可一次算,整個設(shè)計可用計算機(jī)編程來做。
窗函數(shù)的傅式變換W(ejω)的主瓣決定了H(ejω)過渡帶寬。W(ejω)的旁瓣大小和多少決定了H(ejω)在通帶和阻帶范圍內(nèi)波動幅度,常用的幾種窗函數(shù)有:
矩形窗
w(n)=RN(n);
Hanning窗
;
Hamming窗
;
Blackmen窗
;
Kaiser窗。
式中Io(x)為零階貝塞爾函數(shù)。
五.設(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('原女音信號的時域波形');
%給圖形加注標(biāo)簽說明 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)));
%繪制原始語音信號的幅頻響應(yīng)圖 title('FFT變換后聲音的頻譜特性');
%給圖形加注標(biāo)簽說明 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('截取后女音信號的時域波形');
%給圖形加注標(biāo)簽說明 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)));
%繪制原始語音信號的幅頻響應(yīng)圖 title('FFT變換后聲音的頻譜特性');
%給圖形加注標(biāo)簽說明 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采用雙線性變換法設(shè)計IIR濾波器:
人的聲音頻率一般在(1~~4)kHZ之間,則我們只需要設(shè)計一個帶通濾波器即可濾去聲音頻帶以外的無用噪聲,得到比較清晰的聲音。根據(jù)聲音的頻譜圖分析,設(shè)計一個帶通濾波器性能指標(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);
%估計契比雪夫I型濾波器階數(shù) [b,a]=cheby1(N,Ap,wn);
%N指定濾波器階數(shù),wn歸一化
截 %止頻率,Ap通帶波動
[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è)計FFR濾波器
線性相位FIR濾波器通常采用窗函數(shù)法設(shè)計。窗函數(shù)法設(shè)計FIR濾波器的基本思想是:根據(jù)給定的濾波器技術(shù)指標(biāo),選擇濾波器長度N和窗函數(shù)ω(n),使其具有最窄寬度的主瓣和最小的旁瓣。其核心是從給定的頻率特性,通過加窗確定有限長單位脈沖響應(yīng)序列h(n)。工程中常用的窗函數(shù)共有6種,即矩形窗、巴特利特(Bartlett)窗、漢寧(Hanning)窗、漢明(Hamming)窗、布萊克曼(Blackman)窗和凱澤(Kaiser)窗。
這次設(shè)計我采用的是布萊克曼來設(shè)計給定數(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);
%求濾波器頻率響應(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 = 點0 到 N-1之間的理想脈沖響應(yīng) % wc = 截止頻率(弧度)% N = 理想濾波器的長度
tao =(N-1)/2;
% 理想脈沖響應(yīng)的對稱中心位置 n = [0:(N-1)];
% 設(shè)定脈沖響應(yīng)長度 m = n-tao + eps;
% 加一個小數(shù)以避免零作除數(shù)
hd = sin(wc*m)./(pi*m);
% 理想脈沖響應(yīng)
布拉克曼窗1振幅0.500406080時間FIR帶通濾波器幅頻特性500振幅/分貝20FIR帶通濾波器幅頻特性1.51振幅-50-100-15000.5?/10.5000.5?/?1
5.5用IIR濾波器對信號進(jìn)行濾波
用自己設(shè)計的IIR濾波器分別對采集的信號進(jìn)行濾波,在Matlab中,IIR濾波器利用函數(shù)filter對信號進(jìn)行濾波。程序如下: [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)));
%繪制原始語音信號的幅頻響應(yīng)圖 title('濾波前女音的頻譜特性');
%給圖形加注標(biāo)簽說明 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)));
%繪制原始語音信號的幅頻響應(yīng)圖 title('l濾波后女音的頻譜特性');
%給圖形加注標(biāo)簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;
濾波前女音的頻譜特性300幅值/A***030004000頻率/Hz濾波后女音的頻譜特性500060006040幅值/A***0頻率/Hz400050006000
結(jié)果分析:
由上面濾波前后的頻譜圖可看出,濾波器濾除了小于1000Hz和大于3400Hz的頻譜成分。回放語音信號,由于低頻和高頻成分被濾除,聲音變得較低沉。
5.6用FIR濾波器對信號進(jìn)行濾波
用自己設(shè)計的FIR濾波器分別對采集的信號進(jìn)行濾波,在Matlab中,FIR濾波器利用函數(shù)fftfilt對信號進(jìn)行濾波 程序如下:
[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)));
%繪制原始語音信號的幅頻響應(yīng)圖 title('濾波前女音的頻譜特性');
%給圖形加注標(biāo)簽說明 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)));
%繪制原始語音信號的幅頻響應(yīng)圖 title('l濾波后女音的頻譜特性');
%給圖形加注標(biāo)簽說明 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列,指定第一個圖%繪制原始語音信號的幅頻響應(yīng)圖
%給圖形加注標(biāo)簽說明
%添加網(wǎng)格
%讀取聲音
%計數(shù)讀取信號的點數(shù)
%信號的時域采樣點
%采樣間隔
%取信號的一半
%對信號做N點FFT變換
%頻域采樣點
%把畫圖區(qū)域劃分為2行1列,指定第二個圖%繪制原始語音信號的幅頻響應(yīng)圖
%給圖形加注標(biāo)簽說明
%添加網(wǎng)格
axis([0 6000 0 300]);
%改變橫縱坐標(biāo)便于比較頻譜圖
FFT變換后女音的頻譜特性300200幅值/A***00頻率/?FFT變換后男音的頻譜特性***200幅值/A***00頻率/?400050006000
結(jié)果分析:
通過比較上面女音頻譜圖和男音頻譜圖可知,男音的頻譜集中在低頻部分,高頻成分底,譜線較平滑,聲音聽起來低沉。5.8有背景噪聲的信號分析
從硬盤中把一段噪聲(頻譜能量集中在某個小范圍內(nèi))疊加到語音信號中,分析疊加前后信號頻譜的變化,設(shè)計一個合適的濾波器,能夠把該噪聲濾除; 程序如下:
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)));
%繪制原始語音信號的幅頻響應(yīng)圖 title('加噪聲后聲音的頻譜特性');
%給圖形加注標(biāo)簽說明 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)));
%繪制原始語音信號的幅頻響應(yīng)圖 title('濾除噪聲后聲音的頻譜特性');
%給圖形加注標(biāo)簽說明 xlabel('頻率/Hz');ylabel('幅值/A');grid;%添加網(wǎng)格
加噪聲后聲音的頻譜特性3000幅值/A***0008000頻率/Hz濾除噪聲后聲音的頻譜特性***030幅值/A***000頻率/Hz80001000012000
結(jié)果分析
觀察加噪聲后聲音的頻譜圖可知,噪音頻率主要在4000Hz處,只要我們設(shè)計一個,濾波器濾除大概在4000Hz的頻譜即可,回放濾波后的語音信號,可證噪音基本濾除。
六.心得體會:
通過這次課程設(shè)計,讓我對MATLAB的基本應(yīng)用有了更深的了解,還有數(shù)字信號處理在MATLAB中的一些函數(shù)的用法。通過理論推導(dǎo)得出相應(yīng)結(jié)論,并利用MATLAB作為工具進(jìn)行實現(xiàn),從而復(fù)習(xí)鞏固課堂所學(xué)的理論知識,提高對所學(xué)知識的綜合應(yīng)用能力,并從實踐上初步實現(xiàn)對數(shù)字信號的處理。
在這次實驗中,也遇到了很多問題,比如畫信號頻譜的分貝圖時(20*log10(abs(y)))指數(shù)為零時的處理。濾波器的設(shè)計也花了好大的功夫,剛開始不會設(shè)計參數(shù),一頭霧水,通過同學(xué)的指導(dǎo)和討論,得知通過觀察信號的頻譜圖,看噪音頻率集中在那一部分,設(shè)計濾波器把其濾除即可。可反復(fù)設(shè)置參數(shù)直到濾波后語音信號的效果好為止。
七.參考文獻(xiàn):
(1)《MATLAB LabVIEW SystemView》翁劍楓 葉志前 編著, 機(jī)械工業(yè)出版社;
(2)《MATLAB及在電子信息課程中的應(yīng)用》陳懷琛 吳大正 高西全編著,電子工業(yè)出版社;
(3)《MATLAB在數(shù)字信號處理中的應(yīng)用》(弟2版)薛年喜 編著,清華大學(xué)出版社;
(4)《MATLAB擴(kuò)展編程》何強(qiáng) 何英
編著,清華大學(xué)出版社;(5)《MATLAB7簡明教程》吳清 曹輝林 編著,清華大學(xué)出版社;(6)MATLAB5.3精要.編程及高級應(yīng)用》程衛(wèi)國 馮峰 王雪梅 劉藝 編著,機(jī)械工程出版社。
第三篇:數(shù)字信號處理課程總結(jié)(推薦)
數(shù)字信號處理課程總結(jié)
信息09-1班 陳啟祥 金三山 趙大鵬 劉恒
進(jìn)入大三,各種專業(yè)課程的學(xué)習(xí)陸續(xù)展開,我們也在本學(xué)期進(jìn)行了數(shù)字信號處理這門課程的學(xué)習(xí)。
作為信心工程專業(yè)的核心課程之一,數(shù)字信號處理的重要性是顯而易見的。在近九周的學(xué)習(xí)過程中,我們學(xué)習(xí)了離散時間信號與系統(tǒng)的時域及頻域分析、離散傅里葉變換、快速傅里葉變換、IIR及FIR數(shù)字濾波器的設(shè)計及結(jié)構(gòu)等相關(guān)知識,并且在實驗課上通過MATLAB進(jìn)行了相關(guān)的探究與實踐。總體來說,通過這一系列的學(xué)習(xí)與實踐,我們對數(shù)字信號處理的有關(guān)知識和基礎(chǔ)理論已經(jīng)有了初步的認(rèn)知與了解,這對于我們今后進(jìn)一步的學(xué)習(xí)深造或參加實際工作都是重要的基礎(chǔ)。
具體到這門課程的學(xué)習(xí),應(yīng)當(dāng)說是有一定的難度的。課本所介紹的相關(guān)知識理論性很強(qiáng),并且與差分方程、離散傅里葉級數(shù)、傅里葉變換、Z變換等數(shù)學(xué)工具聯(lián)系十分緊密,所以要真正理解課本上的相關(guān)理論,除了認(rèn)真聆聽老師的講解,還必須要花費大量時間仔細(xì)研讀課本,并認(rèn)真、獨立地完成課后習(xí)題。總之,理論性強(qiáng)、不好理解是許多同學(xué)對數(shù)字信號處理這門課程的學(xué)習(xí)感受。
另外,必須要說MATLAB實驗課程的開設(shè)是十分必要的。首先,MATLAB直觀、簡潔的操作界面對于我們真正理解課堂上學(xué)來的理論知識幫助很大;其次,運用MATLAB進(jìn)行實踐探究,也使我們真正意識到,在信息化的今天,研究數(shù)字信號離不開計算機(jī)及相關(guān)專業(yè)軟件的幫助,計算機(jī)及軟件技術(shù)的發(fā)展,是今日推動信息技術(shù)發(fā)展的核心動力;最后,作為信息工程專業(yè)的學(xué)生,在許多學(xué)習(xí)與實踐領(lǐng)域需要運用MATLAB這樣一個強(qiáng)大工具,MATLAB實驗課程的開設(shè),鍛煉了我們的實踐能力,也為我們今后在其他領(lǐng)域運用MATLAB打下了基礎(chǔ)。
課程的結(jié)束、考試的結(jié)束不代表學(xué)習(xí)的結(jié)束,數(shù)字信號處理作為我們專業(yè)的基礎(chǔ)之一,是不應(yīng)當(dāng)被我們拋之腦后的。
最后感謝老師這幾周來的教誨與指導(dǎo),謝謝老師!
2012年5月7日
第四篇:數(shù)字信號處理課程總結(jié)(全)
數(shù)字信號處理課程總結(jié)
以下圖為線索連接本門課程的內(nèi)容:
xa(t)數(shù)字信號前置濾波器A/D變換器處理器D/A變換器AF(濾去高頻成分)ya(t)x(n)
一、時域分析
1. 信號
? 信號:模擬信號、離散信號、數(shù)字信號(各種信號的表示及關(guān)系)? 序列運算:加、減、乘、除、反褶、卷積 ? 序列的周期性:抓定義
njwna、e?(n)(可表征任何序列)cos(wn??)u(n)、? 典型序列:、、RN(n)、?x(n)??x(m)?(n?m)
m???特殊序列:h(n)2. 系統(tǒng)
? 系統(tǒng)的表示符號h(n)? 系統(tǒng)的分類:y(n)?T[x(n)]
線性:T[ax1(n)?bx2(n)]?aT[x1(n)]?bT[x2(n)] 移不變:若y(n)?T[x(n)],則y(n?m)?T[x(n?m)] 因果:y(n)與什么時刻的輸入有關(guān) 穩(wěn)定:有界輸入產(chǎn)生有界輸出
? 常用系統(tǒng):線性移不變因果穩(wěn)定系統(tǒng) ? 判斷系統(tǒng)的因果性、穩(wěn)定性方法 ? 線性移不變系統(tǒng)的表征方法:
線性卷積:y(n)?x(n)*h(n)
NMk差分方程: y(n)??ak?1y(n?k)??bk?0kx(n?k)3. 序列信號如何得來?
xa(t)x(n)抽樣
? 抽樣定理:讓x(n)能代表xa(t)? 抽樣后頻譜發(fā)生的變化? ? 如何由x(n)恢復(fù)xa(t)?
?sin[xa(mT)?T(t?mT)]
xa(t)=?m????T
(t?mT)
二、復(fù)頻域分析(Z變換)
時域分析信號和系統(tǒng)都比較復(fù)雜,頻域可以將差分方程變換為代數(shù)方程而使分析簡化。A. 信號 1.求z變換
?定義:x(n)?X(z)??x(n)zn????n
收斂域:X(z)是z的函數(shù),z是復(fù)變量,有模和幅角。要其解析,則z不能取讓X(z)無窮大的值,因此z的取值有限制,它與x(n)的種類一一對應(yīng)。
? x(n)為有限長序列,則X(z)是z的多項式,所以X(z)在z=0或∞時可能會有∞,所以z的取值為:0?z??;
? x(n)為左邊序列,0?z?Rx?,z能否取0看具體情況;
? x(n)為右邊序列,Rx??z??,z能否取∞看具體情況(因果序列); ? x(n)為雙邊序列,Rx??z?Rx? 2.求z反變換:已知X(z)求x(n)
? 留數(shù)法
? 部分分式法(常用):記住常用序列的X(z),注意左右序列區(qū)別。? 長除法:注意左右序列 3.z變換的性質(zhì):
? 由x(n)得到X(z),則由x(n?m)?z?mX(z),移位性; ? 初值終值定理:求x(0)和x(?);
? 時域卷積和定理:y(n)?x(n)*h(n)?Y(z)?X(z)H(z); ? 復(fù)卷積定理:時域的乘積對應(yīng)復(fù)頻域的卷積; ? 帕塞瓦定理:能量守恒
?
?n???x(n)2?12?????X(ejw)dw2
4.序列的傅里葉變換
?公式:X(ejw)??x(n)en????jwn
x(n)?12?????X(ej?)ej?nd?
注意:X(ejw)的特點:連續(xù)、周期性;X(ejw)與X(z)的關(guān)系 B. 系統(tǒng)
由h(n)?H(z),系統(tǒng)函數(shù),可以用來表征系統(tǒng)。
? H(z)的求法:h(n)?H(z);H(z)=Y(z)/X(z); ? 利用H(z)判斷線性移不變系統(tǒng)的因果性和穩(wěn)定性 ? 利用差分方程列出對應(yīng)的代數(shù)方程
MNMy(n)??ak?1y(n?k)?k?bk?0x(n?k)?kY(z)X(z)?b?k?0Nkz?k
k1??ak?1z?k? 系統(tǒng)頻率響應(yīng)H(ejw):以2?為周期的?的連續(xù)函數(shù)
?
H(e)?jw?h(n)en?????jwn
H(e?jw)??h(n)en???jwn,當(dāng)h(n)為實序列時,則有H(ejw)=H*(e?jw)
三、頻域分析
根據(jù)時間域和頻域自變量的特征,有幾種不同的傅里葉變換對
? 時間連續(xù),非周期?頻域連續(xù)(由時域的非周期造成),非周期(由時域的連續(xù)造成); ?X(j?)??x(t)e????j?tdt
x(t)?12????X(j?)ej?td?
? 時間連續(xù),周期?頻域離散,非周期
X(jk?0)?1T0T0/2?x(t)e?jk?0tdt
?T0/2x(t)??X(jk?0)ejk?0t
? 時間離散,非周期?頻域連續(xù),周期
?
X(e)?jw?x(n)en????jwn
x(n)?12?????X(ej?)ej?nd?,w??T(數(shù)字頻率與模擬頻率的關(guān)系式)
? 時間離散,周期?頻域離散,周期
~X(k)?N?1?n?0~x(n)e?j2?Nkn?~?x(n)W
knNn?0N?11~x(n)?NN?1?n?0~X(k)ej2?Nkn?1NN?1?n?0~?knX(k)WN
? 本章重點是第四種傅里葉變換-----DFS ? 注意:
x(n)和X(k)都是以N為周期的周期序列; 1)~x(n)和X(k)的定義域都為(??,?)
2)盡管只是對有限項進(jìn)行求和,但~;
~~~例如:k?0時,X(0)?N?1?x(n)
n?0~~k?1時,X(1)?N?1?n?0~x(n)e?j2?Nn
2?NNnN?1~k?N時,X(N)?N?1?n?0?j~x(n)e??n?02?N~~x(n)=X(0)
~k?N?1時,X(N?1)?N?1?n?0~x(n)e?j(N?1)n~?X(1)
x(n)也有類似的結(jié)果。x(n)和X(k)一
同理也可看到~可見在一個周期內(nèi),~~一對應(yīng)。
?? 比較X(e)?jw?x(n)en????jwn~和X(k)?N?1?n?0~x(n)e?j2?Nkn?~?x(n)W,當(dāng)x(n)knNn?0N?1x(n)的一個周期內(nèi)有定義時,即x(n)=~x(n),0?n?N?1,則在只在~??N?12?Nj2?Nk時,X(ejw)?X(k)。
?1,k?r?? 0,k?r?~? ?en?0(k?r)nx(n)和X(k)的每個周期值都只是其主值區(qū)間的周期延拓,所以求和? 因為~~在任一個周期內(nèi)結(jié)果都一樣。
? DFT:有限長序列x(n)只有有限個值,若也想用頻域方法分析,它只屬于序列的傅里葉變換,但序列的傅氏變換為連續(xù)函數(shù),所以為方便計算機(jī)處理,也希望能像DFS一樣,兩個域都離散。將x(n)想象成一個周期x(n)的一個周期,然后做DFS,即 序列~
~X(k)?N?1?n?0~x(n)e?j2?NknN?1??n?0x(n)e?j2?Nkn
x(n)只有x(n),不是真正的周期序列,但因為求和只需N注意:實際上~個獨立的值,所以可以用這個公式。同時,盡管x(n)只有N個值,但依上式求出的X(k)還是以N為周期的周期序列,其中也只有N個值獨立,這樣將~X(k)規(guī)定在一個周期內(nèi)取值,成為一個有限長序列,則會引出
N?1?j2?Nkn~DFT X(k)??x(n)en?0RN(k)
x(n)?1NN?1?n?0X(k)ej2?NknRN(n)
比較:三種移位:線性移位、周期移位、圓周移位
三種卷積和:線性卷積、周期卷積、圓周卷積
重點:1)DFT的理論意義,在什么情況下線性卷積=圓周卷積 2)頻域采樣定理:掌握內(nèi)容,了解恢復(fù)
3)用DFT計算模擬信號時可能出現(xiàn)的幾個問題,各種問題怎樣引起?
混疊失真、頻譜泄漏、柵欄效應(yīng)
? FFT:為提高計算速度的一種算法
1)常用兩種方法:按時間抽取基2算法和按頻率抽取基2算法,各自的原理、特點是什么,能自行推導(dǎo)出N小于等于8的運算流圖。2)比較FFT和DFT的運算量; 3)比較DIT和DIF的區(qū)別。
四、數(shù)字濾波器(DF)
一個離散時間系統(tǒng)可以用h(n)、H(z)、差分方程和H(ejw)來表征。問題:
1、各種DF的結(jié)構(gòu)
2、如何設(shè)計滿足要求指標(biāo)的DF?
3、如何實現(xiàn)設(shè)計的DF?
A. 設(shè)計IIR DF,借助AF來設(shè)計,然后經(jīng)S---Z的變換即可得到。
1)脈沖響應(yīng)不變法:思路、特點 2)雙線性變換法:思路、特點、預(yù)畸變 3)模擬濾波器的幅度函數(shù)的設(shè)計 B. 設(shè)計FIR DF 1)線性相位如何得到?條件是什么?各種情況下的特點。2)窗函數(shù)設(shè)計法:步驟、特點 3)頻率抽樣法:步驟、特點 C. 實現(xiàn)DF
M?a
標(biāo)準(zhǔn)形式:H(z)?k?0Nkz?k
bkz?k1??k?1
第五篇:數(shù)字信號處理實驗-FFT的實現(xiàn)
實
驗
報
告
學(xué)生姓名:
學(xué) 號:
指導(dǎo)教師:
一、實驗室名稱:數(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]為復(fù)數(shù)序列的情況下,完全直接運算N點DFT需要(N-1)2次復(fù)數(shù)乘法和N(N-1)次加法。因此,對于一些相當(dāng)大的N值(如1024)來說,直接計算它的DFT所作的計算量是很大的。
FFT的基本思想在于,將原有的N點序列分成兩個較短的序列,這些序列的DFT可以很簡單的組合起來得到原序列的DFT。例如,若N為偶數(shù),將原有的N
22點序列分成兩個(N/2)點序列,那么計算N點DFT將只需要約[(N/2)·2]=N/2次復(fù)數(shù)乘法。即比直接計算少作一半乘法。因子(N/2)2表示直接計算(N/2)點DFT所需要的乘法次數(shù),而乘數(shù)2代表必須完成兩個DFT。上述處理方法可以反復(fù)使用,即(N/2)點的DFT計算也可以化成兩個(N/4)點的DFT(假定N/2為偶數(shù)),從而又少作一半的乘法。這樣一級一級的劃分下去一直到最后就劃分成兩點的FFT運算的情況。
3.基2按時間抽取(DIT)的FFT算法思想:
設(shè)序列長度為N?2L,L為整數(shù)(如果序列長度不滿足此條件,通過在后面補(bǔ)零讓其滿足)。
將長度為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]的后半段本算法推導(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)。
這樣,通過計算兩個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算法思想:
設(shè)序列長度為N?2L,L為整數(shù)(如果序列長度不滿足此條件,通過在后面補(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因為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)推導(dǎo)過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通過上面的推導(dǎo),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 四、實驗?zāi)康模?/p> 離散傅氏變換(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。 六、實驗器材(設(shè)備、元器件): 安裝MATLAB軟件的PC機(jī)一臺,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é)果進(jìn)行對比。 (5)(拓展要求)嘗試對其他快速傅立葉變換算法(如Goertzel算法)進(jìn)行MATLAB編程實現(xiàn),并用它來計算要求的序列的DFT值。并與前面的結(jié)果進(jìn)行對比。 (6)(拓展要求)在提供的DSP實驗板上演示要求的2種序列的FFT算法(基2-DIT),用示波器觀察實際計算出來的頻譜結(jié)果,并與理論結(jié)果對比。 八、實驗數(shù)據(jù)及結(jié)果分析: 程序:(1)對要求的2種序列直接進(jìn)行DFT計算的程序 (2)對要求的2種序列進(jìn)行基2-DIT和基2-DIF FFT算法程序(3)對要求的2種序列用MATLAB中提供的FFT函數(shù)進(jìn)行計算的程序 結(jié)果:(1)對2種要求的序列直接進(jìn)行DFT計算的頻域波形 (2)對2種要求的序列進(jìn)行基2-DIT和基2-DIF FFT算法頻域波形(3)對2種要求的序列用MATLAB中提供的FFT函數(shù)計算的頻域波形。(4)(拓展要求)分析利用上面的方法畫出的信號頻譜與理論計算出來的頻譜之間的差異,并解釋這種差異。 (5)(拓展要求)保持序列點數(shù)不變,改變DFT計算點數(shù)(變?yōu)?024點),觀察頻譜的變化,并分析這種變化,由此討論如何提高頻譜分辨力的問題。 九、實驗結(jié)論: 十、總結(jié)及心得體會: 十一、對本實驗過程及方法、手段的改進(jìn)建議: