数字信号处理实验作业完全版.doc

郑州轻工业学院 通信工程 索xx
实验1：
理想采样信号的序列，幅度谱，相位谱，以及改变参数后的图像。
源程序：
clc;
n=0:50;
A=444.128;
a=50*sqrt(2.0)*pi;
T=0.001;
w0=50*sqrt(2.0)*pi;
x=A*exp(-a*n*T).*sin(w0*n*T);
close all
subplot(3,2,1);stem(x,’.’);
title('理想采样信号序列');
k=-25:25;
W=(pi/12.5)*k;
X=x*(exp(-j*pi/12.5)).^(n'*k);
magX=abs(X);
subplot(3,2,2);stem(magX,’.’);title('理想采样信号序列的幅度谱');
angX=angle(X);
subplot(3,2,3);stem(angX);title('理想采样信号序列的相位谱')
n=0:50;
A=1;a=0.4,w0=2.0734;T=1;
x=A*exp(-a*n*T).*sin(w0*n*T);
subplot(3,2,4);stem(x,’.’);
title('理想采样信号序列');
k=-25:25;
W=(pi/12.5)*k;
X=x*(exp(-j*pi/12.5)).^(n'*k);
magX=abs(X);
subplot(3,2,5);
stem(magX,’.’)
title('理想采样信号序列的幅度谱');

上机实验答案：
分析理想采样信号序列的特性
产生在不同采样频率时的理想采样信号序列Xa(n)，并记录各自的幅频特性，观察频谱“混淆”现象是否明显存在，说明原因。

源程序：
A=444.128;a=50*pi*sqrt(2.0);W0=50*pi*sqrt(2.0);n=-50:1:50;T1=1/1000;
Xa=A*(exp(a*n*T1)).*(sin(W0*n*T1));
subplot(3,3,1);plot(n,Xa);
title('Xa序列');xlabel('n');ylabel('Xa');k=-25:25;
X1=Xa*(exp(-j*pi/12.5)).^(n'*k);
subplot(3,3,2); stem(k,abs(X1),'.');
title('Xa的幅度谱');
xlabel('k');ylabel('·幅度');
subplot(3,3,3);stem(k,angle(X1),'.');
title('Xa的相位谱');xlabel('k');ylabel('相位');
T2=1/300;Xb=A*(exp(a*n*T2)).*(sin(W0*n*T2));
subplot(3,3,4);plot(n,Xb);
title('Xb序列');xlabel('n');ylabel('相位');k=-25:25;
X2=Xb*(exp(-j*pi/12.5)).^(n'*k);
subplot(3,3,5); stem(k,abs(X2),'.');
title('Xb的幅度谱');
xlabel('k');ylabel('·幅度');subplot(3,3,6);stem(k,angle(X2),'.');
title(' Xb的相位谱');xlabel('k');ylabel('相位');
T3=1/200;Xc=A*(exp(a*n*T3)).*(sin(W0*n*T3));
subplot(3,3,7);plot(n,Xc);title('Xc序列');
xlabel('n');ylabel('Xc');k=-25:25;
X3=Xc*(exp(-j*pi/12.5)).^(n'*k);
subplot(3,3,8); stem(k,abs(X3),'.');
title('Xc的幅度谱');
xlabel('k');ylabel('幅度');subplot(3,3,9);
stem(k,angle(X3),'.'); title('Xc的相位谱');xlabel('k');ylabel('相位');

由图可以看出：
当采样频率为1000Hz时，采样序列在折叠频率附近处，无明显混叠。
当采样频率进一步降低时，主瓣宽度逐渐变宽，频率混叠现象也逐渐严重，存在较明显的失真现象。
原因：信号理想采样后的频谱是原信号频谱的周期严拓，其严拓周期等于采样角频率。根据时域采样定理，如果原信号是带限信号，且采样频率高于原信号最高频率的2倍才不会发生频率混叠。否则频谱会发生明显混叠如采样频率为200和300Hz时。
1、 离散信号、系统和系统响应的分析
1) 观察信号Xb(n)和Hb(n)的时域和频谱特性，比较系统输出Y(n)与Hb(n)的时域及幅频特性，绘出图形。
clc
n=1:50;
hb=zeros(1,50);
hb(1)=1;hb(2)=2.5;hb(3)=2.5;hb(4)=1;
subplot(3,3,1);stem(hb);title('系统hb[n]');
m=1:50;T=0.001;
A=444.128;a=sqrt(2.0)*pi;
w0=50*sqrt(2.0)*pi
x=A*exp(-a*m*T).*sin(w0*m*T);
subplot(3,3,2);stem(x);title('输入信号x[n]');
y=conv(x,hb);
subplot(3,3,3);stem(y);title('输出信号y[n]');
k=-25:25;
X=x*(exp(-j*pi/12.5)).^(n'*k);
magX=abs(X);
subplot(3,3,4);stem(magX);title('输入信号幅度谱');
angX=angle(X);
subplot(3,3,5);
stem(angX);title('输入信号相位谱');
Hb=hb*(exp(-j*pi/12.5)).^(n'*k);
magHb=abs(Hb);
subplot(3,3,6);stem(magHb);title('系统响应幅度谱');
angHb=angle(Hb);
subplot(3,3,7);stem(angHb);title('系统响应相位谱');
n=1:99;k=1:99;
Y=y*(exp(-j*pi/12.5)).^(n'*k);
magY=abs(Y);
subplot(3,3,8);stem(magY);title('输出信号幅度谱');
angY=angle(Y);
subplot(3,3,9);stem(angY);title('输出信号相位谱');

2，三角波的序列和频谱
源程序：
for i=1:4;
x(i)=i;
end
for i=5:8
x(i)=9-i;
end
close all;
subplot(2,1,1);stem(x,'.');
title('三角波序列');
subplot(2,1,2);stem(abs(fft(x,16)),'.');
title('三角波序列的频谱');

2，反三角波及其频谱
for i=1:4;
x(i)=5-i;
end
for i=5:8
x(i)=i-4;
end
close all;
subplot(2,1,1);stem(x,'.');
title('反三角波序列');
subplot(2,1,2);stem(abs(fft(x,16)),'.');
title('反三角波序列的频谱');

3，单位冲击信号序列的幅度谱和频率谱
源程序：
clc
n=1:50;
x=zeros(1,50);
x(1)=1;
close all;
subplot(3,2,1);stem(x);
title('单位冲击信号序列');
k=-25:25;
X=x*(exp(-j*pi/12.5)).^(n'*k);
magX=abs(X);
subplot(3,2,2);stem(magX);
title('单位冲击信号的幅度谱');
angX=angle(X);
subplot(3,2,3);stem(angX);
title('单位冲击信号的相位谱');

clc
n=1:5;
x=sign(sign(10-n)+1);
close all;
subplot(3,1,1);stem(x);title('单位冲击信号序列');
k=-25:25;X=x*(exp(-j*pi/25)).^(n'*k);
magX=abs(X);
subplot(3,1,2);stem(magX);title('单位冲击信号的幅度谱');
angX=angle(X);
subplot(3,1,3);stem(angX);title('单位冲击信号的相位谱');

三角波与反三角波的比较：
源程序：
for i=1:4;
x(i)=i;
end
for i=5:8
x(i)=9-i;
end
close all;
subplot(4,2,1);stem(x,'.');
title('三角波序列');
subplot(4,2,2);stem(abs(fft(x,16)),'.');
title('三角波序列的频谱');
for i=1:4;
x(i)=5-i;
end
for i=5:8
x(i)=i-4;
end
subplot(4,2,3);stem(x,'.');
title('反三角波序列');
subplot(4,2,4);stem(abs(fft(x,16)),'.');
title('反三角波序列的频谱');

高斯序列及其参数改变后的图形：
源程序：
n=0:15;
p=8;q=2;x=exp(-1*(n-p).^2/q);
close all;
subplot(3,1,1);
stem(abs(fft(x)),'.');
title('高斯序列一')
p=8;q=4;x=exp(-1*(n-p).^2/q);
subplot(3,1,2);
stem(abs(fft(x)),'.');
title('改变参数后的高斯序列二');
p=8;q=8;x=exp(-1*(n-p).^2/q);
subplot(3,1,3);
stem(abs(fft(x)),'.');
title('改变参数后的高斯序列三');

衰减的正弦序列和频谱及改参后的情况
源程序：
n=0:15;
a=0.1;f=0.0625;x=exp(-a*n).*sin(2*pi*f*n);
close all;
subplot(4,1,1);stem(x,'.');
title('衰减正弦序列');
subplot(4,1,2);stem(abs(fft(x)),'.');
title('衰减正弦序列的频谱');
a=0.1;f=0.05625;x=exp(-a*n).*sin(2*pi*f*n);
subplot(4,1,3);stem(x,'.');
title('改参后衰减正弦序列');
subplot(4,1,4);stem(abs(fft(x)),'.');
title('改参后衰减正弦序列的频谱');

三 思考题
对一个有限长序列进行DFT，等价于将该序列周期延拓后进行DFS展开。因为DFS也只是其中一个周期来运算，所以FFT在一定条件下页可以用以分析周期信号序列。如果实正弦信号sin(2pifn)，f=0.1，用16点的FFT来做DFS运算，得到的频谱是信号本身的真实谱吗？
答：不是，因取的点不是一个完整的周期，所以不是真实谱。

实验3：
1， 源程序：
fp=.3;fs=.2;
Wp=fp*2*pi;
Ws=fs*2*pi;
alphap=.8;
alphas=20;
[N,Wn]=cheb1ord(Wp,Ws,alphap,alphas,'s');
[B,A]=cheby1(N,alphap,Wn,'s');
close all;
figure(1);
freqs(B,A);
[BT,AT]=lp2ph(B,A,Wp);
figure(2);
freqs(BT,AT);

2，源程序：
FS=1;
[n,Wn]=buttord(0.2*2*pi,0.3*2*pi,1,25,'s');
[b,a]=butter(n,Wn,'s');
[bz,az]=impinvar(b,a,FS);
freqz(bz,az,512,FS);

3，源程序：
FS=100;T=1/100;
Wp=2/FS*pi*[20000,30000];
Ws=2/FS*pi*[15000,35000];
Rp=1;As=40;
[N,Wc]=buttord(Wp,Ws,Rp,As,'s');
[BB,AB]=butter(N,Wc,'s');freqs(BB,AB);

二 思考题
1双线性变换是一种非线性变换,在实验中你观察到这种非线性关系了吗？应该怎么从那种数字滤波器扶贫特性曲线中可以观察到这种非线性关系？
答：观察到了这种非线性关系。
我们知道双线性变换是靠频率的严重非线性关系得到的，这种Ω与ω的非线性关系使得数字滤波器与模拟滤波器在响应与频率的对应关系上发生畸变。例如，如果模拟滤波器是最平时延滤波器，它的相位与频率是直线关系。但是通过双线性变换所得到的数字滤波器就不再保持相位与频率的直线关系了，因而这个数字滤波器也就不再是一个最平时 延特性的滤波器了！

实验4
源程序：
clc;
clear all;
close all;
wp=0.3*pi;
ws=0.5*pi;
wc=(ws+wp)/2;
tr_width=ws-wp;
M=ceil(6.6*pi/tr_width)+1;
N=[0:1:M-1];
alpha=(M-1)/2;
n=[0:1:(M-1)];
m=n-alpha+eps;
hd=sin(wc*m)./(pi*m);
w_ham=(hanning(M))';
h=hd.*w_ham;
[H,w]=freqz(h,[1],1000,'whole');
H=(H(1:501))';
w=(w(1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
grd=grpdelay(h,[1],w);
delta_w=2*pi/1000;
Rp=-(min(db(1:1:wp/delta_w+1)));
As=-round(max(db(ws/delta_w+1:1:501)));
subplot(2,2,1);
stem(hd);
title('理想冲击响应');
axis([0 M-1 -0.1 0.3]);
ylabel(hd,'.');
subplot(2,2,2);
stem(w_ham，‘’);
title('汉宁窗');
axis([0 M-1 0 1.1]);
ylabel(w_ham,'.');
subplot(2,2,3);
stem(h,'.');
title('实际冲击响应');
axis([0 M-1 -0.1 0.3]);
ylabel('h[n]');
subplot(2,2,4);
plot(w/pi,db);
title('衰减幅度');
axis([0 1 -100 10]);
ylabel('Decibles');