第一篇：數字信號處理實驗二

北京信息科技大學

數字信號處理實驗報告

題 目：

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

學 院： 信息與通信工程學院 專 業： 通信工程專業 姓 名： 班 級：

學 號： 指導老師：

實驗目的

1、熟悉IIR數字濾波器的設計原理與方法。

2、掌握數字濾波器的計算機軟件實現方法。

3、通過觀察對實際心電圖信號的濾波作用，學習數字濾波器在實際中的應用。

實驗儀器及材料

計算機，MATLAB軟件

實驗內容及要求

1．設計巴特沃斯低通數字濾波器對人體心電信號進行濾波

(1)人體心電圖信號在測量過程中會受到工業高頻干擾，所以必須經過低通濾波處理，才能作為判斷心臟功能的有用信息。以下為一個實際心電圖信號采樣序列x(n)，其中存在高頻干擾，抽樣周期Ts=1秒。在實驗中，以x(n)作為輸入序列，濾除其中干擾成分。

x(n)=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0] 對序列x(n)用FFT做頻譜分析，生成x(n)的頻譜圖。

(2)設計一個巴特沃斯低通IIR數字濾波器H(z)。

設計指標參數為：在通帶內頻率低于0.2π時，最大衰減小于1dB； 在阻帶內 [0.3π, π]頻率區間上，最小衰減大于15dB。

j?|H(e)|。寫出數字濾波器H(z)的表達式，畫出濾波器的幅頻響應曲線

(3)用所設計的濾波器對實際心電圖信號采樣序列x(n)進行濾波處理，編寫程序，求濾波后的序列y(n)，并分別畫出濾波前后的心電圖信號波形圖和頻譜圖。

y(n)= [0,0,0,0, 0,0,0,0,-0.14025,0.40279,-0.56085 ,0.33328,0.023981,-0.18809,0.11843,-0.1038,0.11576,-0.1225,0.099815 ,-0.13769 ,0.095249,-0.0070273,0.018867,0.090543,-0.11257,-0.070884 ,0.17676,-0.55407,0.24813,-0.34732,-0.30428,0.59426,-0.29574,-0.063869,0.34018,-0.73334,1.0293,-0.57107,-0.2461,0.83605,-0.83026,0.45459,0.011551,-0.25667,0.23896,-0.17361,0.20829,-0.28417,0.28765 ,-0.2035,0.02865,0.066164,0.077916,-0.36052, 0.53517,-0.5571]

源程序

clear all，clc

x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];%未經濾波的心電圖信號 L=length(x);l=0:L-1;

y=fft(x,L);Wp=0.2*pi;Ws=0.3*pi;Rp=1;Rs=15;

[N,Wn] = buttord(Wp,Ws,Rp,Rs,'s');[b,a] = butter(N,Wn,'s');[numa,dena]=impinvar(b,a,1);w=linspace(0,pi,1024);h=freqz(numa,dena,w);norm=max(abs(h));numa=numa/norm;[z,p]=tf2zp(b,a);figure(1)

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

xlabel('數字頻率');ylabel('幅度響應dB');figure(2)plot(w,abs(h));grid;

xlabel('數字頻率');

ylabel('幅度響應|H(e^(jw))|');figure(3)zplane(z,p);xx=filter(b,a,x);yy=fft(xx,L);figure(4)subplot(2,1,1)stem(l,x);

title('未經濾波的心電圖信號');xlabel('n');subplot(2,1,2)stem(l,xx);

title('經濾波之后的心電圖信號');xlabel('n');figure(5)subplot(2,1,1)plot(l,abs(y));

title('未經濾波的心電圖信號的頻譜');subplot(2,1,2)plot(l,abs(yy));

title('經濾波處理的心電圖信號的頻譜');

2．用help查看內部函數cheb1ord.m及cheby1.m，了解調用格式。

編程設計教材習題6-2，求模擬濾波器Ha(s)的表達式。

源程序

close all clear all clc

Wp=2*pi*3000;Rp=2;Ws=2*pi*12000;Rs=50;Fs=24000;

w=linspace(0,pi,1024);

[N,Wn]=cheb1ord(Wp,Ws,Rp,Rs,'s');e=sqrt(10^(Rp/10)-1);[b,a]=cheby1(N,e,Wn,'s')[numa,dena]=impinvar(b,a,Fs);h=freqz(numa,dena,w);norm=max(abs(h));

plot(w*Fs/pi,20*log10(abs(h)/norm))title('幅度響應')xlabel('頻率(Hz)')ylabel('幅度(dB)')grid

3．模擬濾波器的數字化

用內部函數impinvar及bilinear實現教材習題6-5，求數字濾波器H(z)的表達式。

源程序

close all clear all clc b1=[0 0 1];a1=[1 1 1];b2=[0 0 1];a2=[2 3 1];

[numa1,dena1]=impinvar(b1,a1,0.5)[numa2,dena2]=bilinear(b1,a1,0.5)[numa3,dena3]=impinvar(b2,a2,0.5)

[numa4,dena4]=bilinear(b2,a2,0.5)

本實驗所用的部分MATLAB函數

? L=length(x)：求序列x長度。

? y=fft(x,L)：將序列x(n)做L點快速傅立葉變換，結果賦給序列y(n)。? [n,Wn] = buttord(Wp,Ws,Rp,Rs,'s')：計算模擬Butterworth濾波器的最小階次n和截止頻率為Wn。

? [b,a] = butter(n,Wn,'s')：設計模擬截止頻率為Wn（rad/s）的n階 Butterworth低通濾波器，返回值為模擬濾波器的系數。

? y=filter(b,a,x): 將序列x(n)通過濾波器濾波后生成序列y(n)，濾波器的分母多項式系數構成a向量，分子多項式系數構成b向量。

? [BZ,AZ] = impinvar(B,A,Fs)：沖激響應不變法，返回值為數字濾波器的系數。? [BZ,AZ] = bilinear(B,A,fs)：雙線性變換，返回值為數字濾波器的系數。? [H w]=freqz(b,a)：由濾波器分母多項式系數構成的a向量和分子多項式系數構成的b向量求系統頻響。

截圖

實驗體會

這次的實驗讓這讓我看到了理論與實踐相結合的優勢與用處，讓我受益匪淺。我認識到了自己理論知識的不足，也認識到了我們學習的基礎知識究竟能運用于什么領域，如何運用。我們在老師的耐心指導下調試電路，直到得到要求的效果。讓我們在學習電路、信號等理論知識的同時，明白如何把這些應用于實際。

第二篇：數字信號處理實驗5

實驗五 FFT 算法的應用

1.進一步加深對離散信號 DFT 的理解。

2.運用其 FFT 算法解決實際問題。

1.微機。

2.Matlab 編程環境。

1.熟悉 Matlab 的編程環境和編程語言。

2.學習教材 P213-227，P249-263，掌握快速傅里葉變換（FFT）的原理。

1.實驗重點、難點、特點

快速傅里葉變換（FFT）的原理及應用。難點在 FFT 的應用以及 Matlab 編程中矩陣乘和數乘的應用。

2.實驗原理

一、實驗目的

二、實驗儀器設備

三、實驗學時 學時

四、預習要求

五、實驗特點及實驗原理簡介

1.已知 2N 點實數序列

N=64。理論計算 X(k)=DFT[x(n)]2N。

六、實驗內容及步驟

N=64;n=0:1:N-1;n2=0:1:2*N-1;k=0:1:2*N-1;xn1=cos(2*pi/N*7*2*n)+1/2*cos(2*pi/N*19*2*n);xn2=cos(2*pi/N*7*(2*n+1))+1/2*cos(2*pi/N*19*(2*n+1));wn=xn1+i*xn2;[wk,kk]=dft_my(wn,N);xk1=1/2*(wk+conj(wk(mod(-kk,N)+1)));xk2=1/2*(wk-conj(wk(mod(-kk,N)+1)))/i;XK1=fft(xn1,N);XK2=fft(xn2,N);%subplot(2,2,1);%stem(kk,abs(xk1));grid on;%subplot(2,2,3);%stem(kk,abs(XK1));grid on;%subplot(2,2,2);%stem(kk,abs(xk2));grid on;%subplot(2,2,4);%stem(kk,abs(XK2));grid on;

xk=[xk1+(exp(-1i*pi/N).^kk).*xk2,xk1-(exp(-1i*pi/N).^kk).*xk2];subplot(2,1,1);stem(k,abs(xk),'r');grid on;xlabel('k');ylabel('2*NFFT');title('n點DFT完成2N點FFT');xn=cos(2*pi/N*7*n2)+1/2*cos(2*pi/N*19*n2);xkk=fft(xn,2*N);subplot(2,1,2);stem(k,abs(xkk),'g');grid on;xlabel('k');ylabel('2*NFFT');title('Matlab2N點FFT');

七、問題思考

1.兩次編程計算與理論計算相比較，結果一致嗎？說明原因。

2.兩次編程計算在編程方面、運算量上的比較。

八、心得體會

第三篇：數字信號處理實驗講稿

邯 鄲 學 院

講 稿

2010 ～2011 學年 第 一 學期

分院(系、部)： 信息工程學院 教 研 室： 電子信息工程 課 程 名 稱： 數字信號處理

授 課 班 級： 07級電子信息工程

主 講 教 師： 王苗苗 職

稱：

助教（研究生）

使 用 教 材： 《數字信號處理》

制 作 系 統：

Word2003

邯鄲學院制

實驗一..Matlab仿真軟件介紹

一、實驗目的

熟悉Matlab仿真軟件

二、實驗設備和元器件

含Matlab仿真軟件的計算機

三、實驗內容和步驟

1、學習Matlab仿真軟件的安裝

2、熟悉Matlab仿真軟件的操作環境

3、直接在Matlab仿真軟件的命令窗口實現數值計算

4、編寫M文件

四、實驗報告要求

按照《Matlab程序設計》模板提交實驗報告

五、預習要求

1、熟悉Matlab仿真軟件

2、參閱Matlab及在電子信息類課程中的應用（第2版）唐向宏 電子工業出版社

實驗二 離散信號和系統分析的Matlab實現

一、實驗目的

1、Matlab實現離散信號和系統分析

2、進一步熟悉Matlab軟件操作

二、實驗設備和元器件

含Matlab仿真軟件的計算機

三、實驗內容和步驟

1、利用Matlab產生離散信號

2、利用Matlab計算離散卷積

3、利用Matlab求解離散LTI系統響應

4、利用Matlab計算DTFT

5、利用Matlab實現部分分式法

6、利用Matlab計算系統的零極點

7、利用Matlab進行簡單數字濾波器設計

四、實驗報告要求

按照《Matlab程序設計》模板提交實驗報告

五、預習要求

預習課本上的相關內容

實驗三 利用Matlab實現信號DFT的計算

一、實驗目的

1、Matlab實現信號DFT的計算

2、進一步熟悉Matlab軟件操作

二、實驗設備和元器件

含Matlab仿真軟件的計算機

三、實驗內容和步驟

1、利用Matlab計算信號的DFT

2、利用Matlab實現由DFT計算線性卷積

四、實驗報告要求

按照《Matlab程序設計》模板提交實驗報告

五、預習要求

預習課本上的相關內容

實驗四 利用Matlab實現濾波器設計

一、實驗目的

1、Matlab實現實現濾波器設計

2、進一步熟悉Matlab軟件操作

二、實驗設備和元器件

含Matlab仿真軟件的計算機

三、實驗內容和步驟

1、利用Matlab實現模擬低通濾波器的設計

2、利用Matlab實現模擬域頻率變換

3、利用Matlab實現脈沖響應不變法

4、利用Matlab實現雙線性變換法

5、利用Matlab實現數字濾波器設計

四、實驗報告要求

按照《Matlab程序設計》模板提交實驗報告

五、預習要求

預習課本上的相關內容

實驗五 利用Matlab實現FIR濾波器設計

一、實驗目的

1、Matlab實現實現濾波器設計

2、進一步熟悉Matlab軟件操作

二、實驗設備和元器件

含Matlab仿真軟件的計算機

三、實驗內容和步驟

1、利用Matlab實現窗函數法

2、利用Matlab實現頻率取樣法

3、利用Matlab實現優化設計

四、實驗報告要求

按照《Matlab程序設計》模板提交實驗報告

五、預習要求

預習課本上的相關內容

實驗六..隨機信號功率譜估計的Matlab實現

一、實驗目的

1、Matlab實現實現濾波器設計

2、進一步熟悉Matlab軟件操作

二、實驗設備和元器件

含Matlab仿真軟件的計算機

三、實驗內容和步驟

1、利用Matlab實現隨機序列

2、利用Matlab計算相關函數的估計

3、利用Matlab進行非參數功率譜估計

4、利用Matlab進行AR模型功率譜估計

四、實驗報告要求

按照《Matlab程序設計》模板提交實驗報告

五、預習要求

預習課本上的相關內容

實驗七..數字濾波器結構的Matlab實現

一、實驗目的

1、Matlab實現實現濾波器設計

2、進一步熟悉Matlab軟件操作

二、實驗設備和元器件

含Matlab仿真軟件的計算機

三、實驗內容和步驟

1、利用Matlab實現數字濾波器直接型設計

2、利用Matlab實現數字濾波器級聯設計

3、利用Matlab實現數字濾波器并聯型設計

4、利用Matlab實現數字濾波器格型設計

四、實驗報告要求

按照《Matlab程序設計》模板提交實驗報告

五、預習要求

預習課本上的相關內容

實驗八....利用Matlab實現信號小波分析

一、實驗目的

1、Matlab實現實現濾波器設計

2、進一步熟悉Matlab軟件操作

二、實驗設備和元器件

含Matlab仿真軟件的計算機

三、實驗內容和步驟

1、小波測試信號

2、分解與重構濾波器組

3、離散小波變換

4、離散小波反變換

5、基于小波的信號去噪

6、基于小波的信號壓縮

四、實驗報告要求

按照《Matlab程序設計》模板提交實驗報告

五、預習要求

預習課本上的相關內容

第四篇：數字信號處理實驗

實驗一 自適應濾波器

一、實驗目的

1、掌握功率譜估計方法

2、會用matlab對功率譜進行仿真

二、實驗原理

功率譜估計方法有很多種，一般分成兩大類，一類是經典譜估計；另一類是現代譜估計。經典譜估計可以分成兩種，一種是BT法，另一種是周期法；BT法是先估計自相關函數，然后將相關函數進行傅里葉變換得到功率譜函數。相應公式如下所示：

1?xx(m)?rN??PBTm???N?|m|?1?n?0x*(n)x(n?m)(1?1)??

(1?2)?xx(m)e?jwnr周期圖法是采用功率譜的另一種定義，但與BT法是等價的，相應的功率譜估計如下所示：

1jw?

Pxx(e)?N其計算框圖如下所示：

觀測數據x(n)FFT?jwnx(n)e?n?0N?120?n?N?1(1?3)

取模的平方1/N Pxx(ejw)?

圖1.1周期圖法計算用功率譜框圖 由于觀測數據有限，所以周期圖法估計分辨率低，估計誤差大。針對經典譜估計的缺點，一般有三種改進方法：平均周期圖法、窗函數法和修正的周期圖平均法。

三、實驗要求

信號是正弦波加正態零均值白噪聲，信噪比為10dB，信號頻率為2kHZ，取樣頻率為100kHZ。

四、實驗程序與實驗結果（1）用周期圖法進行譜估計 A、實驗程序： %用周期法進行譜估計 clear all;N1=128;%數據長度 N2=256;N3=512;N4=1024;f=2;%正弦波頻率，單位為kHZ fs=100;%抽樣頻率,單位為kHZ n1=0:N1-1;n2=0:N2-1;n3=0:N3-1;n4=0:N4-1;a=sqrt(20);%由信噪比為10dB計算正弦信號的幅度 wn1=randn(1,N1);xn1=a*sin(2*pi*f*n1./fs)+wn1;Pxx1=10*log10(abs(fft(xn1).^2)/N1);%周期法求功率譜 f1=((0:length(Pxx1)-1))/length(Pxx1);wn2=randn(1,N2);xn2=a*sin(2*pi*f*n2./fs)+wn2;Pxx2=10*log10(abs(fft(xn2).^2)/N2);f2=((0:length(Pxx2)-1))/length(Pxx2);wn3=randn(1,N3);xn3=a*sin(2*pi*f*n3./fs)+wn3;Pxx3=10*log10(abs(fft(xn3).^2)/N3);f3=((0:length(Pxx3)-1))/length(Pxx3);wn4=randn(1,N4);xn4=a*sin(2*pi*f*n4./fs)+wn4;Pxx4=10*log10(abs(fft(xn4).^2)/N4);f4=((0:length(Pxx4)-1))/length(Pxx4);subplot(2,2,1);plot(f1,Pxx1);xlabel('頻率');ylabel('功率(dB)');title('功率譜Pxx，N=128');subplot(2,2,2);plot(f2,Pxx2);xlabel('頻率');ylabel('功率(dB)');title('功率譜Pxx，N=256');subplot(2,2,3);plot(f3,Pxx3);xlabel('頻率');ylabel('功率(dB)');title('功率譜Pxx，N=512');subplot(2,2,4);plot(f4,Pxx4);xlabel('頻率');ylabel('功率(dB)');title('功率譜Pxx，N=1024');B、實驗仿真結果：

（2）采用漢明窗，分段長度L=32，用修正的周期圖求平均法進行譜估計 A：實驗程序: clear all;N=512;%數據長度 Ns=32;%分段長度

f1=2;%正弦波頻率，單位為kHZ fs=100;%抽樣頻率,單位為kHZ n=0:N-1;a=sqrt(20);%由信噪比為10dB計算正弦信號的幅度 wn=randn(1,N);xn=a*sin(2*pi*f1*n./fs)+wn;w=hamming(32)';%漢明窗

Pxx1=abs(fft(w.*xn(1:32),Ns).^2)/norm(w)^2;Pxx2=abs(fft(w.*xn(33:64),Ns).^2)/norm(w)^2;Pxx3=abs(fft(w.*xn(65:96),Ns).^2)/norm(w)^2;Pxx4=abs(fft(w.*xn(97:128),Ns).^2)/norm(w)^2;Pxx5=abs(fft(w.*xn(129:160),Ns).^2)/norm(w)^2;Pxx6=abs(fft(w.*xn(161:192),Ns).^2)/norm(w)^2;Pxx7=abs(fft(w.*xn(193:224),Ns).^2)/norm(w)^2;Pxx8=abs(fft(w.*xn(225:256),Ns).^2)/norm(w)^2;Pxx9=abs(fft(w.*xn(257:288),Ns).^2)/norm(w)^2;Pxx10=abs(fft(w.*xn(289:320),Ns).^2)/norm(w)^2;Pxx11=abs(fft(w.*xn(321:352),Ns).^2)/norm(w)^2;Pxx12=abs(fft(w.*xn(353:384),Ns).^2)/norm(w)^2;Pxx13=abs(fft(w.*xn(385:416),Ns).^2)/norm(w)^2;Pxx14=abs(fft(w.*xn(417:448),Ns).^2)/norm(w)^2;Pxx15=abs(fft(w.*xn(449:480),Ns).^2)/norm(w)^2;Pxx16=abs(fft(w.*xn(481:512),Ns).^2)/norm(w)^2;Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7+Pxx8+Pxx9+Pxx10+Pxx11+Pxx12+Pxx13+Pxx14+Pxx15+Pxx16)/16);f=((0:length(Pxx)-1))/length(Pxx);plot(f,Pxx);xlabel('頻率');ylabel('功率(dB)');title('加窗平均周期圖法功率譜Pxx，N=512');grid on;B：實驗仿真結果：

五．參考文獻：

[1] 丁玉美，闊永紅，高新波.數字信號處理-時域離散隨機信號處理[M].西安：西安電子科技大學出版社，2002.[2] 萬建偉，王玲.信號處理仿真技術[M].長沙：國防科技大學出版社，2008.實驗二

卡爾曼濾波器的設計

一．實驗目的

1.熟悉并掌握卡爾曼濾波、自適應濾波和譜估計的原理。

2.可以仿真符合要求的卡爾曼濾波器、自適應濾波器和各種譜估計方法。

3.掌握卡爾曼濾波器的遞推公式和仿真方法。4.熟悉matlab的用法。二．實驗原理

卡爾曼濾波是用狀態空間法描述系統的，由狀態方程和測量方程所組成?？柭鼮V波用前一個狀態的估計值和最近一個觀測數據來估計狀態變量的當前值，并以狀態變量的估計值的形式給出。其狀態方程和量測方程如下所示：

xk?1?Akxk?wkyk?Ckxk?vk(1?1)

(1?2)其中，k表示時間，輸入信號wk是一白噪聲，輸出信號的觀測噪聲vk也是一個白噪聲，輸入信號到狀態變量的支路增益等于1，即B=1；A表示狀態變量之間的增益矩陣，可隨時間變化，Ak表示第k次迭代的取值，C表示狀態變量與輸出信號之間的增益矩陣，可隨時間變化，其信號模型如圖1.1所示(k用k?1代替)。

vkCk+ykWk-1+XkZ-1Xk-1Ak-1

圖1.1 卡爾曼濾波器的信號模型

卡爾曼濾波是采用遞推的算法實現的，其基本思想是先不考慮輸入信號wk和觀測噪聲vk的影響，得到狀態變量和輸出信號的估計值，再用輸出信號的估計誤差加權矯正狀態變量的估計值，使狀態變量估計誤差的均方值最小。其遞推公式如下所示：

?x?k?e?0.02x?k?1?Hk(yk?e?0.02xk?1)??1??H?P(P?1)?kkk?

?Pk??e?0.04Pk?1?1?e?0.04??Pk?(I?Hk)Pk??(1?12a)(1?12b)(1?12c)(1?12d)假設初始條件Ak,Ck,Qk,Rk,yk,xk?1,Pk?1已知，其中x0?E[x0],P0?var[x0]，那么遞推流程見圖1.2所示。

Pk?1式(1-5)?Pk'式(1-4)Hk式(1-3)式(1-6)

圖1.2 卡爾曼濾波遞推流程圖

三．實驗要求

一連續平穩的隨機信號x(t),自相關rx(?)?e

??xkPk?，信號x(t)為加性噪聲所干擾，噪聲是白噪聲，測量值的離散值y(k)為已知。Matlab仿真程序如下：

%編卡爾曼濾波遞推程序，估計信號x(t)的波形 clear all;clc;Ak=exp(-0.02);%各系數由前面確定； Ck=1;Rk=0.1;p(1)=20;%各初值； Qk=1-exp(-0.04);

p1(1)=Ak*p(1)*Ak'+Qk;%由p1代表p'； x(1)=0;%設信號初值為0；

H(1)=p1(1)*Ck'*inv(Ck*p1(1)*Ck'+Rk);

zk=[-3.2,-0.8,-14,-16,-17,-18,-3.3,-2.4,-18,-0.3,-0.4,-0.8,-19,-2.0,-1.2,-11,-14,-0.9,0.8,10,0.2,0.5,-0.5,2.4,-0.5,0.5,-13,0.5,10,-12,0.5,-0.6,-15,-0.7,15,0.5,-0.7,-2.0,-19,-17,-11,-14] %zk為測量出來的離散值； N=length(zk);%要測量的點數； for k=2:N

p1(k)=Ak*p(k-1)*Ak'+Qk;%未考慮噪聲時的均方誤差陣； H(k)=p1(k)*Ck'*inv(Ck*p1(k)*Ck'+Rk);%增益方程； I=eye(size(H(k)));%產生和H（k）維數相同的單位矩陣； p(k)=(I-H(k)*Ck)*p1(k);%濾波的均方誤差陣；

x(k)=Ak*x(k-1)+H(k)*(zk(k)-Ck*Ak*x(k-1));%遞推公式； end, x %顯示信號x(k)的數據； m=1:N;n=m*0.02;

plot(n,zk,'-r*',n,x,'-bo');%便于比較zk和x(k)在同一窗口輸出； xlabel('t/s','Fontsize',16);ylabel('z(t),x(t)','fontsize',16);

title('卡爾曼濾波遞推——x(t)的估計波形與z(t)波形','fontsize',16)legend('觀測數據z(t)','信號估計值x(t)',2);grid;四．實驗結果

五．實驗小結

通過卡爾曼濾波估計信號與觀測信號比較知，卡爾曼濾波輸出的估計信號x(t)與實際觀測到的離散值z(t)還是存在一定的誤差，卡爾曼濾波是從初始狀態 就采用遞推方法進行濾波，那么在初值迭代后的一段時間內可能會出現較大的誤差，隨著迭代進行，各參數逐漸趨于穩定，后面的估計值與觀察值的誤差就減少了。

六．參考文獻

[1] 丁玉美，闊永紅，高新波.數字信號處理-時域離散隨機信號處理[M].西安：西安電子科技大學出版社，2002.[2] 萬建偉，王玲.信號處理仿真技術[M].長沙：國防科技大學出版社，2008.

第五篇：數字信號處理實驗-FFT的實現

實

驗

報

告

學生姓名：

學 號：

指導教師：

一、實驗室名稱：數字信號處理實驗室

二、實驗項目名稱：FFT的實現

三、實驗原理：

一．FFT算法思想：

1．DFT的定義：

對于有限長離散數字信號{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，稱為旋轉因子。

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

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

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

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

3．基2按時間抽?。―IT）的FFT算法思想：

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

將長度為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上式中利用了旋轉因子的可約性，即：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點偶數點和奇數點序列的N/2點DFT值，所以，一個N點序列的DFT可以用兩個N/2點序列的DFT組合而成。但是，從上式可以看出，這樣的組合僅表示出了X[k]前N/2點的DFT值，還需要繼續利用X1[k],X2[k]表示X[k]的后半段本算法推導才完整。利用旋轉因子的周期性，有：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為整數（如果序列長度不滿足此條件，通過在后面補零讓其滿足）。

在把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為偶數時：

N/2?1X[2r]??n?0[x[n]?x[n?N2]]?WN,k?0,1,...,N?1 N22rnN/2?1k為奇數時：X[2r?1]?前面已經推導過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]的偶數點值X[2r]和奇數點值X[2r?1]分別可以由組合而成的N/2點的序列來求得，其中偶數點值X[2r]為輸入x[n]的前半段和后半段之和序列的N/2點DFT值，奇數點值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命令：

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

四、實驗目的：

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

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

五、實驗內容：

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

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

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

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

七、實驗步驟：

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

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

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

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

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

八、實驗數據及結果分析：

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

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

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

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

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

九、實驗結論：

十、總結及心得體會：

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