北理工数字信号处理实验报告.docx
- 文档编号:24109260
- 上传时间:2023-05-24
- 格式:DOCX
- 页数:44
- 大小:293.79KB
北理工数字信号处理实验报告.docx
《北理工数字信号处理实验报告.docx》由会员分享,可在线阅读,更多相关《北理工数字信号处理实验报告.docx(44页珍藏版)》请在冰豆网上搜索。
北理工数字信号处理实验报告
实验1利用DFT分析信号频谱
一、实验目的
1.加深对DFT原理的理解。
2.应用DFT分析信号的频谱。
3.深刻理解利用DFT分析信号频谱的原理,分析实现过程中出现的现象及解决方法。
二、实验设备与环境
计算机、MATLAB软件环境。
三、实验原理
1.DFT与DTFT的关系
有限长序列
的离散时间傅里叶变换
在频率区间
的N个等间隔分布的点
上的N个取样值可以由下式表示:
由上式可知,序列x(n)的N点DFTX(k),实际上就是x(n)序列的DTFT在N个等间隔频率点
上样本X(k)。
2.利用DFT求DTFT
法一:
由X(k)恢复出
的方法:
法二:
然而在实际MATLAB计算中,上述插值运算不见得是最好的办法。
由于DFT是DTFT的取样值,其相邻两个频率样本点的间距为2/N,所以如果我们增加数据的长度N,使得到的DFT谱线就更加精细,其包络就越接近DTFT的结果,这样就可以利用DFT来近似计算DTFT。
如果没有更多的数据,可以通过补零来增加数据长度。
3.利用DFT分析连续时间信号的频谱
采用计算机分析连续时间信号的频谱,第一步就是把连续时间信号离散化,这里需要进行两个操作:
一是采样,二是截断。
将利用DFT分析连续非周期信号频谱的步骤归纳如下:
(1)确定时域采样间隔T,得到离散序列x(n);
(2)确定截取长度M,得到M点离散序列
,这里
为窗函数。
(3)确定频域采样点数N,要求
。
(4)利用FFT计算离散序列的N点DFT,得到
。
(5)根据式(2-6)由
计算
采样点的近似值。
采用上述方法计算
的频谱,需要注意如下三个问题:
(1)频谱混叠。
如果不满足采样定理的条件,频谱会出现混叠误差。
对于频谱无限宽的信号,应考虑覆盖大部分主要频率分量的范围。
(2)栅栏效应和频谱分辨率。
使用DFT计算频谱,得到的结果只是N个频谱本值,样本值之间的频谱是未知的,像通过一个栅栏观察频谱,称为“栅栏效应”。
频谱分辨率与记录长度成反比,要提高频谱分表率,就要增加记录时间。
(3)频谱泄露。
对信号截断会把窗函数的频谱引入信号频谱,照成频谱泄露。
解决这个问题的主要办法是采用旁瓣小的窗函数,频谱泄露和窗函数均会引起误差。
因此,要合理选取采样间隔和截取长度,必要时还需考虑加适当的窗。
对于连续时间周期信号,我们在采用计算机进行计算时,也总是要进行截断,序列总是有限长的,仍然可以采用上述方法近似计算。
四、实验内容
1.已知x(n)={1,1,1,2},完成如下要求:
(1)计算其DTFT,并画出[-pi,pi]区间的波形。
源代码:
n=0:
3;
x=[2-111];
w=-pi:
0.01*pi:
pi;
X=x*exp(-1j*n'*w);
subplot(211);
plot(w,abs(X));
xlabel('\Omega/\pi');
title('Magnitude');
axis([-pipi04]);
subplot(212);
plot(w,angle(X)/pi);
xlabel('\Omega/\pi');
axis([-pipi-11]);
title('Phase');
实验结果:
(2)计算4点DFT,并把结果显示在
(1)所画的图形中。
源代码:
n=0:
3;
x=[2-111];
w=-pi:
0.01*pi:
pi;
X=x*exp(-1j*n'*w);
subplot(211);
plot(w,abs(X));
xlabel('\Omega/\pi');
title('Magnitude');
axis([-pipi04]);
holdon;
w1=-pi:
0.5*pi:
0.5*pi;
X1=x*exp(-1j*n'*w1);
stem(w1,abs(X1),'filled');
subplot(212);
plot(w,angle(X)/pi);
holdon;
stem(w1,angle(X1)/pi,'filled');
xlabel('\Omega/\pi');
axis([-pipi-14]);
title('Phase');
实验结果:
(3)对x(n)补零,计算64点DFT,并显示结果。
源代码:
n=0:
3;
n1=0:
63;
x=[2-111];
x1=[2-111zeros(1,60)];
w=-pi:
0.01*pi:
pi;
w1=-pi:
pi/32:
pi-pi/32;
X=x*exp(-1j*n'*w);
X1=x1*exp(-1j*n1'*w1);
subplot(211);
plot(w,abs(X));
xlabel('\Omega/\pi');
title('Magnitude');
axis([-pipi04]);
holdon;
stem(w1,abs(X1),'filled');
subplot(212);
plot(w,angle(X)/pi);
holdon;
stem(w1,angle(X1)/pi,'filled');
xlabel('\Omega/\pi');
axis([-pipi-14]);
title('Phase');
(4)根据实验结果,分析是否可以由DFT计算DTFT,如果可以,如何实现。
可以用DFT计算DTFT,由上图看出,N值越大,DFT的包络就越像DTFT的曲线图。
所以,我们可以很直接的得到这个答案,那就是,当我们取得N足够大的时候,我们便可以用DFT来表示DTFT了。
2.考察序列
(1)
时,用DFT估计
的频谱;将
补零加长到长度为100点序列
用DFT估计
的频谱。
要求画出相应波形。
源代码:
n1=0:
10;
x1=(cos(0.48*pi.*n1)+cos(0.52*pi.*n1));
X1=fft(x1);
subplot(211);
xlabel('\Omega/\pi');
title('Magnitude');
stem(n1,abs(X1),'filled');
n=0:
100;
ni=0:
9;
x=[(cos(0.48*pi.*ni)+cos(0.52*pi.*ni)),zeros([1,91])];
X=fft(x);
subplot(212);
xlabel('\Omega/\pi');
title('Magnitude');
stem(n/10,abs(X),'filled');
(2)
时,用DFT估计x(n)的频谱,并画出波形。
源代码:
n1=0:
100;
x1=cos(0.48*pi.*n1)+cos(0.52*pi.*n1);
X1=fft(x1);
xlabel('\Omega/\pi');
title('Magnitude');
stem(n1/100*2*pi,abs(X1),'filled');
(3)根据实验结果,分析怎样提高频谱分辨率。
加大n,n加的足够大的时候,信号的频谱看起来更像是原来实际信号的频谱。
3.已知信号
,其中
,
,
。
从
的表达式可以看出,它包含三个频率的正弦波,但是,从其时域波形(图E2-1)来看,似乎是一个正弦信号,利用DFT做频谱分析,确定适合的参数,使得到的频谱的频率分辨率符合需要。
源代码:
n=0:
0.1:
4.9;
x=0.15*sin(2*pi*n)+sin(4*pi*n)-0.1*sin(6*pi*n);
X=fft(x);
stem(n,abs(X),'filled');
4.利用DFT近似分析连续时间信号
的频谱(幅度谱)。
分析采用不同的采样间隔和截取长度进行计算的结果,并最终确定适合的参数。
源代码:
t=0:
0.1:
99.9;
w=-pi:
0.01:
pi;
x=exp(0.1*t);
X=x*exp(-1j*t'*w);
plot(w,abs(X));
五、实验中的主要结论、遇到的问题及解决方法,收获和体会。
主要结论:
1)可以用DFT来替代DTFT,但是有个条件,就是,在DFT里面去的值要足够的大,DFT的包络才可以代表DTFT的曲线。
2)当时域信号离散的时候,我们可以加大时间的取样,并且时间取得越长,离散信号的频谱就更接近真实连续信号的频谱。
遇到的问题:
在用MATLAB做DFT的时候,一度时间做不出来,信号出来的时候相位有延时,后来发现是时间取样不对。
收获和体会:
用MATLAB实实在在地看到了DFT和DTFT的区别和联系,对课本上的知识不再仅仅只是只能想象的阶段了,我们还可以,用计算机软件来仿真,这给了我们一种思路,一种从其他方面解决问题的方法。
实验三IIR数字滤波器设计
一、实验目的
1.掌握利用脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理及具体方法。
2.加深理解数字滤波器和模拟滤波器之间的技术指标转化。
3.掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的优缺点及使用范围。
二、实验设备与环境
计算机、MATLAB软件环境
三、实验原理
IIR滤波器具有无限长持续时间脉冲响应,而模拟滤波器一般都具有无限长的脉冲响应,
因此它与模拟滤波器相匹配。
IIR滤波器设计的基本方法就是先设计一个合适的模拟滤波器,
然后利用复值映射把模拟滤波器变换成数字滤波器。
1.数字滤波器和模拟滤波器的一些指标
通带、过渡带、通带响应中的容限、阻带的容限、通带波动、阻带衰减等等。
2.模拟原型滤波器
IIR滤波器设计方法由已有的模拟滤波器得到数字滤波器,我们将这些模拟滤波器称为原型滤波器。
常用的模拟原型滤波器有巴特沃斯(Butterworth)滤波器、切比雪夫(Chebyshev,Ⅰ型和Ⅱ型)滤波器和椭圆(Ellipse)滤波器等。
(1)巴特沃斯滤波器
(2)切比雪夫低通滤波器
a、切比雪夫I低通滤波器
b、切比雪夫II低通滤波器
(3)椭圆滤波器
(4)三种滤波器的比较
如果给定相同的设计指标,选用椭圆滤波器所要求的阶数N最低,切比雪夫滤波器次之,巴特沃斯滤波器最高;如果要求的阶数相同,切比雪夫滤波器的过渡带比巴特沃斯滤波器陡,椭圆滤波器的过渡带又比切比雪夫滤波器陡。
然而,从通带的相位响应来看,椭圆滤波器虽然提供了最优的幅度平方响应,但通带上的相位响应非线性较大,而巴特沃斯滤波器在通带上具有相当的线性相位,切比雪夫滤波器的相位特征介于两者之间。
所以在实际设计中,选用何种滤波器应视实际用途和指标要求而定。
3.模拟滤波器到数字滤波器的变换
从模拟滤波器到数字滤波器的变换就是要由Ha(s)进一步求得H(z),也就是由s平面到z平面的变换,这种变换要求数字滤波器能模仿模拟滤波器的特性。
我们最常用的从模拟滤波器到数字滤波器的变换方法有两种:
脉冲响应不变法和双线性变换法。
3.1脉冲响应不变法
基本原理:
从时域响应出发,使数字滤波器的单位脉冲响应h(n)模仿模拟滤波器的单位冲激响应ha(t),h(n)等于ha(t)的取样值。
3.2双线性变换法
基本原理:
从频率响应出发,直接使数字滤波器的频率响应逼近模拟滤波器的频率响应,进而求得数字滤波器的系统函数H(z)。
4.直接利用MATLAB函数设计IIR数字滤波器
前面我们介绍了butter函数、cheby1函数、cheby2函数、ellip函数设计模拟滤波器的方法,这些函数也可直接用于设计IIR数字滤波器,由于双线性变换法的优点,这些函数都采用双线性变换。
这些函数的具体使用方法如下:
[b,a]=butter(N,Wn)设计低通巴特沃斯数字滤波器
[b,a]=cheby1(N,Rp,Wn)设计低通切比雪夫Ⅰ型数字滤波器
[b,a]=cheby2(N,As,Wn,)设计低通切比雪夫Ⅱ型数字滤波器
[b,a]=ellip(N,Rp,As,Wn)设计低通椭圆数字滤波器
[b,a]=butter(N,Wn,'ftype')设计巴特沃斯数字滤波器
[b,a]=cheby1(N,Rp,Wn,'ftype')设计切比雪夫Ⅰ型数字滤波器
[b,a]=cheby2(N,As,Wn,'ftype')设计切比雪夫Ⅱ型数字滤波器
[b,a]=ellip(N,Rp,As,Wn,'ftype')设计椭圆数字滤波器
'ftype'表示数字滤波器的类型,可选选项有'high'、'low'、'stop',分别表示高通、低通、带阻。
上述函数中的参数N和Wn可以采用函数buttord、函数cheb1ord、函数cheb2ord、函数ellipord求得,这些函数能在已知设计指标时给出滤波器的阶数N和截止频率Wn。
这些函数的具体使用方法如下:
[N,Wn]=buttord(Wp,Wst,Rp,As)根据数字滤波器设计指标给出数字巴特沃斯滤波器的阶数N和截止频率Wn;
[N,Wn]=cheb1ord(Wp,Wst,Rp,As)根据数字滤波器设计指标给出数字切比雪夫Ⅰ型滤波器的阶数N和截止频率Wn;
[N,Wn]=cheb2ord(Wp,Wst,Rp,As)根据数字滤波器设计指标给出数字切比雪夫Ⅱ型滤波器的阶数N和截止频率Wn;
[N,Wn]=ellipord(Wp,Wst,Rp,As)根据数字滤波器设计指标给出数字椭圆滤波器的阶数N和截止频率Wn。
其中,Wp和Wst分别表示通带和阻带截止频率,Rp表示通带波动,As表示阻带衰减。
四、实验内容
设频率为
,设计数字低通滤波器,要求:
要求分别设计巴特沃斯、切比雪夫I型、切比雪夫II型和椭圆模拟原型滤波器,并分别结合脉冲响应不变法和双线性变换法进行设计。
结合实验结果,分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求,分析脉冲响应不变法和双线性变换法设计IIR数字滤波器的优缺点及使用范围。
(1)巴特沃斯
wp=2*pi*1000;
ws=2*pi*1500;
Rp=1;
As=15;
Fs=10000;
OmegaP=wp;
OmegaS=ws;
ep=sqrt(10^(Rp/10)-1);
Ripple=sqrt(1/(1+ep*ep));
Attn=1/(10^(As/20));
[N,OmegaC]=buttord(OmegaP,OmegaS,Rp,As,'s');
[z0,p0,k0]=buttap(N);
p=p0*OmegaC;z=z0*OmegaC;
k=k0*OmegaC^N;
ba=k*real(poly(z));
aa=real(poly(p));
[bd,ad]=impinvar(ba,aa,Fs);
[H,w]=freqz(bd,ad,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(bd,ad,w);
subplot(221);
plot(w/pi,mag);
title('A');
ylabel('|H|');
gridon;
subplot(2,2,3);plot(w/pi,db);
title('A(dB)');
xlabel('f(pi)');
ylabel('dB');
axis([0,1,-100,50]);
gridon;
subplot(222);
plot(w/pi,pha/pi);
title('Phase');
xlabel('');
ylabel('Phase(pi)');
axis([0,1,-1,1]);
gridon;
subplot(224);
plot(w/pi,grd);
title('t');
xlabel('pi)');
ylabel('x');
axis([0,1,0,20]);
gridon;
(2)切比雪夫I型
wp=0.2*pi;
ws=0.3*pi;
Rp=1;
As=15;
Fs=10000;
OmegaP=wp*Fs;
OmegaS=ws*Fs;
ep=sqrt(10^(Rp/10)-1);
Ripple=sqrt(1/(1+ep*ep));
Attn=1/(10^(As/20));
[N,OmegaC]=cheb1ord(OmegaP,OmegaS,Rp,As,'s');
[z0,p0,k0]=cheb1ap(N,Rp);
p=p0*OmegaC;z=z0*OmegaC;
k=k0*OmegaC^N;
ba=k*real(poly(z));
aa=real(poly(p));
[bd,ad]=impinvar(ba,aa,Fs);
[H,w]=freqz(bd,ad,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(bd,ad,w);
subplot(2,2,1);
plot(w/pi,mag);
title('|A|');
xlabel('');
ylabel('|H|');
axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,0.4,0.6,1]);
gridon;
set(gca,'XTickMode','manual','YTick',[0,Attn,Ripple,1]);
subplot(2,2,3);
plot(w/pi,db);
title('|A|(dB)');
xlabel('Frep(pi)');
ylabel('dB');
axis([0,1,-100,50]);
gridon;
subplot(2,2,2);
plot(w/pi,pha/pi);
title('Phase');
xlabel('');
ylabel('Phase(pi)');
axis([0,1,-1,1]);
gridon;
subplot(2,2,4);
plot(w/pi,grd);
title('layoff');
xlabel('Frep(pi)');
ylabel('x');
axis([0,1,0,20]);
gridon;
(3)切比雪夫II型
wp=0.2*pi;
ws=0.3*pi;
Rp=1;
As=15;
Fs=10000;
T=1/Fs;
OmegaP=(2/T)*tan(wp/2);
OmegaS=(2/T)*tan(ws/2);
ep=sqrt(10^(Rp/10)-1);
Ripple=sqrt(1/(1+ep*ep));
Attn=1/(10^(As/20));
[N,OmegaC]=cheb2ord(OmegaP,OmegaS,Rp,As,'s');
[z0,p0,k0]=cheb2ap(N,As);
p=p0*OmegaC;
z=z0*OmegaC;
k=k0*OmegaC^N;
ba0=real(poly(z0));
ba0=k0*ba0;
aa0=real(poly(p0));
ba=real(poly(z));
ba=k*ba;
aa=real(poly(p));
[bd,ad]=bilinear(ba,aa,Fs);
[bd1,ad1]=bilinear(ba0,aa0,Fs/OmegaC);
subplot(1,1,1)
[H,w]=freqz(bd,ad,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(bd,ad,w);
subplot(2,2,1);plot(w/pi,mag);
title('A');
xlabel('');
ylabel('|A|');
set(gca,'XTickMode','manual','XTick',[0,0.4,0.6,1]);
gridon;
set(gca,'XTickMode','manual','YTick',[0,Attn,Ripple,1]);
subplot(2,2,3);
plot(w/pi,db);
title('|A|(dB)');
xlabel('freq(rad/s)');
ylabel('dB');
gridon;
subplot(2,2,2);
plot(w/pi,pha/pi);
title('Phase');
xlabel('');ylabel('pi');
gridon;
subplot(2,2,4);
plot(w/pi,grd);
title('layoff');
xlabel('freq(rad/s)');
ylabel('Ñù±¾');
gridon;
(4)椭圆模拟原型
wp=0.2*pi;
ws=0.3*pi;
Rp=1;
As=15;
Fs=10000;
T=1/Fs;
OmegaP=(2/T)*tan(wp/2);
OmegaS=(2/T)*tan(ws/2);
ep=sqrt(10^(Rp/10)-1);
Ripple=sqrt(1/(1+ep*ep));
Attn=1/(10^(As/20));
[N,OmegaC]=ellipord(OmegaP,OmegaS,Rp,As,'s');
[z0,p0,k0]=ellipap(N,Rp,As);
p=p0*OmegaC;
z=z0*OmegaC;
k=k0*OmegaC^N;
ba0=real(poly(z0));
ba0=k0*ba0;
aa0=real(poly(p0));
ba=real(poly(z));
ba=k*ba;
aa=real(poly(p));
[bd,ad]=bilinear(ba,aa,Fs);
[bd1,ad1]=bilinear(ba0,aa0,Fs/OmegaC);
subplot(1,1,1)
[H,w]=freqz(bd,ad,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(bd,ad,w);
subplot(2,2,1);plot(w/pi,mag);
title('A');
xlabel('');
ylabel('|A|');
set(gca,'XTickMod
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北理工 数字信号 处理 实验 报告