数字信号处理实验报告

数字信号处理实验报告
⼀、
课程设计(综合实验)的⽬的与要求
⽬的与要求：
1．掌握《数字信号处理基础》课程的基本理论； 2．掌握应⽤MATLAB 进⾏数字信号处理的程序设计；实验内容：
已知低通数字滤波器的性能指标如下：
0.26p ωπ=，0.75dB p R =，0.41s ωπ=，50dB s A =
要求：
1. 选择合适的窗函数，设计满⾜上述指标的数字线性相位FIR 低通滤波器。

⽤⼀个图形窗
⼝，包括四个⼦图，分析显⽰滤波器的单位冲激响应、相频响应、幅频响应和以dB 为纵坐标的幅频响应曲线。

2. ⽤双线性变换法，设计满⾜上述指标的数字Chebyshev I 型低通滤波器。

⽤⼀个图形窗⼝，
包括三个⼦图，分析显⽰滤波器的幅频响应、以dB 为纵坐标的幅频响应和相频响应。

3. 已知模拟信号
1234()2sin(2)5sin(2)8cos(2)7.5cos(2)x t f t f t f t f t ππππ=+++
其中10.12f kHz =，2 4.98f kHz =，3 3.25f kHz =，4 1.15f kHz =，取采样频率
10s f kHz =。

要求：
(1) 以10s f kHz =对()x t 进⾏取样，得到()x n 。

⽤⼀个图形窗⼝，包括两个⼦图，分别显
⽰()x t 以及()x n (0511n ≤≤)的波形；
(2) ⽤FFT 对()x n 进⾏谱分析，要求频率分辨率不超过5Hz 。

求出⼀个记录长度中的最少点
数x N ，并⽤⼀个图形窗⼝，包括两个⼦图，分别显⽰()x n 以及()X k 的幅值； (3) ⽤要求1中设计的线性相位低通数字滤波器对()x n 进⾏滤波，求出滤波器的输出1()y n ，
并⽤FFT 对1()y n 进⾏谱分析，要求频率分辨率不超过5Hz 。

求出⼀个记录长度中的最少点数1y N ，并⽤⼀个图形窗⼝，包括四个⼦图，分别显⽰()x n （01x n N ≤≤-）、()X k 、
1()y n （101y n N ≤≤-）和1()Y k 的幅值；
(4) ⽤要求2中设计的Chebyshev 低通数字滤波器对()x n 进⾏滤波，求出滤波器的输出
2()y n ，并⽤FFT 对2()y n 进⾏谱分析，要求频率分辨率不超过5Hz 。

求出⼀个记录长度
中的最少点数2y N ，并⽤⼀个图形窗⼝，包括四个⼦图，分别显⽰()x n （01x n N ≤≤-）、()X k 、2()y n （201y n N ≤≤-）和2()Y k 的幅值。

备注：
(1) 要求编写⼀个主程序，完成上述问题。

(2) 要求各结果图均标出图题以及横纵轴的名称。

(3) 要求给主要语句加上标注。

设计思路：
1．在任务⼀中，最⼩阻带衰减As=50dB，则就此可以确定所要选择的窗函数为哈明窗。

2. 对于任务⼆，根据通带频率wp，阻带频率ws,通带衰减Rp，阻带衰减As，计算出所设计
的切⽐雪夫滤波器阶数，极点等重要参数。

3．对连续时间信号进⾏采样，得到()
x n，根据要求绘制出⼦图。

4．频率分辨率等于采样频率除以采样点数，据此可以计算出⼀个记录中的最少点数，执⾏函数fft，可以将x（n）经过傅⾥叶变换求出X(k)。

5．⽤线性相位低通数字滤波器对()
x n进⾏滤波，在设计低通数字滤波器时已经求出了单位冲击响应序列h（n），那么将采样得到的序列x（n）与h（n）进⾏卷积，就可以得到
滤波器的输出
1()
y n，经过fft变换后，可以对
1()
y n进⾏频谱分析。

6．⽤Chebyshev低通数字滤波器对()
x n进⾏滤波时，由于设计时能够将切⽐雪夫滤波器系统函数（z域，双线性变换变换的表达式）的零极点均已求出，那么就可以调⽤filter函
数求出
2()
y n，最后将得到滤波器的输出进⾏fft变换，对
2()
y n进⾏频谱分析，按要求画出相应的⼦图。

实验结果：
数字线性相位FIR低通滤波器
数字Chebyshev I型低通滤波器
以10s f kHz 对()x t 进⾏取样
⽤FFT 对()x n 进⾏谱分析
x n进⾏滤波
线性相位低通数字滤波器对()
x n进⾏滤波
Chebyshev低通数字滤波器对()
程序代码：
主程序：
clear;
%任务⼀：
%窗⼝函数的设计
wp=0.26*pi;
ws=0.41*pi;
tr_width=ws-wp;%带宽
N1=ceil(6.6*pi/tr_width)+1;%阶数
n=0:N1-1;
wc=(ws+wp)/2;%理想低通的截⽌频率
hd=ideal_lp(wc,N1);%理想低通的冲击响应
w_ham=(hamming(N1))';
h=hd.*w_ham;%FIR滤波器的冲击响应
[db,mag,pha,grd,w]=freqz_m(h,[1]);
delta_w=2*pi/1000;
%作图
figure(1);
subplot(2,2,1);
stem(n,h,'r');
axis([0 46 -0.1 0.4]);
title('单位冲击响应');
subplot(2,2,2);
plot(w/pi,pha,'r');
axis([0 1 -4 4]);
title('相频响应');
subplot(2,2,3);
plot(w/pi,mag,'r');
axis([0 0.8 0 1.1]);
title('幅频响应');
subplot(2,2,4);
plot(w/pi,db,'r');
axis([0 0.8 -100 10]);
title('幅频响应(dB)');xlabel('以\pi为单位的频率');ylabel('对数幅度/dB'); %任务⼆; %切⽐雪夫滤波器的设计
Rp=0.75;
As=50;
T=1;
Fs=1/T;
omegap=(2/T)*tan(wp/2);
omegas=(2/T)*tan(ws/2);
[cs,ds]=cheb(omegap,omegas,Rp,As);%调⽤函数求解拉普拉斯变换的分⼦与分母多项式的系数数组
[b,a]=bilinear(cs,ds,Fs);%做双线性变换
[db,mag,pha,grd,w]=freqz_m(b,a);%频域上采样并记录幅频和相频响应
%作图
figure(2);
subplot(3,1,1);
plot(w/pi,mag);
title('幅度响应');grid;
axis([0 0.6 0 1.1]);ylabel('幅度');xlabel('以\pi为单位的频率');
subplot(3,1,2);
plot(w/pi,db);title('幅频响应(dB)');grid;
axis([0 0.6 -60 10]);xlabel('以\pi为单位的频率');ylabel('对数幅度/dB'); subplot(3,1,3);
plot(w/pi,pha);title('相频响应');grid;
axis([0 0.6 -4 4]);xlabel('以\pi为单位的频率');ylabel('相位');
%任务三：
fs=10000;
f(1)=120;
f(2)=4980;
f(3)=3250;
f(4)=1150;
Nq=512;
nq=0:Nq-1;
xg=2*cos(2*pi*f(1)*nq/fs)+5*sin(2*pi*f(2)*nq/fs)+8*cos(2*pi*f(3)*nq/f s)+7.5*cos(2*pi*f(4)*nq/fs);%模拟信号采样的表达式%作图
figure(3);
subplot(2,1,1);
fplot('xat',[0,511],'c');
title('x(t)');
axis([0,512,-24,24]);
subplot(2,1,2);
stem(nq,xg,'r');
title('x(n)');
axis([0,512,-25,25]);
%任务四：
Nx=fs/5;
nq1=0:Nx-1;
xg1=2*cos(2*pi*f(1)*nq1/fs)+5*sin(2*pi*f(2)*nq1/fs)+8*cos(2*pi*f(3)*n q1/fs)+7.5*cos(2*pi*f(4)*nq1/fs);
fprintf('\n***最少点数Nx=%2.0f\n',Nx);
Xkc=fft(xg1,Nx);%fft
Akc=abs(Xkc);%幅度
%作图
figure(4);
subplot(2,1,1);
stem(nq1,xg1);
axis([0,Nx-1,-25,25]);
title('x(n)');
subplot(2,1,2);
stem(nq1,Akc,'r');
title('X(k)');
%任务五：
Ny1=fs/5;
fprintf('\n***y1(n)最少点数Ny1=%2.0f\n',Ny1); [y1]=conv(h,xg1);%求卷积
Y0=fft(y1);%fft
Y1=abs(Y0);
%作图
figure(5);
subplot(4,1,1);
stem(nq1,xg1,'r');
axis([0,Ny1-1,-25,25]);
title('x(n)');
subplot(4,1,2);
stem(nq1,Akc,'pc');
title('X(k)');
subplot(4,1,3);
n2=1:N1+Ny1-1;
stem(n2,y1);
axis([0,2100,-12,12]);
title('y1(n)');
subplot(4,1,4);
stem(n2,Y1,'c');
axis([0,2100,0,7500]);
title('Y1(k)');
%任务六：
Ny2=fs/5;
n3=1:Ny2;
fprintf('\n***y2(n)最少点数Ny2=%2.0f\n',Ny2); y2=filter(b,a,xg1);%滤波
Y2=abs(fft(y2));
%作图
figure(6);
subplot(4,1,1);
stem(nq1,xg1,'r');
axis([0,Ny1-1,-25,25]);
title('x(n)');
subplot(4,1,2);
stem(nq1,Akc,'pc');
title('X(k)');
subplot(4,1,3);
stem(n3,y2);
axis([0,2000,-12,12]);
title('y2(n)');
subplot(4,1,4);
stem(n3,Y2,'c');
title('Y2(k)');
调⽤的函数：
1.function [M,x]=cheb(wp,ws,Rp,As)
ep=sqrt(10^(Rp/10)-1);
omegac=2*tan(wp/2);
omegar=ws/wp;
N=ceil(acosh(sqrt(10^(As/10)-1)/ep)/acosh(omegar)); af=ep^(-1)+sqrt(ep^(-2)+1); a=0.5*(af^(1/N)-af^(-1/N));
b=0.5*(af^(1/N)+af^(-1/N));
A=-a*omegac;
B=b*omegac;
D=1;
for i=0:N-1
p(i+1)=A*sin((2*i+1)/(2*N)*pi);
q(i+1)=B*cos((2*i+1)/(2*N)*pi);
w(i+1)=p(i+1)+j*q(i+1);
D=D*w(i+1);
end
fprintf('\n***切⽐雪夫滤波器的阶数N=%2.0f\n',N);
fprintf('\n***切⽐雪夫滤波器的极点：\n');
w
for k=1:N
u(k)=abs(w(k));
end
u
if max(u)<1
fprintf('\n***系统稳定\n');
end
if(rem(N,2)==1)
L=abs(D);
else
L=abs(1/sqrt(1+ep^2)*D);
end
fprintf('\n***切⽐雪夫滤波器的系数：\n');
L
M=[zeros(1,N-1),L];
x=real(poly(w));
2.function [db,mag,pha,grd,W]=freqz_m(b,a)
[H,W]=freqz(b,a,1000,'whole');
H=(H(1:1:501))';
W=(W(1:1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
grd=grpdelay(b,a,W);
grd=grd';
3.function hd=ideal_lp(wc,M)
alpha=(M-1)/2;
n=0:M-1;
m=n-alpha+eps;
hd=sin(wc*m)./(pi*m);
4.function y=xat(t)
y=2*cos(2*pi*120*t)+5*sin(2*pi*4980*t)+8*cos(2*pi*3250*t)+7.5*cos(2*p i*1150*t);
三、课程设计（综合实验）总结或结论
通过本次的课程设计使我对FIR与IIR滤波器有了更加深⼊地了解：IIR滤波器的系统函数是有理分式，分母对应于反馈⽀路，它是递归结构的系统，所以只有当所有极点都在单位圆内时滤波器才稳定。

FIR滤波器的系统函数是多项式，是⾮递归结构，它只在原点处有
⼀个N阶极点，因⽽系统稳定。

另外在本次的课设过程中我学会了MATLAB的使⽤，掌握了MATLAB内部函数的调⽤⽅法，通过查询MATLAB函数⼿册使我能够在遇到问题时准确清楚地了解函数中各个参数代表的的实际意义，这在设计滤波器过程中⼤⼤提⾼了设计效率。

最后在本次的实践过程中让我对数字信号处理部分知识点理解也更加透彻，总体收获颇丰。

四、参考⽂献
数字信号处理姚天任，江太辉第三版华中科技⼤学出版社，2007。