数字信号处理第三版用MATLAB上机实验
实验二：时域采样与频域采样
一、 时域采样
1. 用MATLAB 编程如下:
%1时域采样序列分析fs=1000
A=444.128; a=222.144; w=222.144; ts=64*10^(-3); fs=1000;T=1/fs; n=0:ts/T-1; xn=A*exp((-a)*n/fs).*sin(w*n/fs); Xk=fft(xn); subplot(3,2,1);stem(n,xn);xlabel('n,fs=1000Hz'); ylabel('xn');title('xn');
subplot(3,2,2);plot(n,abs(Xk));xlabel('k,fs=1000Hz'); title('|X(k)|'); %1时域采样序列分析fs=200
A=444.128; a=222.144; w=222.144; ts=64*10^(-3); fs=200;T=1/fs; n=0:ts/T-1; xn=A*exp((-a)*n/fs).*sin(w*n/fs);
Xk=fft(xn);subplot(3,2,3);stem(n,xn);xlabel('n,fs=200Hz'); ylabel('xn');title('xn');
subplot(3,2,4);plot(n,abs(Xk));xlabel('k,fs=200Hz'); title('|X(k)|'); %1时域采样序列分析fs=500
A=444.128; a=222.144; w=222.144; ts=64*10^(-3); fs=500;T=1/fs; n=0:ts/T-1; xn=A*exp((-a)*n/fs).*sin(w*n/fs); Xk=fft(xn);
subplot(3,2,5);stem(n,xn);xlabel('n,fs=500Hz'); ylabel('xn');title('xn');
subplot(3,2,6);plot(n,abs(Xk));xlabel('k,fs=500Hz'); title('|X(k)|');
2. 经调试结果如下图:
xn
1000500
20
4060n,fs=1000Hz
xn
80
00
20
|X(k)|
x
n
4060k,fs=1000Hz
|X
(k)|
80
200100
510n,fs=200Hz
xn
15
00
510k,fs=200Hz|X(k)|
15
x n
500
x n
10
2030n,fs=500Hz
40
00
10
2030k,fs=500Hz
40
实验结果说明:对时域信号采样频率必须大于等于模拟信号频率的两倍以上, 才 能使采样信号的频谱不产生混叠.fs=200Hz时，采样信号的频谱产生了混叠，fs=500Hz和fs=1000Hz时，大于模拟信号频率的两倍以上，采样信号的频谱不产生混叠。
二、 频域采样
1. 用MATLAB 编程如下:
%频域采样定理验证 M=26;N=32;n=0:M; n1=0:13;x1=n1+1; n2=14:26;x2=27-n2;
x=[x1,x2];Xk=fft(x,512); X32k=fft(x,32);
k=0:511;w=(pi/512)*k;
subplot(321);stem(n,x);xlabel('n'); ylabel('xn');axis([0,31,0,15]);
subplot(322);plot(w,abs(Xk));xlabel('k'); ylabel('|X(k)|');axis([0,1,0,200]) X16k=X32k(1:2:N);
x32n=ifft(X32k);x16n=ifft(X16k,16); k1=0:31;k2=0:15;
subplot(323);stem(k1,abs(X32k));xlabel('k'); ylabel('X32k');axis([0,31,0,200]);
subplot(325);stem(k2,abs(X16k));xlabel('k'); ylabel('|X(k)|');axis([0,15,0,200]) n=0:31;
subplot(324);stem(n,abs(x32n));xlabel('n'); ylabel('|x(n)|');axis([0,31,0,15]) n1=0:15;
subplot(326);stem(n1,abs(x16n));xlabel('n'); ylabel('|x(n)|');axis([0,31,0,15])
2. 经调试结果如下图:
200
|X (k ) |
10
n
20
30
x n
1000
0.5k
1
X 32k
|x (n ) |
k
n
|X (k ) |
|x (n ) |
k
10
n
20
30
实验结果说明:对频域采样点数N 必须大于等于时域离散信号的长度M, 才能使时
域不产生混叠.
实验三：用FFT 对信号作频谱分析
一、 x1-x4变换区间N 为8，x5、x6变换区间N 为16
1. 用MATLAB 编程如下:
x1=[1 1 1 1 0 0 0 0];
x2=[1 2 3 4 4 3 2 1];
x3=[4 3 2 1 1 2 3 4]; x4=cos(0.25*pi*n);
N=8;n=0:7;k=0:7; X1k=fft(x1,N);
subplot(2,2,1);stem(n,x1,'.'); xlabel('n');ylabel('|x1(n)|'); subplot(2,2,2);stem(k,abs(X1k),'.'); xlabel('k');ylabel('|X1(k)|'); X2k=fft(x2,N);
subplot(2,2,3);stem(n,x2,'.'); xlabel('n');ylabel('|x2(n)|'); subplot(2,2,4);stem(k,abs(X2k),'.'); xlabel('k');ylabel('|X2(k)|'); figure(2)
X3k=fft(x3,N);
subplot(2,2,1);stem(n,x3,'.'); xlabel('n');ylabel('|x3(n)|');
subplot(2,2,2);stem(k,abs(X3k),'.'); xlabel('k');ylabel('|X3(k)|'); X2k=fft(x4,N);
subplot(2,2,3);stem(n,x4,'.'); xlabel('n');ylabel('|x4(n)|');
subplot(2,2,4);stem(k,abs(X2k),'.'); xlabel('k');ylabel('|X4(k)|'); fs=64;N=16;
n=0:N-1;k=n; x5=cos(n*pi/4)+cos(n*pi/8);
x6=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); X5k=fft(x5,N); X6k=fft(x6,N); figure(3)
subplot(2,2,1);stem(n,x5,'.'); xlabel('n');ylabel('|x5(n)|'); subplot(2,2,2);stem(abs(X5k),'.'); xlabel('k');ylabel('|X5(k)|'); subplot(2,2,3);stem(n,x6,'.'); xlabel('n');ylabel('|x6(n)|'); subplot(2,2,4);stem(abs(X6k),'.'); xlabel('k');ylabel('|X6(k)|');
2. 经调试结果如下图:
|
) (n 1|x n
|
) (n 2|x 0
2
46
8
n
|
) (n 3|x 0
2
46
8
n
n
|
) (k 1X |k
k
k
|
) (k 4X |k
|X 2(k ) |
|X 3(k ) |
|x 4(n ) |
86
|X 5(k )
|
5
n
10
15
|x 5(n ) |
420
k
15
|X 6(k ) |
5
n
10
15
|x 6(n )
|
10
5
05
10k
1520
二、 x1-x4变换区间N 为16，x5、x6变换区间N 为32
1. 用MATLAB 编程如下:
x1=[1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0]; x2=[1 2 3 4 4 3 2 1 0 0 0 0 0 0 0 0]; x3=[4 3 2 1 1 2 3 4 0 0 0 0 0 0 0 0]; x4=cos(0.25*pi*n); N=16;n=0:15;k=0:15; X1k=fft(x1,N);
subplot(2,2,1);stem(n,x1,'.'); xlabel('n');ylabel('|x1(n)|');
subplot(2,2,2);stem(k,abs(X1k),'.'); xlabel('k');ylabel('|X1(k)|'); X2k=fft(x2,N);
subplot(2,2,3);stem(n,x2,'.'); xlabel('n');ylabel('|x2(n)|');
subplot(2,2,4);stem(k,abs(X2k),'.'); xlabel('k');ylabel('|X2(k)|'); figure(2)
X3k=fft(x3,N);
subplot(2,2,1);stem(n,x3,'.'); xlabel('n');ylabel('|x3(n)|');
subplot(2,2,2);stem(k,abs(X3k),'.'); xlabel('k');ylabel('|X3(k)|'); X2k=fft(x4,N);
subplot(2,2,3);stem(n,x4,'.'); xlabel('n');ylabel('|x4(n)|');
subplot(2,2,4);stem(k,abs(X2k),'.'); xlabel('k');ylabel('|X4(k)|'); fs=64;N=32; n=0:N-1;k=n;
x5=cos(n*pi/4)+cos(n*pi/8);
x6=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); X5k=fft(x5,N); X6k=fft(x6,N); figure(3)
subplot(2,2,1);stem(n,x5,'.'); xlabel('n');ylabel('|x5(n)|'); subplot(2,2,2);stem(abs(X5k),'.'); xlabel('k');ylabel('|X5(k)|'); subplot(2,2,3);stem(n,x6,'.'); xlabel('n');ylabel('|x6(n)|'); subplot(2,2,4);stem(abs(X6k),'.'); xlabel('k');ylabel('|X6(k)|');
2. 经调试结果如下图:
|
) (n 1x |n
|
) (n 2x |n
|
) k (1X |k
k
|X 2(k ) |
|X 3(k ) |
n
|x 3(n ) |
k
|X 4(k ) |
n
|x 4(n ) |
k
2015
|X 5(k )
|
10
20n
30
40
|x 5(n ) |
1050
k
2015
|X 6(k ) |
10
20n
30
40
|x 6(n ) |
1050
k
实验结果说明:频谱分辨率与FFT 的变换区间有关，当N 较大时，离散谱的包络才能逼近连续谱，因此N 要选大一些。
实验三：用双线性变换法设计IIR 数字滤波器
一、双线性变换法设计butter worth filter
1. 编写程序如下：
Rp=1;As=15;wp=0.2*pi;ws=0.3*pi;T=1;Fs=1;
Omegap=(2/T)*tan(wp/2);Omegas=(2/T)*tan(ws/2); ksp=((10^(0.1*Rp)-1)/(10^(0.1*As)-1))^0.5 Asp=Omegas/Omegap
N=-(log(ksp)/log(Asp)) N=ceil(N)
wc=Omegas/((10^(As/10)-1)^(1/(2*N))) [z,p,k]=buttap(N) [b,a]=zp2tf(z,p,k) [bt,at]=lp2lp(b,a,wp)
[bb,ab]=bilinear(bt,at,Fs) [h,w]=freqz(bb,ab,501); figure(1);
plot(w/pi,abs(h));t itle('butterworth滤波器的幅频特性') figure(2);
plot(w/pi,angle(h)); title('butterworth滤波器的相频特性')
2. 运行结果如下： ksp =0.0920 Asp = 1.5682 N = 5.3044 N =6
wc = 0.7662 z =[]
p =-0.2588 + 0.9659i -0.2588 - 0.9659i -0.7071 + 0.7071i -0.7071 - 0.7071i -0.9659 + 0.2588i -0.9659 - 0.2588i k =1
b = 0 0 0 0 0 0 1
a = 1.0000 3.8637 7.4641 9.1416 7.4641 3.8637 bt = 0.0615
at=1.0000 2.4276 2.9467 2.2676 1.1633 0.3784 bb=0.0003 0.0017 0.0043 0.0058 0.0043 0.0017 ab=1.0000 -3.6543 5.8693 -5.2192 2.6894 -0.7574
1.0000 0.0615 0.0003 0.0908
3. 运行图形如下：
butterworth 滤波器的幅频特性
butterworth 滤波器的相频特性
二、心电图滤波
1. 编写程序如下：
wp = 0.2*pi; ws = 0.3*pi; Rp = 1; As = 15;
wp1 = (2/T)*tan(wp/2); % Prewarp Prototype Passband freq ws1= (2/T)*tan(ws/2);
[N,wc]=buttord(wp1,ws1,Rp,As,'s') [B,A]=butter(N,wc,'s') [b,a] = bilinear(B,A,T);
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]; k=1; close all; figure(1)
subplot(2,2,1)
n=0:55; stem(n,x,'.'); axis([0 56 -100 50]); hold on;
n=0:60; m=zeros(61);
plot(n,m); xlabel('n'); ylabel('x(n)'); title('心电图信号采样序列x(n)'); y=filter(b,a,x); subplot(2,2,3)
n=0:55; stem(n,y,'.'); axis([0 56 -15 5]);
hold on; n=0:60; m=zeros(61);
plot(n,m); xlabel('n'); ylabel('y(n)'); title('三级滤波后的心电图信号'); axis([0 56 -100 50]); [H,w]=freqz(b,a,100); mag=abs(H);
db=20*log10((mag+eps)/max(mag)); subplot(2,2,2) plot(w/pi,db);
axis([0,0.5,-50,10]); title('滤波器的幅频响应');
2. 运行结果如下：
N = 6
wc =0.7662
B = 0 0 0 0 0 0 0.2024 A = 1.0000 2.9605 4.3822 4.1124 2.5728 1.0205 0.2024
3. 运行图形如下：
心电图信号采样序列x(n)0
滤波器的幅频响应
x (n )
-20
-40
20
40
0.1
0.2
0.3
0.4
n
三级滤波后的心电图信号
y (n )
020
n
40
实验四： 用窗函数法设计FIR 数字滤波器
一、hanning window 用ideal 求hd,h=hd.*(B)' boxcar hamming blackman
1. 编写程序如下： %实验4子函数：产生hd
function hd=ideal(w,N); alpha=(N-1)/2; n=[0:(N-1)]; m=n-alpha+eps; hd=sin(w*m)./(pi*m); wc=0.25*pi;N=33; n=0:N-1; B=boxcar(N) hd=ideal(wc,N);%得到理想低通滤波器 h=hd.*(B)';%得到数字FIR 滤波器
[H,w]=freqz(h,[1],1024,'whole');
mag=abs(H); db=20*log10((mag+eps)/max(mag));
pha=angle(H); subplot(221) stem(n,h);title('h(n)');xlabel('n') subplot(222) plot(w/pi,mag) title('幅频特性');xlabel('w/pi')
axis([0,1,0,1.5]) subplot(223); wk=2*[0:1023]/1024; plot(w/pi,db);grid; title('衰减特性');xlabel('w/pi')
axis([0,1,-150,0]) subplot(224); plot(w,pha) axis([0,1,-4,4]) title('相频特性');xlabel('w')
2. 运行图形如下：
h(n)
10
2030n
衰减特性
40
0.501.51
幅频特性
0.5w/pi相频特性
1
0-50-100-150
420-2
0.5w/pi
1
-40
0.5w
1
二、汉宁窗
1. 编写程序如下：
wp=0.3*pi;ws=0.5*pi;As=40;Bt=ws-wp;N0=ceil(6.2*pi/Bt); N=N0+mod(N0+1,2);%确保N 为奇数
wc=(wp+ws)/2/pi;n=0:N-1;hn=fir1(N-1,wc,hanning(N))
fh0=fft(hn,1024);fh=20*log10(abs(fh0));wk=2*[0:1023]/1024; subplot(221);stem(n,hn) subplot(222);plot(wk,fh);grid; title('衰减特性');xlabel('w/pi')
axis([0,1,-150,0]) pha=angle(fh0); subplot(223) plot(wk,abs(fh0)) title('幅频特性');xlabel('w/pi')
axis([0,1,0,1.5]) subplot(224); plot(wk,pha) axis([0,1,-4,4]) title('相频特性');xlabel('w')
axis([0,1,-4,4]) 2. 运行结果如下： hn =
-0.0000 -0.0008 -0.0012 0.0023 0.0061 -0.0000 -0.0135
-0.0117 0.0160 0.0349 -0.0000 -0.0646 -0.0571 0.0900 0.2998 0.3999 0.2998 0.0900 -0.0571 -0.0646 -0.0000 0.0349 0.0160 -0.0117 -0.0135 -0.0000 0.0061 0.0023
-0.0012 -0.0008 -0.0000
3. 运行图形如下：
10
20
30
幅频特性
1.510.50
0.51
w/pi
三、 remez 函数设计FIR 低通滤波器
1. 编写程序如下： f=[2000,5500]; m=[0,1];Fs=16000; rp=1;rs=75;
dat1=(10^(rp/20)-1)/(10^(rp/20)+1); dat2=10^(-rs/20); rip=[dat2,dat1];
[M,fo,mo,w]=remezord(f,m,rip,Fs); N=M+1;
hn=remez(M,fo,mo,w); n=0:N-1;
wk=2*[0:1023]/1024; fh=fft(hn,1024);
fh1=20*log10(abs(fh)); subplot(121);stem(n,hn) axis([0,15,-0.6,0.6])
subplot(122);plot(wk,fh1);grid; title('|H(w)|');xlabel('w/pi') axis([0,1,-120,0])
衰减特性
0-50
-100-1500
0.51
w/pi相频特性
420-2
-40
0.51
w
运行图形如下：
5
10
15
|H(w)|
-20
-40-60-80-100-120
0.51
w/pi
2.