数字信号处理实验.docx
《数字信号处理实验.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验.docx（23页珍藏版）》请在冰豆网上搜索。
数字信号处理实验
数字信号处理实验
电信1班3120503005詹筱珊
1-11给出如下差分方程：
（1）计算并画出单位脉冲响应
。
（2）由此
规定的系统是否稳定？
解：
MATLAB实现程序：
a=[1,-1,0.9]；b=[1]；
x=impseq（0,-20,120）；%输入x为单位脉冲序列
n=[-20:
120]；
h=filter（b,a,x）；%系统输出为单位脉冲响应
stem（n,h,'.'）；axis（[-20,120,-1.1,1.1]）
title（'脉冲响应'）；text（125,-1.1,'n'）；ylabel（'h（n）'）
sum（abs（h））%对单位脉冲响应的模值求和
程序运行结果：
ans=
14.8785
系统是稳定的。
1-12已知某模拟信号
，将它分别用不同的采样频率进行采样得到离散时间信号，试分析在以下两种采样情况下对信号频率的影响。
（1）采样频率
；
（2）采样频率
解：
（1）MATLAB实现程序：
dt=0.00005;t=-0.005:
dt:
0.005;
xa=exp（-1000*abs（t））;
fs=5000,Ts=1/fs;n=-25:
1:
25;
x=exp（-1000*abs（n*Ts））;
N=500;k=0:
1:
N;w=pi*k/N;
X=x*exp（-j*n'*w）;X=real（X）;
W=[-fliplr（w）,w（2:
N+1）];X=[fliplr（X）,X（2:
N+1）];
subplot（2,1,1）；plot（t*1000,xa）；xlabel（'时间：
ms'）;ylabel（'x1（n）'）
title（'离散信号'）；holdon
stem（n*Ts*1000,x）;holdoff
subplot（2,1,2）；plot（w/pi,X）;xlabel（'以pi为单位的频率'）；ylabel（'X1（w）'）;title（'离散信号傅里叶变换'）
运行结果：
（2）MATLAB实现程序：
dt=0.00005;t=-0.005:
dt:
0.005;
xa=exp（-1000*abs（t））;
fs=1000,Ts=1/fs;n=-25:
1:
25;
x=exp（-1000*abs（n*Ts））;
N=500;k=0:
1:
N;w=pi*k/N;
X=x*exp（-j*n'*w）;X=real（X）;
W=[-fliplr（w）,w（2:
N+1）];X=[fliplr（X）,X（2:
N+1）];
subplot（2,1,1）；plot（t*1000,xa）；xlabel（'时间：
ms'）;ylabel（'x1（n）'）
title（'离散信号'）；holdon
stem（n*Ts*1000,x）;holdoff
subplot（2,1,2）；plot（w/pi,X）;xlabel（'以pi为单位的频率'）；ylabel（'X1（w）'）;title（'离散信号傅里叶变换'）
运行结果：
2-25已知某系统的系统函数为
求其零、极点并绘出零、极点图。
解：
MATLAB实现程序：
b=[0.30.10.30.10.2];a=[1-1.21.5-0.80.3];
r1=roots（a）
r2=roots（b）
zplane（b,a）
运行结果：
r1=
0.1976+0.8796i
0.1976-0.8796i
0.4024+0.4552i
0.4024-0.4552i
r2=
0.3236+0.8660i
0.3236-0.8660i
-0.4903+0.7345i
-0.4903-0.7345i
2-26已知因果系统
，绘出
的幅频和相频特性曲线。
解：
MATLAB实现程序：
b=[1,0];a=[1,-0.9]
[H,w]=freqz（b,a,100,'whole'）;magH=abs（H）;phaH=angle（H）;
subplot（2,1,1）,plot（w/pi,magH）;grid
xlabel（''）;ylabel（'幅度'）;title（'幅频响应'）
subplot（2,1,2）;plot（w/pi,phaH/pi）;grid
xlabel（'频率（单位：
pi）'）;ylabel（'相位（单位:
pi）'）;title（'相频响应'）
运行结果：
2-27一个线性时不变系统，描绘它的差分方程为
（1）在
之间求得并画出系统的脉冲响应，从脉冲响应确定系统的稳定性；
（2）画出该系统的幅频、相频特性；
（3）如果此系统的输入为
，在
间求系统的输出
；
（4）讨论当输入改为
时，输出波形如何变化？
为什么？
试根据系统的幅频特性解释。
解：
MATLAB实现程序：
b=[1,2,1];a=[1,-0.5,0.25];
x=impseq（0,0,100）;n=0:
100;
h=filter（b,a,x）;
stem（n,h,'.'）;title（'脉冲响应波形'）;
sum（abs（h））
figure;
[H,w]=freqz（b,a,200）;
magH=abs（H）;phaH=angle（H）;
subplot（2,1,1）,plot（w/pi,magH）;grid;
xlabel（''）;ylabel（'幅度'）；title（'幅频响应'）
subplot（2,1,2）;plot（w/pi,phaH/pi）;grid;
xlabel（'频率（单位：
pi）'）;ylabel（'相位（单位：
pi）'）;title（'相频响应'）;
figure;
n1=0:
200;
x1=5+3*cos（0.2*pi*n1）+4*sin（0.3*pi*n1）;
y1=filter（b,a,x1）;
subplot（2,1,1）;plot（n1,x1）;title（'输入x1波形'）;subplot（2,1,2）;plot（n1,y1）;title（'输出y1波形'）；
figure;
n1=0:
200;
x2=5+3*cos（0.2*pi*n1）+4*sin（0.8*pi*n1）;
y2=filter（b,a,x2）;
subplot（2,1,1）;plot（n1,x2）;title（'输入x2波形'）;subplot（2,1,2）;plot（n1,y2）;title（'输出y2波形'）
运行结果：
ans=
6.5714
3-4
是一4点序列：
（1）计算离散时间傅里叶变换（DTFT）即
并画出它的幅度和相位。
（2）计算
的4点DFT。
解：
（1）MATLAB实现程序：
x=[1,1,1,1];w=[0:
1:
500]*2*pi/500;
[H]=freqz（x,1,w）;
magH=abs（H）;phaH=angle（H）;
subplot（2,1,1）;plot（w/pi,magH）;grid
xlabel（''）;ylabel（'|X|'）;title（'DTFT的幅度'）
subplot（2,1,2）;plot（w/pi,phaH/pi*180）;grid
xlabel（'以pi为单位的频率'）;ylabel（'度'）;title（'DTFT的相角'）
运行结果：
（2）MATLAB实现程序：
N=4;k=0:
N-1;
X=dft（x,N）;
magX=abs（X）,phaX=angle（X）*180/pi
subplot（2,1,1）;plot（w*N/（2*pi）,magH,'--'）;
axis（[-0.1,4.1,0,5]）;holdon
stem（k,magX）;
ylabel（'|X（k）|'）;title（'DFT的幅度：
N=4'）；text（4.3,-1,'k'）
holdoff
subplot（2,1,2）;plot（w*N/（2*pi）,phaH*180/pi,'--'）;
axis（[-0.1,4.1,-200,200]）;holdon
stem（k,phaX）;
ylabel（'度'）；title（'DFT相角：
N=4'）；text（4.3,-200,'k'）
运行结果：
3-12
，利用MATLAB程序求如下
、
。
（1）取
的前10点数据，求N=10点的
、
并作图；
（2）将
补零至100点，求N=100点的
、
并作图。
解：
（1）MATLAB实现程序：
n=[0:
1:
9];x=cos（0.48*pi*n）+cos（0.52*pi*n）;
w=[0:
1:
500]*2*pi/500;
X=x*exp（-j*n'*w）;
magX=abs（X）;
x1=fft（x）;magX1=abs（x1（1:
1:
10））;
k1=0:
1:
9;w1=2*pi/10*k1;
subplot（3,1,1）;stem（n,x）;title（'x（n）,0<=n<=9'）;xlabel（'n'）
axis（[0,10,-2.5,2.5]）;line（[0,10],[0,0]）;
subplot（3,1,2）;plot（w/pi,magX）;title（'DFT幅度'）;xlabel（'频率（单位：
pi）'）;axis（[0,1,0,10]）
subplot（3,1,3）;stem（w1/pi,magX1）;title（'DFT幅度'）；
xlabel（'频率（单位：
pi）'）;axis（[0,1,0,10]）
运行结果：
（2）MATLAB实现程序：
n=[0:
1:
9];y=cos（0.48*pi*n）+cos（0.52*pi*n）;
n1=[0:
1:
99];x=[y（1:
1:
10）zeros（1,90）];
w=[0:
1:
500]*2*pi/500;
X=x*exp（-j*n1'*w）;magx=abs（X）;
x1=fft（x）;magx1=abs（x1（1:
1:
50））;
k1=0:
1:
49;w1=2*pi/100*k1;
subplot（3,1,1）;stem（n,y）;title（'x（n）,0<=n<=90zeros'）;xlabel（'n'）
axis（[0,100,-2.5,2.5]）;line（[0,100],[0,0]）;
subplot（3,1,2）;plot（w/pi,magx）;title（'DTFT幅度'）；xlabel（'频率（单位：
pi）'）;axis（[0,1,0,10]）
subplot（3,1,3）;stem（w1/pi,magx1）;title（'DFT幅度'）;xlabel（'频率（单位：
pi）'）;axis（[0,1,0,10]）
运行结果：
5-7利用脉冲响应不变法设计一个数字巴特沃斯低通滤波器，滤波器的技术要求为：
通带截止频率
rad/s，阻带最小衰减
。
研究不同采样频率对所设计数字滤波器频率响应的影响，设采样频率
分别取1kHz,2kHz,4kHz。
解：
（1）fs=1000Hz
MATLAB实现程序：
wp=200*pi;ws=600*pi;Rp=3;Rs=12;
[n,wn]=buttord（wp,ws,Rp,Rs,'s'）
[b,a]=butter（n,wn,'s'）
[db,mag,pha,w]=freqs_m（b,a,500*2*pi）;
plot（w/（2*pi）,db）;axis（[0,500,-20,1]）;holdon
fs=1000;[bz,az]=impinvar（b,a,fs）;
[db,mag,pha,grd,w]=freqz_m（bz,az）;
plot（0.5*fs*w/pi,db）;axis（[0,500,-20,1]）;holdoff
运行结果：
（2）fs=2000Hz
MATLAB实现程序：
wp=200*pi;ws=600*pi;Rp=3;Rs=12;
[n,wn]=buttord（wp,ws,Rp,Rs,'s'）
[b,a]=butter（n,wn,'s'）
[db,mag,pha,w]=freqs_m（b,a,500*2*pi）;
plot（w/（2*pi）,db）;axis（[0,500,-20,1]）;holdon
fs=2000;[bz,az]=impinvar（b,a,fs）;
[db,mag,pha,grd,w]=freqz_m（bz,az）;
plot（0.5*fs*w/pi,db）;axis（[0,500,-20,1]）;holdoff
运行结果：
（3）fs=4000Hz
MATLAB实现程序：
wp=200*pi;ws=600*pi;Rp=3;Rs=12;
[n,wn]=buttord（wp,ws,Rp,Rs,'s'）
[b,a]=butter（n,wn,'s'）
[db,mag,pha,w]=freqs_m（b,a,500*2*pi）;
plot（w/（2*pi）,db）;axis（[0,500,-20,1]）;holdon
fs=4000;[bz,az]=impinvar（b,a,fs）;
[db,mag,pha,grd,w]=freqz_m（bz,az）;
plot（0.5*fs*w/pi,db）;axis（[0,500,-20,1]）;holdoff
运行结果：
5-9已知
利用双线性变换法设计一个ChebyshevI型数字高通滤波器。
解：
MATLAB实现程序：
Rp=1;Rs=20;T=0.001;fp=300;fs=200;
wp=2*pi*fp*T;ws=2*pi*fs*T;
wp1=（2/T）*tan（wp/2）;ws1=（2/T）*tan（ws/2）;
[n,wn]=cheb1ord（wp1,ws1,Rp,Rs,'s'）;[b,a]=cheby1（n,Rp,wn,'high','s'）;
[bz,az]=bilinear（b,a,1/T）;
[db,mag,pha,grd,w]=freqz_m（bz,az）;plot（w/pi,db）;axis（[0,1,-30,2]）
运行结果：
6-5根据下列技术指标，设计一个FIR低通滤波器。
通带截止频率
rad，通带允许波动
；
阻带截止频率
rad，阻带衰减
。
解：
MATLAB实现程序：
Wp=0.2*pi;Ws=0.3*pi;Rp=0.25;Rs=50;
N=80;n=[0:
1:
N-1];
Wc=（Ws+Wp）/2;
hd=ideal_lp（Wc,N）;
w_han=（hamming（N））';
h=hd.*w_han
[db,mag,pha,grd,w]=freqz_m（h,1）;
subplot（2,2,1）;plot（w/pi,db）;xlabel（'\omega^pi'）;ylabel（'dB'）;
subplot（2,2,2）;plot（w/pi,pha）;xlabel（'\omega^pi'）;ylabel（'pha'）;
subplot（2,2,3）;stem（n,h,'.'）;xlabel（'n'）;ylabel（'h（n）'）;
subplot（2,2,4）;stem（n,w_han,'.'）;xlabel（'n'）;ylabel（'Hammingwindow'）;
运行结果：
6-6利用频率采样法去设计一个线性相位FIR低通滤波器，已知
（1）采样点数
；
（2）采样点数
；设置一个过渡点
。
（3）采样点数
；设置两个过渡点
，
。
解：
（1）MATLAB实现程序：
N=33;alpha=（N-1）/2;k=0:
N-1;wk=（2*pi/N）*k;
Hk=[ones（1,9）,zeros（1,16）,ones（1,8）];
angH=-alpha*（2*pi）/N*k;H=Hk.*exp（i*angH）;
h=real（ifft（H,N））;
[db,mag,pha,grd,w]=freqz_m（h,1）;
[Hr,ww,a,L]=Hr_Type1（h）;
subplot（3,2,1）;
plot（ww/pi,Hr,wk（1:
33）/pi,Hk（1:
33）,'o'）;title（'幅度特性'）;
axis（[0,1,-0.5,1.2]）;
subplot（3,2,2）;plot（w/pi,db）;axis（[0,1,-100,10]）;title（'幅度特性'）;
运行结果：
（2）MATLAB实现程序：
N=33;alpha=（N-1）/2;k=0:
N-1;wk=（2*pi/N）*k;
Hk=[ones（1,9）,0.39,zeros（1,14）,0.39,ones（1,8）];
angH=-alpha*（2*pi）/N*k;H=Hk.*exp（i*angH）;
h=real（ifft（H,N））;
[db,mag,pha,grd,w]=freqz_m（h,1）;
[Hr,ww,a,L]=Hr_Type1（h）;
subplot（3,2,1）;
plot（ww/pi,Hr,wk（1:
33）/pi,Hk（1:
33）,'o'）;title（'幅度特性'）;
axis（[0,1,-0.5,1.2]）;
subplot（3,2,2）;plot（w/pi,db）;axis（[0,1,-100,10]）;title（'幅度特性'）;
运行结果：
（3）MATLAB实现程序：
N=34;alpha=（N-1）/2;k=0:
N-1;wk=（2*pi/N）*k;
Hk=[ones（1,9）,0.5925,0.1099,zeros（1,13）,-0.1099,-0.5925,-ones（1,8）];
angH=-alpha*（2*pi）/N*k;H=Hk.*exp（i*angH）;
h=real（ifft（H,N））;
[db,mag,pha,grd,w]=freqz_m（h,1）;
[Hr,ww,a,L]=Hr_Type2（h）;
subplot（3,2,1）;
plot（ww/pi,Hr,wk（1:
34）/pi,Hk（1:
34）,'o'）;title（'幅度特性'）;
axis（[0,1,-0.5,1.2]）;
subplot（3,2,2）;plot（w/pi,db）;axis（[0,1,-100,10]）;title（'幅度特性'）;
运行结果：
升级会员 