数字信号处理实验，滤波器设计

更新时间：2023-12-30 05:50:01 阅读量：教育文库文档下载

说明：文章内容仅供预览，部分内容可能不全。下载后的文档，内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的，是否完整无缺。

数字信号处理实验指导书（修订版1）

数字信号处理（本科）实验指导书

（修订版1）

戴虹编

2010年6月

数字信号处理实验指导书（修订版1）

实验一信号、系统及系统响应实验二 FFT频谱分析实验三 IIR数字滤波器设计实验四 FIR数字滤波器设计

实验五离散系统的时域及复频域分析

综合实验双音多频（DTMF）通信设计的MATLAB仿真附录实验报告格式

数字信号处理实验指导书（修订版1）

实验一信号、系统及系统响应

一、实验目的

1. 掌握典型序列的产生方法。

2. 掌握DFT的实现方法，利用DFT对信号进行频域分析。

3. 熟悉连续信号经采样前后频谱的变化，加深对时域采样定理的理解。

4. 分别利用卷积和DFT分析信号及系统的时域和频域特性，验证时域卷积定理。二、实验环境

1． Windows2000操作系统 2． MATLAB6.0 三、实验原理 1. 信号采样

对连续信号xa(t)=Aesin(Ω0t)u(t)进行采样，采样周期为T,采样点0≤n<50，得采样序列xa(n)= 。 2. 离散傅里叶变换(DFT) 设序列为x(n),长度为N,则

X(e

2πMjωk

N?1-at

)=DFT[x(n)]=?n?0x(n) e-jωk

jωk

其中ωk=k(k=0,1,2,…,M-1),通常M>N,以便观察频谱的细节。|X(e)|----x(n)的幅频谱。

4.连续信号采样前后频谱的变化

^Xa(jΩ)=

^1T??Xm???a[j(??m?s)]

即采样信号的频谱Xa(jΩ)是原连续信号xa(t)的频谱Xa(jΩ)沿频率轴，以周期重复出现，幅度为原来的倍。 5. 采样定理

由采样信号无失真地恢复原连续信号的条件，即采样定理为：。 6．时域卷积定理

设离散线性时不变系统输入信号为x(n),单位脉冲响应为h(n),则输出信号y(n)= ；由时域卷积定理，在频域中， Y(ejω)=FT[y(n)]= 。四、实验内容 1.分析采样序列特性 (1)程序输入

产生采样序列xa(n)= Ae程序 shiyan11.m

clear % clc %

-anT

sin(Ω0nT)u(n)(0≤n<50),其中A=44.128,a=502π,Ω0=502π

采样频率fs（可变），T=1/fs。(要求写%程序注释)

数字信号处理实验指导书（修订版1）

A=444.128;

a=50*sqrt(2)*pi; % w0=50*sqrt(2)*pi; fs=input('输入采样频率fs='); T=1/fs; N=50; n=0:N-1;

xa=A*exp(-a*n*T).*sin(w0*n*T); % subplot(221);stem(n,xa,'.');grid; % M=100;

[Xa,wk]=DFT(xa,M); % f=wk*fs/(2*pi); % subplot(222);plot(f,abs(Xa));grid; % DFT子函数:DFT.m

function [X,wk]=DFT(x,M)

N=length(x); % n=0:N-1;

for k=0:M-1 wk(k+1)=2*pi/M*k;

X(k+1)=sum(x.*exp(-j*wk(k+1)*n)); % end

(2)实验及结果分析

d. a,b,c中，哪几种情况出现了频谱混叠现象？；出现频谱混叠的原因是。 2.时域离散信号和系统响应分析

(1)hb(n)= δ(n)+2.5δ(n-1)+2.5δ(n-2)+δ(n-3) 程序语句为hb=[1,2.5,2.5,1];

(2)卷积语句：y=conv(x,h)

其中x----输入序列x(n)；h----单位脉冲响应h(n); y-------输出序列y(n)。 3.卷积定理验证

(1)编程实现y(n)=xa(n)*hb(n),其中

xa(n)= Ae-anTsin(Ω0nT)u(n)(0≤n<50), A=1,a=0.4,Ω0=2.0734,T=1 hb(n)=δ(n)+2.5δ(n-1)+2.5δ(n-2)+δ(n-3)

及Y(ejωk)=DFT[y(n)](M=100)，(利用DFT.m)；绘出|Y(ejωk)|波形。 (2)编程实现Xa(e)=DFT[Xa(n)](M=100)及Hb(e)=DFT[hb(n)](M=100)；计算Y(ejωk)= Xa(ejωk) Hb(ejωk)；绘出|Y(ejωk)|波形。

问：(1)和(2)中的|Y(e)| 波形一致吗？；为什么？。

jωkjωk

jωk

数字信号处理实验指导书（修订版1）

4. 正弦序列的产生

设正弦序列x(n)=sin(ωn),取样频率 fs=64Hz,信号频率 f=5Hz, 样点n=0~N-1,(N=64),ω=2πf/fs,编程产生x(n)并绘图。程序： %shiyan141.m clear clc

N=64;fs=64;f=5; % n=0:1:N-1; % w=2*pi*f/fs; % x=sin(w*n); % stem(x,'.') % 实验要求:

1)在%后的空格内填入注释 2)运行上述程序，绘出结果波形。

3)设双音多频信号\为x(n)=sin(ω1n)+sin(ω2n),其中f1=697Hz,f2=1336Hz, 取样频率fs=8000Hz,编程产生x(n)并绘出x(n)的波形。 (n=0~799)

5.拓展题

数字信号处理算法之一 ----时间平均

时间平均可用来消除周期信号中的随机噪声。

对m个周期的振动信号x(n)=s(n)+u(n)[其中s(n)--故障信号，为余弦序列；

u(n)---随机噪声]做时间平均,即把x(n)按周期n分段，将每个周期的对应点相加再做平均，计算时间平均前后的信噪比。 %时间平均程序 program11.m clear clc

n=10; % m=input('m='); % load sip % x=sip;

s=zeros(1,n);

for i=1:m % s=s+x(1+n*(i-1):i*n); end s=s/m;

k=[0:n-1];

subplot(321);plot(k,x(1:n));grid; subplot(322);plot(k,s);grid; i=0:1:n-1;

s0=cos(2*pi*i/n); ps=sum(s0.^2)/n;

pu=1;

snr0=10*log10(ps/pu) %原信噪比

py=sum((s-s0).^2)/n;

snr=10*log10(ps/py) %时间平均后的信噪比

数字信号处理实验指导书（修订版1）

实验三 IIR数字滤波器设计

一、实验目的

(1)掌握用双线性变换法设计IIR数字低通和高通滤波器。 (2)设计低通滤波器对实际心电图信号进行滤波。

(3)设计低通滤波器对含有啸叫噪声的音乐信号进行消噪。

(4)设计IIR数字低通和高通滤波器对某个DTMF(双音多频)信号进行频带分离。二、实验环境

1.Windows98以上操作系统 2.安装MATLAB6.0以上版本三、实验原理

1．选频型数字滤波器的种类有、、和滤波器。 2. 从实现方法上，数字滤波器通常分为和滤波器。

3．IIR滤波器的设计目的是根据技术指标，找到；IIR滤波的MATLAB语句为y= ;

4．双线性变换法设计低通IIR滤波器的步骤是

(1) (2) (3) 四、实验内容

1. 人体心电图信号在测量过程中往往受到工业高频干扰，必须经过低通滤波处理后才能作为判断心脏功能的有用信息。给出一实际心电图信号采样序列样本x(n)，其中存在高频干扰。试以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} 。低通滤波器设计指标：ωp=0.2πrad,ωs=0.3πrad,ap=1dB,as=15dB。

H(z)?0.0007378(1?z)(1?1.268z?1?16?2?1?0.705z?2)(1?1.0106z?1?0.3583z)(1?0.904z?0.215z?2)已设计出H(z) (p300)

3??K?1Hk(z)x(n)H1(z)y1(n)H2(z)H(z)y2(n)H3(z)y3(n)＝y(n)Hk(z)?A(1?2z1?Bkz?1?1?z?2)?Ckz?2,k?1,2,3其中

A=0.09036 ；B1=1.2686,C1=-0.7051 ；B2=1.0106,C2=-0.3583；B3=0.9044,C3=-0.2155 *IIR滤波的Matlab语句：y=filter(b,a,x)

b,a----Hk(z)分子/分母系数；x---输入信号x(n)；y---滤波结果y(n)。 **求频响特性的Matlab语句：[H,w]=freqz(b,a,N)

b,a----Hk(z)分子/分母系数；N---频率点数；H—滤波器的频响特性H(w); w---数字频率轴w。

(1)用IIR低通滤波器(原理设计)对实际的心电信号进行滤波 (IIRECG.m)

数字信号处理实验指导书（修订版1）

%程序：(IIRECG.m)

clear 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];%心电信号x(n)[含高频噪声] A=0.09036; %IIR低通滤波器系数 B1=1.2686;C1=-0.7051; B2=1.0106;C2=-0.3583; B3=0.9044;C3=-0.2155; b=[A 2*A A]; a1=[1 -B1 -C1]; N=128;

[H1,w]=freqz(b,a1,N);%频响特性H1(w)

y1=filter(b,a1,x); %y1(n)是x(n)经H1(z)滤波的结果 a2=[1 -B2 -C2];

[H2,w]=freqz(b,a2,N);%频响特性H2(w)

y2=filter(b,a2,y1); %y2(n)是y1(n)经H2(z)滤波的结果 a3=[1 -B3 -C3];

[H3,w]=freqz(b,a3,N);%频响特性H3(w)

H=H1.*H2.*H3; %总的频响特性H(w)=H1(w)H2(w)H3(w) mag=abs(H); %滤波器的幅频特性

db=20*log10((mag+eps)/max(mag));%幅频特性(dB)

y3=filter(b,a3,y2); %y3(n)是y2(n)经H3(z)滤波的结果 X=fft(x,N); %对x(n)做N点fft,得其频谱X

wx=2*pi*(0:N/2-1)/N; %将坐标轴从频率点k转换为数字频率wx(wx=2*pi*k/N) %%%绘图%%%

subplot(221);plot(x);grid on;title('x(n)');

subplot(222);plot(wx/pi,abs(X(1:N/2)));grid on;title('|X(wx)|'); subplot(223);plot(w/pi,db);grid on;title('|H(w)|(db)'); subplot(224);plot(y3);grid on;title('y3(n)'); 要求：

绘出并分析结果图形。

(2)用现有的Matlab函数设计上述IIR低通滤波器,对实际的心电信号进行滤波，同样绘出题1的结果图形。 Matlab函数:

*根据技术指标ωp ,ωs ,αp ,αs计算巴特沃思滤波器的阶数N和3dB截止频率ωc语句：[N,wc]=buttord(wp/pi,ws/pi,ap,as) 注：数字频率以π为单位。

**根据N,ωc确定数字滤波器H(z)分子/分母多项式的系数b,a 低通 [b,a]=butter(N,wc)

数字信号处理实验指导书（修订版1）

高通 [b,a]=butter(N,wc,’high’) %程序IIRECGbt.m clear 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]; %心电信号x(n)[含高频噪声] wp=0.2*pi;ws=0.3*pi;ap=1;as=15; %低通滤波器技术指标 N=128;

[N1,wc]=buttord(wp/pi,ws/pi,ap,as)% [b,a]=butter(N1,wc) % [H,w]=freqz(b,a,N); % y3=filter(b,a,x); % mag=abs(H); % db=20*log10((mag+eps)/max(mag));%幅频特性(dB)

X=fft(x,N); % wx=2*pi*(0:N/2-1)/N; % %%%绘图%%%

subplot(221);plot(x);grid on;title('x(n)');

subplot(222);plot(wx/pi,abs(X(1:N/2)));grid on;title('|X(wx)|'); subplot(223);plot(w/pi,db);grid on;title('|H(w)|(db)'); subplot(224);plot(y3);grid on;title('y3(n)'); 1)在空格中填写程序注释。 2)运行上述程序，得：

滤波器阶数N= ；3dB截止频率wc= ；

b= ；a= 。 3)修改as=50,得：

N= ；wc= ；

b= ；a= 。

2.文件”yinyue.wav”是含有啸叫噪声的音乐信号，首先对其进行谱分析，观察信号和噪声的频带范围，再设计适当的IIR数字低通滤波器对其进行消噪处理恢复原信号，将结果存入音乐文件”yinyuexiaozao.wav”并用耳机监听消噪前后的音乐信号。 %程序:IIRxiaozaoy.m

clear clc

N1=81920; %做fft的点数

y1=wavread('yinyue.wav'); %读入含噪音乐信号 sound(y1) %播放y1 pause(10) %暂停10s y1=y1*6;

数字信号处理实验指导书（修订版1）

Y1=fft(y1,N1); % PSD1=Y1.*conj(Y1)/N1; % fs=8192;

f=8192/N1*(0:N1/2-1); % %%%绘图%%% figure(1)

subplot(211);plot(real(y1));grid on;

subplot(212);plot(f,PSD1(1:N1/2));axis([0 5000 0 500]);grid on; 实验要求：

(1)运行上述程序，绘出结果波形并监听含噪音乐信号。

(2)从第二张子图可见音乐信号的频带范围在 Hz,啸叫噪声频率为 Hz。 (3)根据%后的要求编写程序，将这些程序连在上述程序之后。 %%IIR低通滤波器%%

fp= ; %通带截止频率fp=1000Hz fs1= ; %阻带截止频率fs1=1200Hz wp=2*pi*fp/fs; %化为数字频率wp和ws

ws= ;

Rp= ;As= ; %Rp=1dB,As=50dB [N,wn]= ;

%根据技术指标计算巴特沃思滤波器的阶数N和3dB截止频率ωc [b,a]= ;

%根据N,ωc确定数字滤波器H(z)分子/分母多项式的系数b,a %%%%%%%%%%%%%%%%

[H,w]=freqz(b,a,fs,'whole'); %频响特性H(w) mag= ; %幅频特性

db= ; %幅频特性(以dB为单位)

y= ; %用所设计的滤波器对y1进行低通滤波消噪,得y %%%%绘图%%%%

figure(2)

subplot(211);plot(w*fs/(2*pi),db);axis([0 fs/2 -400 100]);grid on; subplot(212);plot(real(y));grid on;

sound(real(y)) %播放y

wavwrite(real(y),'yinyuexiaozao.wav');%将结果存入音乐文件：\(4)绘出结果波形并监听已消噪的音乐信号。

3.DTMF(双音多频)信号是按键电话拨号音，其中”2”键是由下列低频音和高频音所构成：x(n)=sin(2πf1/fs×n)+ sin(2πf2/fs×n), f1=697Hz, f2=1336Hz,取样频率fs=8000Hz。

(n=0~1599)

编写程序，完成以下功能：

(1)产生x(n)并绘图，用”sound”语句播放x(n),用耳机监听,并将x(n)存入音乐文件 ”DTMF2.wav”。

(2)设计适当的低通和高通滤波器对x(n)进行滤波，将低频音和高频音分离,得分离后的信号: x1(n)=sin(2πf1/fs×n), x2(n)=sin(2πf2/fs×n), 画出x1(n)和x2(n)，并分别将x1(n)和x2(n)存入音乐文件 ”DTMFdi.wav”和”DTMFgao.wav”，用耳机监听这两种声音。

数字信号处理实验指导书（修订版1）

%程序：DTMF.m

波形：

数字信号处理实验指导书（修订版1）

实验四 FIR数字滤波器设计

一、实验目的

1.掌握用窗函数法设计FIR低通数字滤波器。 2.了解各种窗函数对滤波特性的影响。

3.设计FIR低通滤波器对含噪音乐信号进行消噪。

4.用二维低通FIR滤波器连接图像中文字的断裂部分。二、实验环境

1.Windows98以上操作系统 2.安装MATLAB6.0以上版本三、实验原理

1． FIR滤波器的设计目的是根据技术指标，设计。FIR滤波的MATLAB语句为y= 。 2.窗函数

(1)矩形窗的表达式为： (2)汉宁窗的表达式为： (3)哈明窗的表达式为：四、实验内容

1. 人体心电图信号在测量过程中往往受到工业高频干扰，必须经过低通滤波处理后才能作为判断心脏功能的有用信息。给出一实际心电图信号采样序列样本x(n)，其中存在高频干扰。试以x(n)作为输入序列，用窗函数法设计一个FIR低通滤波器滤除其中的干扰成分。 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} 。低通滤波器设计指标：ωp=0.2πrad,ωs=0.3πrad,ap=1dB,as=50dB。选哈明窗，即w(n)=[0.54-0.46cos(2πn/(N-1))]RN(n)。 *FIR滤波的Matlab语句：y=conv(h,x)

h----滤波器单位脉冲响应h(n)；x---输入信号x(n)；y---滤波结果y(n)。程序：

%用FIR低通滤波器对实际的心电信号进行滤波 (FIRECG.m) clear clc

wp=0.2*pi;ws=0.3*pi;% tr=ws-wp; % N=ceil(8*pi/tr)+1 % n=[0:1:N-1];

wc=(ws+wp)/2; % m=n-(N-1)/2+eps;

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

w_ham=(0.54-0.46*cos(2*pi*n/(N-1)));% h=hd.*w_ham; % [H,w]=freqz(h,[1],1000,'whole'); % mag=abs(H);

db=20*log10((mag+eps)/max(mag)); %

数字信号处理实验指导书（修订版1）

delta_w=2*pi/1000;

ap=-(min(db(1:1:wp/delta_w+1))) %技术指标验证ap,as As=-round(max(db(ws/delta_w+1:1:501))) %%%绘图%%% figure(1)

subplot(221);stem(n,hd,'.');grid on;title('hd(n)');

subplot(222);stem(n,w_ham,'.');grid on;title('哈明窗w_ham(n)'); subplot(223);stem(n,h,'.');grid on;title('h(n)');

subplot(224);plot(w(1:501)/pi,db(1:501));grid on;title('H(w)'); %%%%%%%%%%

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];% 含噪的ECG信号

y=conv(x,h); % X=fft(x,128); % Y=fft(y,128); % wx=2*pi*(0:N/2-1)/N; % %%%绘图%%%

figure(2)

subplot(221);plot(x);axis([0 60 -100 20]);grid on;title('x(n)'); subplot(222);stem(wx/pi,abs(X(1:N/2)),'.');grid on;title('|X(w)|');

subplot(223);plot(y);axis([(N-1)/2 (N-1)/2+60 -100 20]);grid on;title('y(n)'); subplot(224);stem(wx/pi,abs(Y(1:N/2)),'.');grid on;title('|Y(w)|'); %%%绘图%%%

要求：

1) 绘出结果图形，并观察图2的第1和第3张小图，可见，y(n)和x(n)相比，延迟了点。

2) 运行上述程序，得：

滤波器阶数N= ；实际的ap= ；实际的as= ；

可见，在设计相同技术指标的滤波器时，FIR滤波器的阶数 IIR滤波器的阶数。 3) 修改窗函数为汉宁窗，得：

滤波器阶数N= ；实际的ap= ；实际的as= ；

数字信号处理实验指导书（修订版1）

2.文件”noisy.wav”是含有高频噪声的音乐信号，首先对其进行谱分析，观察信号和噪声的频带范围，再设计适当的FIR低通滤波器对其进行消噪处理恢复原信号，将结果存入音乐文件”FIRxiaozao.wav”并用耳机监听消噪前后的音乐信号。 %程序: yinyueFIR.m clear clc

fs=16000;

x=wavread('noisy.wav');%读入声音含噪音乐文件到x Nx=length(x); n=0:Nx-1; Ts=1/fs; t=n*Ts;

L=ceil(log2(Nx)); N=2^L;

X=fft(x,N); %对x做频谱X f=fs/N*(0:N/2-1);

figure(1)

subplot(211);plot(t,x); %画x

subplot(212);plot(f,abs(X(1:N/2)));%画|X(f)|(N/2点）实验要求：

(1)运行上述程序，绘出结果波形并监听含噪音乐信号。

(2)从第二张子图可见音乐信号的频带范围在 Hz,噪声频带在 Hz。 (3)根据%后的要求编写程序，将这些程序连在上述程序之后。 %%%%%%FIR低通滤波器%%%%

fp= ; %通带截止频率fp=1500Hz fs1= ; %阻带截止频率为fs1=1800Hz wp=2*pi*fp/fs; %技术指标wp和ws1 ws1=2*pi*fs1/fs;

tr= ; %过渡带tr=ws1-wp NFIR= ;%滤波器的阶数NFIR n=[0:1:NFIR-1]; %滤波器的时间范围n

wc= ; %理想低通滤波器的截止频率 m= ; %m=n-(NFIR-1)/2+eps; hd= ; %理想低通滤波器hd(n)

w_ham= ;%哈明窗

h= ; %加哈明窗之后的实际低通滤波器h(n) [H,w]=freqz(h,[1],N,'whole'); %低通滤波器的频响特性 mag= ; %幅频特性mag db=20*log10((mag+eps)/max(mag)); %幅频特性(dB)

subplot(221); ;grid on;title('hd(n)');%画hd(n)

subplot(222); ; %画哈明窗w_ham(n) grid on;title('哈明窗w_ham(n)');

subplot(223); ; %画h(n) grid on;title('h(n)');

数字信号处理实验指导书（修订版1）

subplot(224);plot(w(1:N/2)/pi,db(1:N/2)); %画db

grid on;title('H(w)');

%%%%%%%%%%%%%%低通滤波%%%%%%%%%%%%%%%% y= ; %低通滤波消噪y(n)=x(n)*h(n) Y=fft(y,N); ny=0:length(y)-1;

ty=ny*Ts; figure(3)

subplot(211);plot(ty,y); %画y

subplot(212);plot(f,abs(Y(1:N/2)));%画|Y(f)|(N/2点)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% sound(x,fs); %以fs=16000Hz播放x pause(5) %停顿5s

sound(y,fs); %以fs=16000Hz播放y wavwrite(y,'FIRxiaozao.wav');

%将y存入声音文件\

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% (4)绘出结果波形并监听已消除高频噪声的音乐信号。

?11?...4.二维低通滤波器可进行图像平滑。用卷积核h(m,n)?N*N???1.........1??...?1??

N?N对图像f(m,n)进行低通滤波，连接图像中文字的断裂部分。

%程序：wenzi.m

clear clc

f=imread(‘ea.jpg’); %读入文字图像’ea.jpg’到二维数组f f=double(f); %将f转换为双精度型 N=11; %卷积核的尺寸N

h=1/(N*N)*ones(N,N); %卷积核h

y=imfilter(f,h); %用h对f进行二维低通滤波(平滑),得输出图像y H=fft2(h,32,32); %对卷积核h做二维fft subplot(211);imshow(f,[ ]); %显示f subplot(223);mesh(abs(fftshift(H))); %显示H的幅频谱 subplot(224);imshow(y,[ ]); %显示y 实验要求：

(1) 运行上述程序，绘出结果波形。

(2) 改变N=5或N=11,问：输出图像有何改变?

数字信号处理实验指导书（修订版1）

实验五离散系统的时域及复频域分析

一、实验目的

1.掌握Matlab编程求解离散时间系统差分方程的两种方法：迭代法和filter函数法。 2.掌握利用Z变换对系统进行复频域分析。

3.掌握系统零极点的绘制方式及Z反变换的求解方法。 4.熟悉Z变换的应用。二、实验环境

1.Windows2000以上操作系统 2.安装MATLAB6.0以上版本三、实验原理

1. 系统的输出y(n)与输入x(n)及单位脉冲响应h(n)的关系: 2. 一个N阶常系数线性差分方程表达式 :

3.系统函数H(z)与X(z),Y(z)的关系为：H(z)= , N阶常系数线性差分方程对应的系统函数H(z)= 。

4.H(z)的零点定义为：；H(z)的极点定义为：。 5.本实验涉及的Matlab函数 1) 求解差分方程：y=filter(b,a,x)

其中：x---输入信号；y---差分方程的解。b,a----输入/输出信号前的系数。 2) 画零极点图:zplane(b,a,n) 其中： b,a----H(z)分子/分母系数。 3) 求z反变换

[r,p,k]=residuez(b,a)

b---H(z)的分子多项式系数；a---H(z)的分母多项式系数；输出：r,p,k的含义如下：

“residuez”可将H(z)分解为简单有理分式之和： H(z)?r(1)1?p(1)z?1?...?r(n)1?p(n)z?1?k(1)?k(2)z?1?.... ?1?

式<1>中p(1),p(2)…p(n)的集合为列向量p;r(1),r(2)…r(n)的集合为列向量r; k(1),k(2)…的集合为行向量k;

已知r,p,k即可手工写出h(n)的表达式：

h(n)={r(1)[p(1)]n + r(2)[p(2)]n+…r(n)[p(n)]n}u(n)+k(1)δ(n)+k(2)δ(n-1)+…… 四、实验内容

1. 已知：系统的差分方程为：

y(n)-y(n-1)+0.9y(n-2)=x(n)

1)分别利用迭代法和filter函数法求此系统的单位脉冲响应h(n),并绘图。 [1] 迭代法： %程序：shiyan511.m clear

clc N=100;

x=[1 zeros(1,N-1)]; %产生x(n)=δ(n)

数字信号处理实验指导书（修订版1）

1) 查表(频率点矩阵sm)将频率点转换为对应数字键。

数组c中不等于0的下标就是各信号的频率点，KL,KH,查表sm，即可将各DTMF信号还原为相应的数字键。 *用到的函数：

i)find(c)----找出c中≠0的数据的下标。 ii)nnz(c)----找出c中≠0的数据的个数。 2)查找过程。

从sm的第一行开始查，查到,则数字键AN=sm的下标减1,并跳出本级循环。 FFT解码程序： %DTMF信号的检测

A=wavread('D2.wav'); subplot(212); plot(A); N=256;

for s=1:8*n

R=out(200*(s-1)+1:200*s); y=fft(R,N);

c(s,:)=abs(y(1:64));r(s,:)=c(s,:); z=find(c(s,:)<40);

c(s,z)=zeros(size(z)); end

for g=1:n

figure(g+1)

subplot(211);plot(r(8*(g-1)+1,:)) subplot(212);plot(c(8*(g-1)+1,:)) end

%还原为数字键

sm=[31 44;23 40;23 44;23 48;26 40;26 44;26 48;28 40;28 44;28 48]; for i3=1:8*n b=nnz(c(i3,:)); if b==2

q1=find(c(i3,:)); for i4=1:10

if q1==sm(i4,:)

AN(i3)=i4-1;break; end end else

AN(i3)=NaN; end end AN

要求：给上述程序加上注释 %...... 2. DFT算法

数字信号处理实验指导书（修订版1）

基本思想：用FFT算法解码每帧信号共涉及256个频率分量，故每帧信号要算N=256点FFT, 但实际上，组成所有DTMF信号只用到8个频率分量,(且已知),于是，可直接利用DFT定义式进行频谱分析，且每帧信号只算8点DFT,以避开FFT中许多无意义的计算，且同样达到解码的目的。解码过程:

(1)接收DTMF信号，并画图。

(2)对信号作DFT（每帧信号作8点）,画幅频谱,从中找出代表各信号的特征字。

N?1X(k)?DFT[x(n)]?j2?N??x(n)Wn?0knN k?0,1,...N-1 （1）其中WN?eWN?ek?j2πNk称为旋转因子，则?j2 ?fkfs?e?e?jωk?cosωk?jsinωk (2)(Hz)这里，fk?????频率点k对应的实际频率 fs?????取样频率(Hz) ωk?????k对应的数字频率（证：?fk?kF??(2)式成立。kfsN?kN?fkfs,即2? kN?rad)2? fkfs?ωk

当x(n)中包含频率为fk的分量时，则x(n)在fk上DFT的幅值|X(k)|与它在其它频率上DFT的幅值相比，是最大的。于是，以200点为一帧，对x(n)在8个特定频率(设为频率1~8）上作DFT,并画幅频谱图，从中找出幅值最大的两条谱线，在解决频谱泄漏现象之后，这两条谱线的横坐标就是代表各信号的特征字(频率1~8中的某两个)。查下表(sm),将各DTMF信号还原为相应的数字键。

数字键对应特征字组合表

2 3

提示：

4 5 1 4 7 * 6 2 5 8 0 7 3 6 9 # 8 A B C D 1）信号的接收(与FFT法一致)

2）直接用DFT公式检测信号(信号矩阵A共8*n行，每行200点)

数字信号处理实验指导书（修订版1）

N?1DFT公式： X(k1)?N?1?x(m)?(cos? m?0k1k1?jsin? k1)m??x(m)?[cos(m? m?0)?jsin(m? k1)] ?1?其中：N?200,m?0~199,k1?1~8(8个频率点)i)建立? k1的矩阵ww?2*pi/8000*[697,770,...1633];ii)Matlab命令：sum(求和)例:A?[1,1];B?[2,3];则：sum(A.*B)?1?2?1?3?5本课题： ?1?中，令d?[cos(m? k1)?jsin(m? k1)](m?0~199)a2?x(m)???某一帧(第l1帧)的信号(200点)a2?A((l1?1)*200?1:l1*200)?X(k1)?sum(d.*a2) 注：k1?1~8,l1?1~8*n另一种方法：(利用DFT公式的矩阵形式,可节省一个循环))[X(1),...X(8)]?a2*w1 (矩阵乘法，与点乘不同其中： w1?exp(?j*2*?/fs).^(m'*fk) m?0~199;fk?[697,770,...1633]iii)作幅频谱，存入r(8n行，8列)

iv)选取适当的阀值，消除频谱泄漏，将处理后的幅频谱存入c(8n行，8列) 3)画图(r,c)(与FFT法一致) 1) 还原为数字键(改sm)

要求：根据以上原理设计DFT法解码程序。

3．卷积法

当用DFT算法求x(n)的频谱X(k)时，发现这也是一卷积过程，下面说明怎样用卷积法解码。

(1)接收DTMF信号，并绘图。

(2)用卷积法检测DTMF信号，找出代表各信号的特征字。卷积法解码的框图如下图所示：

????hk(n)?WNu(n)???yk(n)?x(n)*hk(n)????X(k) 上图方框是一个离散时间系统，广义来说，它就是一个数字滤波器，其单位脉冲响应为：

hk(n)?WNu(n),长度为N的因果序列x(n)通过该系统后，输出为yk(n)，而x(n)在某一频

?knx(n)??kn?n?N-1率点k上的DFT:

X(k)=yk(N-1) <1>

数字信号处理实验指导书（修订版1）

N?1证： yk(n)?x(n)*hk(n)?N?200(一帧信号)N?1?x(m)Wm?0?k(N?1?m)N2?N?k(n?m)N(卷积公式)N?1则：yk(N?1)?其中： W-k(N-1)N?x(m)Wm?0N??1??x(m)Wm?0?k(N?1)N?WNkm????WN?1?kNN?e?j?kN?e?j2?k?cos2? k-jsin2? k?1?yk(N?1)??证明成立。?x(m)Wm?0kmN?X(k)注：(1)由于Matlab中，下标从1开始，所以(2)hk(m)?WN其中：? k1???kmy(N)?yk(N?1)u(m)?cos(? k1m)?jsin(? k1m) m?0~N-1fs?2? k1N-kn

2? fk1 X(k)就是x(n)与单位脉冲响应为hk(n)=WNu(n)的滤波器（系统）的卷积在n=N-1上的输出值。而对不同的k, hk(n)就不同。于是k取DTMF信号用到的8个频率点，且每次卷积x(n)和hk(n)均取N=200点，余下的解码过程同DFT算法。提示：卷积在MATLAB中用“conv”函数实现。要求：根据以上原理设计卷积法解码程序。 4.迭代法

为使解码过程更接近于硬件的实现，可在卷积法的基础上，找到hk(n)满足的差分方程，这样，以接收到的DTMF信号x(n)为输入，用迭代法解此差分方程，则在n=N-1时刻的输出yk(N-1)就是X(k),再画幅频谱|X(k)|或功率谱|X(k)|2,即可对信号进行频谱分析，从而解码。迭代法解码过程如下：

(1) 接收DTMF信号，并绘图。

(2) 用迭代法对DTMF信号进行检测，找出代表各信号的特征字。理论：

[一]求hk(n)满足的差分方程。

1)由hk(n)????Hk(z)???hk(n)的差分方程???knNZ变换由卷积法得?11?W?kNhk(n)?W?Yk(z)X(z)?knNu(n),则其z变换为：Hk(z)??k?1?Wn?0z?n??(WNzn?0?k?1?k?1)nz?1 ，?（1?WNz）Yk(z)?X(z), 即Yk(z)?WNzYk(z)?X(z),hk(n)满足的差分方程为：yk(n)?WNyk(n?1)?x(n)?k对上式两边进行z反变换，得

数字信号处理实验指导书（修订版1）

1)由hk(n)????Hk(z)???hk(n)的差分方程由卷积法得?11?W?kNZ变换(迭代法1)???knNhk(n)?W?Yk(z)X(z)?knNu(n),则其z变换为：Hk(z)??k?1?Wn?0z?n??(WNzn?0?k?1?k?1)nz?1 ，?（1?WNz）Yk(z)?X(z), 即Yk(z)?WNzYk(z)?X(z),hk(n)满足的差分方程为：yk(n)?WNyk(n?1)?x(n)?k对上式两边进行?kz反变换，得即： yk(n)?WNyk(n?1)?x(n)注：可用压缩空间的迭代法求解上式。(2个空间存放yk即可)即： yk(n)=WN-kyk(n-1)+x(n) <1>

则yk(N-1)=X(k)。[其实，yk(N-1)=y(2),即在每个频率上迭代的最后一个值]

[二]Hk(z)分母有理化。

<1>式中，由于WN-k是复数，因此在计算yk(n)时，会碰到复数运算，不方便，为避免复数运算，首先可对Hk(z)进行分母有理化，其基本原理为：用一对复共轭极点代替单极点滤波器的思想。则：Hk(z)?? W?kN11?W?ej2πN?kNkz?1?2πN1?WNz(1?Wk?kNk?1kN?1z?1)(1?Wz)?1?WNz?kkk?1?11?(WN?WN)z?z?2?WkN?e?1-j?2cos(2πN k)?2cosωk ,?Hk(z)?1?WNz1?2cosωk zk-1?z?2 ?2?

Hk(z)分母有理化前后在z平面上的零极点分布如下图所示：

图 1 Hk(z)零极点分布图

从图中可得：Hk(z)原来只有一个极点d1=WN（在单位圆上），经分母有理化后，Hk(z)的极Hk(z)?1?WNzk?1-1-k

1?2cos ωk z?z ??2Yk(z)X(z), ?（1?2cos ωk z?z）Yk(z)?（1?WNz）X(z)-1?2k?1对上式两边进行 z 反变换，得yk(n)?2cos ωkyk(n?1)?yk(n?2)?x(n)?WNx(n?1)即：yk(n)?x(n)?WNx(n?1)?2cos ωkyk(n?1)?yk(n?2) ?3?点变为2个，即d1=WN-k、 d2=WNk,且又出现一个零点c2=WNk,因为d2和c2相互抵消，所以

kkHk(z)的表达式不变，且Hk(z)分母中没有出现复数。由<2>式得

数字信号处理实验指导书（修订版1）

Hk(z)?1?WNzk?1-11?2cosωk z?zk ??2Yk(z)X(z), ?（1?2cosωk z?z）Yk(z)?（1?WNz）X(z)k-1?2k?1对上式两边进行z反变换，得yk(n)?2cosωkyk(n?1)?yk(n?2)?x(n)?WNx(n?1)此时，

即：yk(n)?x(n)?WNx(n?1)?2cosωkyk(n?1)?yk(n?2) ?4?hk(n)系统结构图如下：

图2 hk(n)系统结构图2

式<4>为二阶差分方程，而图中却有三个单位延时环节，所以上图不是标准形式的二阶系统结构图，可通过正则化，将它变换为只有两个延时环节的标准形式(称为正准型)。正则化步骤如下： ?2?式中，令H1(z)?H2(z)?1?WNzk?111?2cosωkz?1?z?2?V(z)X(z) ?2'??Yk(z)V(z) ?2''?,则 Hk(z)?H1(z)?H2(z) v(n)?x(n)?2cosωkv(n?1)?v(n?2)　　?4'?则H1(z)满足的差分方程为H2(z)满足的差分方程为由卷积定理：频域相乘yk(n)?v(n)?WNv(n?1)　　?4''?对应于时域相卷得：hk(n)?h1(n)*h2(n) ?5?相同的单位延时环节，k因此可将hk(n)拆成h1(n)和h2(n)相串联的形式，再合并得下图：

图3 hk(n)系统结构图(正准型)

显然，上图是正准型的二阶系统结构图。综上所述，hk(n)满足的差分方程为 (4’)和（4’’）式，即：

数字信号处理实验指导书（修订版1）

(4')??v(n)?x(n)?2cosωkv(n?1)?v(n?2)　　??ky(n)?v(n)?Wv(n?1)　　 (4'')N?k?且X(k)?yk(N)?v(N)?WNv(N?1) ?6?k

画|X(k)|2,对信号进行频谱分析。

(4’)式中，n=0,1…N-1,共进行N次实数乘法，而(4’’)中，仅在n=N时刻计算1次乘法，但

WNK还是复数。由于解码算法中，不关心X(k),只要求出|X(k)|或|X(k)|2即可进行频谱分析，这里，算|X(k)|。由式<6>得，

|X(k)|?|v(N)?WNv(N?1) | (其中WN?cos?k?jsin?k)?|v(N)?cos?kv(N?1)?jsin?kv(N?1)|2222222k2k2

?[v(N)?cos?kv(N?1)]?sin?kv(N?1)?v(N)?2cos?kv(N)v(N?1)?v(N?1) ?7?

这样解码过程中就完全避免了复数运算。余下解码过程同DFT算法相应部分。要求：根据以上原理设计迭代法解码程序。

数字信号处理实验指导书（修订版1）

附录：实验报告格式