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

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

（修订版1）

戴虹编

2010年6月

目录

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

实验二FFT频谱分析

实验三IIR数字滤波器设计

实验四FIR数字滤波器设计

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

综合实验双音多频（DTMF）通信设计的MATLAB仿真

附录实验报告格式

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

一、实验目的

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

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

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

分别利用卷积和DFT分析信号及系统的时域和频域特性，验证时域卷积定理。

二、实验环境

I.Windows2000操作系统

2.MATLAB6.0

三、实验原理

1.信号采样

2.对连续信号xa⑴=Ac-atsin(COt)u⑴进行采样，采样周期为T,采样点0Wn<50,得采样序

歹I」xa(n)=o

3.离散傅里叶变换(DFT)

设序列为x(n),长度为N,则

N-1

X©3k尸DFT[x(n)]=Zx(n)e-ju，kn,

〃=0

271

其中3k=——k(k=0,l,2，…,M-1),通常M>N,以便观察频谱的细节。|X@")卜___x(n)的幅频谱。

M

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

Xa(j。彳£X/j(Q一〃Q)1

3.即采样信号的频谱a(j。)是原连续信号xa(t)的频谱Xa(j。)沿频率轴，以周期

重复出现，幅度为原来的倍。

4.采样定理

由采样信号无失真地恢复原连续信号的条件，即采样定理为：

6.时域卷积定理

设离散线性时不变系统输入信号为x(n),单位脉冲响应为h(n),则输出信号

y(n)=：由时域卷积定理，在频域中,

Y©*>FT[y(n)];。

四、实验内容

1.分析采样序列特性

(1)程序输入

产生采样序列Xa(n)=Ae-anTsin(QonT)u(n)(OWnv5O),其中A=44.128,a=5072K,Qo=50y/lit

采样频率fs(可变)，T=l/fs。(要求写％程序注释)

程序shiyan1l.m

clear%.

clc%.

A=444.128;

a=50*sqrt(2)*pi;%

w0=50*sqrt(2)*pi;

fs=input('输入采样频率fs=");

T=l/fs;

N=50;

n=O:N-l;

xa=A*exp(-a*n*T).*sin(wO*n*T);%.

subplot(221);stem(n,xa,'.');grid:%

M=JOO;

|XaAvk]=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-l;

fork=0:M-l

wk(k+l)=2*pi/M*k;

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

end

(2)实验及结果分析

a.取fs=1000(Hz),绘出x式n)及|Xa(丁的波形。

b.取fs=300(Hz),绘出xa(n；及|Xa(ej3k)|的波形。

c.取fs=200(Hz),绘出xa(n)及|Xa(ej3k)|的波形。

dab,c中，哪几种情况出现了频谱混叠现象？;出现频谱混叠的原因是

2.时域离散信号和系统响应分析

(1)hb(n)=3(n)+2.5S(n-1)+2.56(n-2)+6(n-3)

程序语句为临二[1,2.525』];

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

其中x--输入序列x(n)；h--单位脉冲响应h(n);y.....瑜出序列y(n)。

3.卷积定理验证

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

anT

xa(n)=Aesin(QonT)u(n)(O^n<5O),A=l,a=0.4,Qo=2.O734,T=l

hb(n)=5(n)+2.5S(n-1)+2.5S(n-2)+8(n-3)

及Y(cj3k尸DFT[y(n)](M=l()0),(利用DFT.m)；绘出|Y(ej3k)|波形。

wk

(2)编程实现Xa©3k)=DFT[Xa(n)](M=100)及HbCe")=DFT[hb(n)](M=100)；

计算Y(dMk)=Xa-)HbM)；绘出|Y(/k)|波形。

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

4.正弦序列的产生

设正弦序列x(n)=sin(3n),取样频率fs=64Hz,信号频率f=5Hz,

样点n=0-N-I,(N=64),w=2nf/fs,编程产生x(n)并绘图。

程序：

%shiyan141.m

clear

clc

N=64;fs=64;f=5;%____________________________________________________________________

n=0:%____________________________________________________________________

w=2*pi*f7fs;%____________________________________________________________________

x=sin(w*n);%_____________________________________________________________________

stem(x,'.')%_____________________________________________________________________

实验要求：

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

3)设双音多频信号“2"为x(n)=sin(3]n)+sin(32n),其中fl=6971Iz,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分段，将每个周期的对•应点相加再做平均，计

算时间平均前后的信噪比，

先时间平均程序programi1.m

clear

clc

n=10;%____________________________________________________________

m=input('m=');%___________________________________________________________

loadsip%____________________________________________________________

x=sip;

s=zcros(1,n);

fori=l:m

s=s+x(1+n*(i-1):i*n);

end

s=s/m;

k=l():n-l];

subplot(32l);plot(k,x(l:n));grid;

subplot(322);plot(k,s);grid;

i=0:l:n-l;

s0=cos(2*pi*i/n);

ps=sum(sO.A2)/n;

pu=l；

snr()=l()*logl()(ps/pu)%原信噪比

py=sum((s-sO).A2)/n;

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

%数据文件sip的产生程序shu.m

clear

clc

n=10;m=100();

i=0:l:n-l;

s=cos(2*pi*i/n);

x=zeros(l.n*m);

u=randn(size(1:n*m));%_______________________________________________

fbrj=l:m%_______________________________________________

x(1+n*(j-1):j*n)=s+u(1+n*(j-l):j*n);

end

savesipx-ascii%_______________________________________________

实验要求：

1)在空格中填入注释。

2)先运行shu.m,产生数据文件sip,再运行programil.m,分别绘出当m=l0,100,500,100001,输

入/输出信号的波形。

实验二FFT频谱分析

一、实验目的

I.理解FFT算法的编程思想。

2.熟练掌握利用FFT对信号作频谱分析。

包括止确地进行参数选择、作频谱图以及读频谱图。

3.了解FFT的应用。

二、实验环境

l.Windows98以上操作系统

2.安装MATLAB6.0以I：版本

三、实验原理

1.谱分析参数选择

I)设信号x⑴最高频率为fc,对其进行取样得x(n),根据取样定理，取样频率fs必须满足:

2)设谱分辨率为F,则最小记录时间tpmi后：取样点数

N二;为使用快速傅里叶变换(FFT)进行谱分析,N还须满足：

N=o

2.用FFT计算信号x(n)的频谱。［设x(n)为实信号］

1)对信号x(n)作N点FFT,得频谱X(k)(k=O~N-1)

X(k尸XR(k)+jXI(k)(k=0~N/2-l),XR(k)—X(k)的实部；XI(k)—X(k)的虚部。

Matlab语句:Y=fft(x,N)

其中:x----x(n)Y----X(k)

2)幅频谱:|X(k)|二，由于x(n)为实信号，因此|X(k)|对称，

Matlab语句:abs(Y)

iii)功率谱:PSD(k)=|X(k)F/N=X(k)X*(k)/N

Matlab语句:PSD=Y.*conj(Y)/N

其中:conj(Y)-X*(k)［X(k)的共趣］

3.读频谱图

频谱图中任意频率点k对应实际频率为:f"o

三、用FFT实现线性卷积运算

用FFT实现y(n)=x(n)*h(n)的步骤为：

1)设x(n)及h(n)的长度分别为N1和N2。为使循环卷积等于线性卷积，用补0的方法使x(n),h(n)

长度均为N,则N须满足N2；为用FFT计算DFT,则N还须满足N二。

2)用FFT计算X(k),H(k);(N点)

3)Y(k)=________________

4)y(n)=“

四、实验内容

理解DIT-FFT算法原理程序，并用它计算X(k)=FFT［R4(n儿分别取N=4,8,16和64,

绘出幅频谱|X(k)|。

%程序DIT.m

clear

cic

x=input('x=');%

N=input('N=');%

x(lcngth(x)+l:N)=zcros(l,N-length(x));%补0x(l:N)

l=log2(N);

xl=zeros(l,N);

forjl=l:N%倒序

x1(j1)=x(bin2dec(fliplr(dec2bin(j1-1,1)))+1);

end

%%FFT(DIT)%%

M=2;

while(M<=N)

W=exp(-2*j*pi/M);%旋转因子W

V=l;

fork=0:l:M/2-l%k为每级蝶形运算旋转因子的个数

fori=0:M:N-1%1为各群的首序号

p=k+i;

q=p+M⑵

A=xl(p+1);

B=xl(q+l)*V;

xl(p+l)=A+B;%本级蝶形运算,xl最终存放X(k)

xl(q+l)=A-B;

end

V=V*W;%旋转因子W的变化

end

M=2*M;%第乂级

end

%%%%%%

subplot(21l);stem(x,'/);gridon;%__________________

title('x(n),);%__________________

subplot(212);stem(abs(x1),'.');gridon;%__________________

title('|X(k)r);%__________________

2.FFT谱分析

设信号为x(t)=sin(2nflt)+sin(2nf2t)+随机噪声，fl=50Hz,f2=120Hz,以取样频率

fs=IkHz对x⑴进行取样，样本长度tp=O.25s,得x(n),对x(n)作256点FFT得频谱X(k),

画原信号x(n),幅频谱|X(k)|以及功率谱PSD(k),对信号进行谱分析。

%程序pufenxi.m

clear

cic

fs=1000;

t=0:1/fs:0.25;1

N=256;%

fl=50;f2=120;%

s=sin(2*pi*fl*t)+sin(2*pi*f2*t);%

x=s+randn(size(t));3信号+噪声x(n)

Y=fft(x,N);%

PSD=Y.*conj(Y)/N;%

f=fs/N*(0:N/2-l);%

subplot(311);plot(x);%

subplot(312);plot(f,abs(Y(1:N/2)));考

subplot(313);plot(f,PSD(1:N/2));%

1)画出图形窗口显示的图形，并注名每个图形的含义。

2)回答下列问题:

i)观察幅频谱图，可以发现，信号x(n)含有的两个频率分量分别是

Hz和Hzo

ii)在该程序中的“f=fs/N%0:N/2-l);”下添加“k=0:N/2-1;”，

aa

plot(f,abs(Y(l:N/2)));"改为plot(kzabs(Y(l:N/2)));”

重新运行该程序并观察幅频谱图，图中两峰值对应的下标分别是

和。它们的含义为0

再将该程序中的N改为512,重新运行该程序并观察幅频谱图，这时图中

两峰值对应的下标分别是和。结果是否和上面的相同？

为什么？。

iii)本例的频谱分辨率F是Hz,改变E2-60HZ,问：在幅频谱中，能否分辨“和

f2

对应的频率分量？。为什么？

再改变f2=52Hz,问：在嗝频谱中，能否分辨fl和f2对应的频率分量？

为什么？。

再改变f2=600Hz,在幅频谱中，f2对应的频率分量出现在Hz；

问：在fs=1000Hz的情况下，能否正确检测f2灼应的频率分量？

3.为什么？。

4.为了正确检测f2对应的频率分量，则fs至少取多少Hz?Hzo在该程序中改变

ts,验证你的结论。

5.iv)比较幅频谱和功率谱，可以发现功率谱具有的特性。

6.FFT实现任意两个序列的快速卷积。

%程序fftjuanji.m

clear

clc

xl=input('xl=*);x2=input("2=')%_______________________

Nl=length(xl);N2=length(x2);%序列x1(n),x2(n)的长度

E=ceil(log2(N1+N2-1));%ceil向+8方向取整

N=2AE;%_______________________

xl=[xl,zeros(lzN-Nl)];%_______________________

x2=[x2,zeros(1,N-N2)];

Xl=fft(xlzN);%_______________________

X2=fft(x2,N);

Y=X1.*X2;%_______________________

y=ifft(Y,N)%

结果分析：

1)回至ljMATLAB窗口，键入：

xl=[l1ILX2=[l2],回车。

结果:y=______________________

2)间：可用Matlab中的什么函数验算上述卷积结果？

4.利用谱分析观察太阳黑子周期性。

以100年中记录到的太阳黑子出现次数为信号x(n),对x(n)作功率谱，从中观察太阳黑子周期

性。

%程序taiyangheizi.m

clear

clC!

x=[101826635317209215412585683823102483

132131118906760474121166471434454348

4228108201512143546413024167428

17365062677148288135712213810386633724

111540629812496666454392L7423559496...

77594447301673774);乳00年中太阳黑子出现的次数

subplot(211);plot(x)彩画x(n)

N-128;fs-1;%fs-lHz,N-128^

s=x-mean(x);告对X作零均值化处理(去除直流分量)

Y=；%对5做N点fft

PSD=;%做功率谱PSD

f-；需将频率定标为实际频率£

subplot(212);;%画功率谱(N/2点)

1)填写空格中的画图语句并绘出结果图形。

2)从s的功率谱观察到，其幅度最高处对应的横坐标f=Hz,

则太阳黑子每隔年出现一次最高峰。

3)在对s做FFTH'J,为何可取ts=lHz,N=128点？

实验三HR数字滤波器设计

一、实验目的

(1)掌握用双线性变换法设计IIR数字低通和高通滤波器.

(2)设计低通滤波器对实际心电图信号进行滤波。

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

(4)设计HR数字低通和高通滤波器对某个DTMFC双音多频)信号进行频带分离。

二、实验环境

l.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}。

低通滤波器设计指标：3P=0.2nrad,9s=0.3nrad,ap=1dB,as=15dBo

=_________________________0.0007378(1+zT,

""(1-1.268Z-,+0.705Z_2)(1-1.0106Z_|+0.3583Z_2)(1-0.904Z-,+0.215Z-2)

已设计出H(z)(p300)

=Kn7Hk(z)

A(l4-2z-'+z-2)

HJz)=，4二1,2,3

其中

A=0.09036；B1=1.2686,C।=-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,aIh(z)分了/分母系数：N频率点数：H滤波器的频响特性Il(\v)；

W----数字频率轴Wo

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

%程序:(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=-O.7051;

B2=1.0106;C2=-0.3583;

B3=0.9044;C3=-0.2155;

b=[A2*AA];

al=[l-Bl-Cl];

N=128;

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

yl=filtcr(b,al,x);%yl(n)是x(n)经Hl(z)滤波的结果

a2=[l-B2-C2];

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

y2=filter(b,a2,yl);%y2(n)是yl(n)经H2(z)滤波的结果

a3=[l-B3-C3];

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

H=H1.*H2.*H3;%总的频响特性H(w)=Hl(w)H2(w)H3(w)

mag=abs(II);%滤波器的幅频特性

db=20*1og10((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-l)/N;先将坐标轴从频率点k转换为数字频率wx(wx=2*pi*k/N)

嬲%绘图嬲

subplot(221);plot(x);gridon;title(,x(n),);

subplot(222);plot(wx/pi,abs(X(l:N/2)));gridon;title。X(wx)|;

subplot(223);plot(w/pi.db);gridon;title(*|H(w)(db)*);

subplot(224);plot(>f3);gridon;title('y3(n)');

要求：

绘出并分析结果图形。

⑵用现有的Matlab函数设计上述1IR低通滤波器,对实际的心电信号进行滤波，同样绘出

题1的结果图形。

Matlab函数：

*根据技术指标3p,cos，与，as计算巴特沃思滤波器的阶数N和3dB截止频率3c语句：

[N,wc]=buttore!(wp/pi,ws/pi,ap,as)

注：数字频率以兀为单位，

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

低通[b,a]=butter(N,wc)

高通[b,a]=butter(N,wc.^igh,)

%程序IIRECGbt.m

clear

clc

x=[-4,-2,0,-4,-6,-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;\vs=0.3*pi;ap=l;as=15;%低通滤波器技术指标

N二128;

[Nl,wc]=buttord(wp/pi,ws/pi,ap,as)%

[b,a]=butter(Nl,wc)%______________________________________________

[H,w]=freqz(b,a,N);%______________________________________________

y3=fi1ter(b,a,x);%_______________________________________________

mag=abs(H);%_______________________________________________

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

X=fft(x,N);%_______________________________________________

wx=2*pi*(0:N/2-l)/N;%______________________________________________

舟溉绘图与嬲

subplot(221);plot(x);gridon;title('x(n)');

subplot(222);plot(wx/pi,abs(X(l:N/2)));gridon;title(，X(wx)|');

subplot(223);plot(w/pi.db);gridon;titleC|H(w)(db),);

subplot(224);plot(y3);gridon;title('y3(n)');

1)在空格中填写程序注释，

2)运行上述程序，得：

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

b=；a=o

3)修改as=50,得:

N=；wc=；

b=；a=o

2.文件"yinyue.wav”是含有啸叫噪声的音乐信号，首先对其进行谱分析，观察信号和噪声的

频带范围，再设计适当的IIR数字低通漉波器对其进行消噪处理恢复原信号，将结果存入音

乐文件"yinyuexiaozao.wav"并用耳机监听消噪前后的音乐信号。

%程序:HRxiaozaoy.m

clear

clc

N1=81920;%做£代的点数

yI=wavread('yinyue.wav');%读入含噪音乐信号

sound(y1)%播放yl

pause(10)%暂停10s

yl=yl*6;

Yl=fft(yl,Nl);%.

PSDl=YI.*conj(Yl)/Nl;%

fs=8192;

f=8l92/NI*(0:Nl/2-l);%

%%%绘图％%%

figure(1)

subplot(21l);plot(real(y1));gridon;

subplot(212);plot(f,PSD1(1:N1/2));axis([050(X)050()J);gridon;

实验要求：

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

(2)从第二张子图可见音乐信号的频带范围在Hz,啸叫噪声频率为Hz。

(3)根据％后的要求编写程序，将这些程序连在上述程序之后。

%%I1R低通滤波器％%

fp=;%通带截止频率fp=1000Hz

fsl=;%阻带截止频率fsl=1200Hz

wp=2*pi*fp/fs;%化为数字频率wp和ws

Rp=;As=;%Rp=ldB,As=50dB

[N,wn]=;

%根据技术指标计算巴特沃思滤波器的阶数N和3dB截止频率3c

[b,a]=；

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

%%%%%%%%%%%%%%%%

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

mag=;%幅频特性

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

y=;%用所设计的滤波器对yl进行低通滤波消噪，得y

%%%%绘图％%%%

figure⑵

subplot(211);plot(w*fs/(2*pi),db);axis([0fs/2-400l(X)]);gridon;

subplot(212);plot(real(y));gridon;

sound(real(y))%播放y

wavwrite(real(y),'yinyuexiaozao.wav');%将结果存入音乐文件:"yinyuexiaozao.wav"

(4)绘出结果波形并监听己消噪的音乐信号。

3.DTMF(双音多频)信号是按键电话拨号音，其中“2”键是由下列低频音和高频音所构成:

x(n)=sin(2nfl/fsXn)+sin(2nf2/fsXn),fl=697Hz,f2=1336Hz,取样频率fs=8(X)0Hzo

(n=0-1599)

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

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

DTMF2.wav”。

(2)设计适当的低通和高通滤波器对x(n)进行滤波，将低频音和高频音分离，得分离后的信号:

xl(n)=sin(2nfl/fsXn),x2(n)=sin(2nf2/fsXn),画出xl(n)和x2(n),并分别将xl(n)和x2(n)存

入音乐文件"DTMFdi.wav"和“DTMFgao.wav"，用耳机监听这两种声音。

%程序:DTMF.m

波形:

实验四FIR数字滤波器设计

一、实验目的

1.掌握用窗函数法设计FIR低通数字滤波器。

2.了解各种窗函数对滤波特性的影响。

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

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

二、实验环境

LWindows98以上操作系统

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.2nrad,s=().3nrad,ap=1dB,as=5()dB。

选哈明窗，BPw(n)=[0.54-0.46cos(2nn/(N-l))]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)+l%_________________________________

n=[0:1:N-1];

wc=(ws+wp)/2;%_________________________________

m=n-(N-l)/2+eps;

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

w_ham=(0.54-0.46*cos(2*pi*n/(N-l)));%

h=hd.*w_ham;%

[H,w]=freqz(h,[1],1000：'who1e');%

mag=abs(H);

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

delta_w=2*pi/I000;

ap=-(min(db(l:1:wp/delta_w+l)))%技术指标验证ap,as

As=-round(max(db(ws/delta_w+l:1:501)))

%%%绘图燃%

figure(l)

subplot(221);stem(n,hd：'.');gridon;title('hd(n)');

subplot(222);stem(n,w_ham,'.');gridon;titie('哈明窗w_ham(n),);

subplot(223);stem(n,hJ.');gridomtitleC,h(n)');

subplot(224);plot(w(l:501)/pi,db(l:501));gridon;title(*II(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-l)/N;%____________________________________________

烬战绘图搬先

figure⑵

subplot(221);plot(x);axis([060-10020]);gridon;x(n),);

subplot(222);stem(wx/pi,abs(X(l:N/2)),*.*);gridon;title(*|X(w)');

subplot(223);plot(y);axis([(N-l)/2(N-l)/2+60-10020]);gridon;title(*y(n)*);

subplot(224);stem(wx/pi,abs(Y(1:N/2))»*.*);gridon;title(*|Y(w)

%%%绘图燃%

1)要求：

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

_________点。

运行上述程序，得：

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

2)可见，在设计相同技术指标的滤波器时，FIR滤波器的阶数HR滤波器的阶

数。

修改窗函数为汉宇窗，得:

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

2.文件"noisy.wav”是含有高频噪声的音乐信号，首先对其进行谱分析，观察信号和噪声的

频带范围，再设计适当的FIR低通滤波器对其进行消噪处理恢复原信号，将结果存入音乐文

件"FIRxiaozao.wavM并用耳机监听消噪前后的音乐信号。

%程序:yinyueFIR.m

clear

cic

fs=16000:

x=wavreadCnoisy.wav，);％读入声音含噪音乐文件到x

Nx=length(x);

n=0:Nx-l;

Ts=l/fs;

t=n*Ts;

L=ccil(log2(Nx));

N=2AL;

X=fft(x,N);%对x做频谱X

f=fs/N*(0:N/2-l);

figured)

subplot(21l);plot(t,x);%画x

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

实验要求：

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

(2)从第一张子图可见音乐信号的频带范围在Hz,噪声频带在HZo

(3)根据％后的要求编写程序，将这些程序连在上述程序之后。

%%%%%%FIR低通滤波器％%%%

fp=;%通带截止频率fp=1500Hz

fsl=;%阻带截止频率为fsl=18(X)Hz

wp=2*pi*fp/fs;%技术指标wp和wsl

wsl=2*pi*fsl/fs;

tr=;%过渡带tr=wsl-wp

NFIR=;%述波器的阶数NFIR

n=[():l:NFIR-l];%淀波器的时间范围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*logl0((mag+eps)/max(mag));%幅频特性(dB)

subplot(221);;gridon;title('hd(n)');%iSjhd(n)

subplot(222);;%画哈明窗w_ham(n)

gridon;titlc(，哈明窗w_ham(n)');

subplot(223);;%画h(n)

gridon:title('h(n)');

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

gridon;title('H(w)');

%%%%%%%%%%%%%%低通滤波％%%%%%%%%%%%%%%%

y=;%低通滤波消噪y(n)=x(n)*h(n)

Y=fft(y,N);

ny=O:length(y)-l;

ty=ny*Ts;

figure(3)

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sound(x,fs);%以fs=16000IIz播放x

pause(5)%停顿5s

sound(y,fs);%以fs=16000Hz播放y

wavwrite(y,'FIRxiaozao.wav');

%将y存入声音文件"FIRxiaozao.wav"

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

4.二维低通滤波器可进行图像平滑。用卷积核出,〃)二..........

J…UNXN

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

%程序:wenzi.m

clear

clc

f=imrcad(；ea.jpg,);%读入文字图像'ea.jpg，到二维数组f

f=double(f);%将£转换为双精度型

N=ll;%卷积核的尺寸N

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

y=imfilter(f,h);%用h对f进行二维低通滤波(平滑)，得输山图像y

H=fft2(h,32,32);%对卷积核h做二维fft

subplot(2il);imshow(f,[]);%显示f

subplot(223);mesh(abs(fftshift(H)));%显示H的幅频谱

subplot(224);imshow(y,[]);%显示y

(1)实验要求：

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

改变N=5或N=ll,问：输出图像有何改变?

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

一、实验目的

1.掌握Matlab编程求解离散时间系统差分方程的两种方法：迭代法和filter函数法。

2.掌握利用Z变换对系统进行更频域分析。

3.掌握系统零极点的绘制方式及Z反变换的求解方法。

4.熟悉Z变换的应用。

二、实验环境

I.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)=o

4.H(z)的零点定义为：；H(z)的极点定义为：。

5.本实验涉及的Matlab函数

1)求解差分方程:产filter(b,a,x)

其中：x—输入信号；y-差分方程的解。ba—输入/输出信号前的系数。

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)分解为简单有理分式之和：

rd)

H(z)=<1>

1-/ADZ-11-P(n)z~l

式＜l＞中p(l),p(2)…p(n)的集合为列向量p;r(l),i•⑵…r(n)的集合为列向量r;

k(l),k(2)…的集合为行向量k;

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

h(n)={r(l)[p(l)]n+r(2)fp(2)ln+...r(n)[p(n)]n}u(n)+k(1)8(n)+k⑵8(n-l)+……

四、实验内容

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

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

1)分别利用迭代法和filler函数法求此系统的单位脉冲响应h(n),并绘图。

[1]迭代法：

%程序：shiyan5ll.m

clear

clc

N=100;

x=[lzeros(l,N-l)];%产生x(n)=5(n)

y(l)=O;%y(-2)=0

y(2)=0;%y(-l)=O

y(3)=i;%y(O)=i

forn=4:N%迭代法求解该差分方程，得