数字信号处理.docx
- 文档编号:4851383
- 上传时间:2022-12-10
- 格式:DOCX
- 页数:23
- 大小:138.24KB
数字信号处理.docx
《数字信号处理.docx》由会员分享,可在线阅读,更多相关《数字信号处理.docx(23页珍藏版)》请在冰豆网上搜索。
数字信号处理
前言
错误!
未找到目录项。
“数字信号处理”是一门理论和实践密切结合的课程,为了深入地掌握课程内容,应当在学习理论的同时,作习题和上机实验。
上机实验不仅可以帮助读者深入地理解和消化基本理论,而且能锻炼初学者的独立解决问题的能力。
本课程根据课程重点编写了五个实验,供学生使用或参考。
由于数字信号处理实验的主要目的是验证数字信号处理的有关理论,进一步理解巩固所学理论知识,所以,对实验用算法语言不作任何限制。
为了提高实验效率,我们提倡学生选用编程效率比C语言高好几倍的MATLAB语言,按照指导书的要求,上机编程完成实验。
实验一用FFT进行谱分析实验
一、实验目的:
(1)用FFT进行谱分析,了解fft.m文件的各参数及使用方法;
(2)学习提高频率分辨率的方法,加深对栅栏效应和频谱泄漏等概念的理解。
二、实验设备:
计算机,MATLAB软件。
三、实验原理:
离散傅里叶变换(DFT)可以用快速傅里叶变换(FFT)算法来计算。
在MATLAB信号处理工具箱中,提供了函数fft()、ifft()分别求解离散傅里叶变换与逆变换。
调用格式如下:
Xk=fft(x)
Xk=fft(x,N)
表示计算信号x的快速离散傅里叶变换Xk。
当x的长度N为2的整数次方时,采用基2算法,否则采用较慢的分裂基算法。
当length(x)>N时,截断x,否则补零。
x=ifft(Xk)
x=ifft(Xk,N)
表示计算Xk的逆离散傅里叶变换。
1.用FFT进行谱分析
用FFT的结果分析x(t)=cos(2π×50t)+0.5cos(2π×150t)+0.3cos(2π×250t)的频谱。
t=0:
0.02/64:
0.04;
f1=50;
y1=cos(2*pi*f1*t)+0.5*cos(2*pi*3*f1*t)+0.3*cos(2*pi*5*f1*t);
subplot(311);plot(t,y1);
t=0:
0.02/16:
0.02-0.02/16;
f=cos(2*pi*f1*t)+0.5*cos(2*pi*3*f1*t)+0.3*cos(2*pi*5*f1*t);
F_1024=2*abs(fft(f,16))/16;
k=0:
1:
15;
subplot(312);stem(k,abs(F_1024));%由于栅栏效应,只能看到16条谱线
axis([0,16,0,1.5])
F_1024=2*abs(fft(f,1024))/16;%补零减小栅栏效应,可以得到连续频谱
L=0:
1023;
subplot(313);plot(L/1023,abs(F_1024));
set(gca,'xtick',[0,0.0625,0.125,0.1875,0.25,0.3125,0.375,0.4375,0.5,0.5625,0.625,0.6875,0.75,0.8125,0.875,0.9375,1])%频率刻度为归一化频率
运行该程序,结果显示如下:
2.用FFT进行谱分析中的观测时间的选取
改变观测时间可以提高频率分辨率。
对x(t)=cos(2π×50t)+0.5cos(2π×75t)进行DFT。
t=0:
0.02/64:
0.04;
f1=50;f2=75;
y1=cos(2*pi*f1*t);
y2=0.5*cos(2*pi*f2*t);
subplot(411);plot(t,y1);
holdon;
subplot(411);plot(t,y2);
holdon;
k=0:
1:
7;
f=cos(pi/4*k)+0.5*cos(3*pi/8*k);
F_1024=2*abs(fft(f,1024))/8;
L=0:
1023;
subplot(412);plot(L/1023,abs(F_1024));
%axis([0,1024,0,2])
set(gca,'xtick',[0,0.125,0.25,0.375,0.5,0.625,0.75,0.875,1])
holdon;
k=0:
1:
15;
f=cos(pi/4*k)+0.5*cos(3*pi/8*k);
F_1024=2*abs(fft(f,1024))/16;
L=0:
1023;
subplot(413);plot(L/1024,abs(F_1024));
set(gca,'xtick',[0,0.125,0.25,0.375,0.5,0.625,0.75,0.875,1])
holdon;
k=0:
1:
31;
f=cos(pi/4*k)+0.5*cos(3*pi/8*k);
F_1024=2*abs(fft(f,1024))/32;
L=0:
1023;
subplot(414);plot(L/1024,abs(F_1024));
set(gca,'xtick',[0,0.125,0.25,0.375,0.5,0.625,0.75,0.875,1])
holdon;
运行该程序,结果显示如下:
四、实验内容与步骤:
1.将例题程序输入计算机,运行并对各图中的波形形状进行解释,并说明理由。
2.修改用FFT进行谱分析的程序,将采样点数增加为32(截取2个周期)个,重新用FFT进行谱分析。
3.修改用FFT进行谱分析的程序,将采样点数增加为64(截取4个周期)个,重新用FFT进行谱分析。
实验二用脉冲响应不变法设计IIR数字滤波器
一、实验目的:
(1)了解巴特沃思模拟低通滤波器、切比雪夫Ⅰ型、Ⅱ型模拟低通滤波器函数的使用方法;
(2)学习脉冲响应不变法设计IIR数字滤波器的设计方法。
二、实验设备:
计算机,MATLAB软件。
三、实验原理:
MATLAB提供了设计巴特沃思模拟低通滤波器、切比雪夫Ⅰ型、Ⅱ型模拟低通滤波器、椭圆模拟低通滤波器的函数,调用格式为
[z,p,k]=buttap(n)
[z,p,k]=cheb1ap(n,Rp)
[z,p,k]=cheb2ap(n,Rs)
[z,p,k]=ellipap(n,Rp,Rs)
其中,n为滤波器的阶次;z、p、k分别为滤波器传递函数的零点、极点和增益;Rp为通带波纹,Rs为阻带衰减。
MATLAB还提供了模拟滤波器的频率变换函数,可以借助模拟低通滤波器的系统函数,经过适当的频率变换,得到高通、带通、带阻滤波器的系统函数。
调用格式为
(1)低通到低通的变换
[bt,at]=lp2lp(b,a,wp)
[At,Bt,Ct,Dt]=lp2lp(A,B,C,D,wp)
(2)低通到高通的变换
[bt,at]=lp2hp(b,a,wp)
[At,Bt,Ct,Dt]=lp2hp(A,B,C,D,wp)
(3)低通到带通的变换
[bt,at]=lp2bp(b,a,w0,Bw)
[At,Bt,Ct,Dt]=lp2bp(A,B,C,D,w0,Bw)
(4)低通到带阻的变换
[bt,at]=lp2bs(b,a,w0,Bw)
[At,Bt,Ct,Dt]=lp2bs(A,B,C,D,w0,Bw)
其中,b、a和bt、at分别为变换前和变换后的系统函数的分子和分母系数向量;第二种格式是系统状态空间形式。
wp为通带频率,w0为中心频率,Bw为带宽。
冲激响应不变法设计IIR数字滤波器
MATLAB提供了使用脉冲响应不变法设计IIR数字滤波器的函数impinvar(),调用格式为
[bz,az]=impinvar(b,a,fs)
[bz,az]=impinvar(b,a,fs,tol)
将模拟滤波器(b,a)变换成数字滤波器(bz,az)。
其中,fs表示采样频率,单位为Hz,默认值为1。
Tol表示区分多重极点的程度,默认值为0.1%。
MATLAB提供了freqz()函数可方便地绘制出系统的频率特性,调用格式为
freqz(b,a)
[h,w]=freqz(b,a,n)
[h,f]=freqz(b,a,n,Fs)
h=freqz(b,a,w)
h=freqz(b,a,f,Fs)
[h,w]=freqz(b,a,n,’whole’)
[h,f]=freqz(b,a,n,’whole’,Fs)
其中,b、a分别为系统函数的分子、分母系数向量;n为频率的计算点数,常取2的整数次幂;横坐标为数字角频率ω,范围为0到π。
[h,w]=freqz(b,a,n)自动设定n个频率点来计算频率特性h,n个频率点均匀地分布在0到π,这n个频率值记录在w中。
[h,f]=freqz(b,a,n,Fs)自动设定n个频率点来计算频率特性h,n个频率点均匀地分布在0~Fs/2中,这n个频率值记录在w中。
[h,w]=freqz(b,a,n,’whole’)表示在0到2π中均匀选取n个点计算频率特性。
[h,f]=freqz(b,a,n,’whole’,Fs)表示在0到Fs中均匀选取n个点计算频率特性。
h=freqz(b,a,w)计算在向量w中指定的频率处的频率特性。
不带输出变量的freqz函数,将在当前图形窗口中绘制出幅频和相频曲线。
例:
采用脉冲响应不变法设计一个低通切比雪夫Ⅰ型数字滤波器,技术指标为:
通带频率是300Hz,阻带频率是500Hz,采样频率是1000Hz,通带波纹
dB,阻带衰减
dB。
程序如下:
wap=2*pi*300;was=2*pi*500;%通带、阻带截止频率
rp=0.3;rs=60;%通带、阻带衰减
fs=1200;%采样频率
[N,wn]=cheb1ord(wap,was,rp,rs,'s');%选择滤波器的最小阶数
[z,p,k]=cheb1ap(N,rp);%创建切比雪夫Ⅰ型模拟低通滤波器,z、p、k分别为滤波器传递函数的零点、极点和增益
[b,a]=zp2tf(z,p,k);%系统零极点增益模型转换成系统函数模型Ha(p)
[bt,at]=lp2lp(b,a,wn);%低通到低通的变换即Ha(p)去归一化转换为Ha(s)
[bz,az]=impinvar(bt,at,fs);%脉冲响应不变法将模拟滤波器转换为数字滤波器
[h,f]=freqz(bz,az,512,fs);%求幅频响应
hdb=20*log10(abs(h));
subplot(1,2,1);plot(f,abs(h))
xlabel('频率/Hz');ylabel('幅值')
subplot(1,2,2);plot(f,hdb)
xlabel('频率/Hz');ylabel('幅值(db)')
执行结果如图所示。
四、实验内容与步骤:
1.将上例输入计算机运行,分析结果
2.将上例中的滤波器换成巴特沃思模拟低通滤波器,写出程序并运行,分析结果。
实验三用双线性变换法设计IIR数字滤波器
一、实验目的:
(1)了解butter.m文件的各参数及使用方法;
(2)学习双线性变换法设计IIR数字滤波器的设计方法。
二、实验设备:
计算机,MATLAB软件。
三、实验原理:
在MATLAB中,butter.m文件可用来直接设计巴特沃思数字滤波器,实际上它是把buttord、buttap、lp2lp及bilinear等文件都包含了进去,从而使设计过程更简捷,其调用格式是:
①[B,A]=butter(N,Wn);
②[B,A]=butter(N,Wn,’high’);
③[B,A]=butter(N,Wn,’stop’);
④[B,A]=butter(N,Wn,’s’);
格式①~③用来设计数字滤波器,所以B,A分别是H(z)分子、分母多项式的系数向量,Wn是通带截止频率,范围在0~1之间,1对应抽样频率的一半。
若Wn是标量,则格式①用来设计低通数字滤波器;若Wn是1x2的向量,则格式①用来设计数字带通滤波器,格式②用来设计数字高通滤波器,格式③用来设计数字带阻滤波器,显然,这时的Wn是1x2的向量;格式④用来设计模拟滤波器。
例:
试用双线性变换法设计一数字低通滤波器,给定的技术指标为fp=75Hz,αp=3dB,fs=225Hz,αs=20dB,采样频率为600Hz,指定模拟滤波器采用巴特沃思低通滤波器。
解由于2π对应600Hz,所以ωp=2π×75/600=0.25π,ωs=2π×225/600=0.75π。
(1)将数字滤波器的技术指标转换为模拟滤波器的技术指标。
由于在变换过程中,系数2/T被约掉,实际上变换结果与T无关,为了简便,由式(7-48)计算技术指标时省去系数2/T,得模拟频率为
rad/s
rad/s
(2)设计巴特沃思低通滤波器。
确定阶数N
取N=2,查表7-1得
去归一化得
(3)用双线性变换法求H(z)。
用Matlab设计上例要求的IIR滤波器的程序如下:
fp=75;fs=225;%通带、阻带截止频率
f=600;%采样频率
rp=3;rs=20;%通带、阻带衰减
wp=2*pi*fp/f;ws=2*pi*fs/f;%通带、阻带截止数字频率
wap=tan(wp/2)%通带截止模拟频率,相当于T取2
was=tan(ws/2)%阻带截止模拟频率
[n,wn]=buttord(wap,was,rp,rs,'s');%’s’是确定巴特沃思模拟滤波器阶次和3dB截止模拟频率
[z,p,k]=buttap(n);%设计归一化巴特沃思模拟低通滤波器,z极点,p零点和k增益
[bp,ap]=zp2tf(z,p,k)%转换为Ha(p)表示,bp分子系数,ap分母系数
[bs,as]=lp2lp(bp,ap,wap)%Ha(p)去归一化转换为Ha(s)表示,bs分子系数,as分母系数
[bz,az]=bilinear(bs,as,1/2)%双线性变换为H(z),bz分子系数,az分母系数,采样频率取1/2
freqz(bz,az,32,600)%画数字滤波器的频率响应和相位响应
运行结果如下:
wap=0.41421356237310
was=2.41421356237309
bp=001
ap=1.000000000000001.414213562373091.00000000000000
bs=0.17157287525381
as=1.000000000000000.585786437626900.17157287525381
.0976********
az=1.00000000000000-0.942809041582060.33333333333333
可以看出计算结果与上例一致。
四、实验内容与步骤:
1.将上例输入计算机运行,分析结果
2.将上例中的滤波器换成切比雪夫Ⅰ型模拟低通滤波器,写出程序并运行,分析结果。
实验四常用窗函数的特性和FIR数字滤波器设计实验
一、实验目的:
(1)了解常用窗函数的特性;
(2)学习FIR数字低通滤波器的设计方法。
二、实验设备:
计算机,MATLAB软件。
三、实验原理与步骤:
窗函数法的一般设计步骤为:
(1)给定要求的频率响应函数
;
(2)计算单位样值响应
;
(3)根据过渡带宽及阻带最小衰减的要求,选定窗的大小N,N可通过多次尝试后进行最优确定;
(4)根据所选择的合适的窗函数
来修正
,得到所设计的FIR滤波器的单位样值响应
,n=0,1,…,N-1。
1.窗函数特性
MATLAB提供的窗函数主要有:
w=boxcar(n)%矩形窗
w=triang(n)%三角窗
w=hanning(n)%汉宁窗
w=hamming(n)%汉明窗
w=blackman(n)%布莱克曼窗
w=kaiser(n,beta)%凯塞窗
w=chebwin(n,r)%切比雪夫窗
其中n是窗函数的长度;w是由窗函数的值组成的n阶向量。
画三角窗、汉宁窗、汉明窗、布莱克曼窗窗函数的matlab程序:
N=31;t=(0:
N-1);
w0=boxcar(N);%矩形窗
w1=bartlett(N);%三角窗
w2=hanning(N);%汉宁窗
w3=hamming(N);%汉明窗
w4=blackman(N);%布莱克曼窗
%若画离散窗函数,用语句stem(w)
figure
plot(t,w1,'-k',t,w2,'-ok',t,w3,'-*k',t,w4,'-+k');
legend('三角窗','汉宁窗','汉明窗','布莱克曼窗');
画三角窗、汉宁窗、汉明窗、布莱克曼窗幅度特性的matlab程序:
figure
[h0,f]=freqz(w0,1,512,2);
[h1,f]=freqz(w1,1,512,2);
[h2,f]=freqz(w2,1,512,2);
[h3,f]=freqz(w3,1,512,2);
[h4,f]=freqz(w4,1,512,2);
subplot(221);
H1=20*log10(abs(h1)/max(h1))
plot(f,H1);grid
axis([0,1,-100,0]);
title('三角窗');
subplot(222);
H2=20*log10(abs(h2)/max(h2))
plot(f,H2);grid
axis([0,1,-100,0]);
title('汉宁窗');
subplot(223);
H3=20*log10(abs(h3)/max(h3))
plot(f,H3);grid
axis([0,1,-100,0]);
title('汉明窗');
subplot(224);
H4=20*log10(abs(h4)/max(h4))
plot(f,H4);grid
axis([0,1,-100,0]);
title('布莱克曼窗');
2.设计一个FIR低通滤波器
基于窗函数的FIR(有限冲激响应)滤波器设计——标准频率响应。
格式:
b=fir1(n,Wn)
b=fir1(n,Wn,'ftype')
b=fir1(n,Wn,Window)
b=fir1(n,Wn,'ftype',Window)
说明:
fir1函数以经典方法实现加窗线性相位FIR数字滤波器设计,它可设计出标准的低通、带通、高通和带阻滤波器(具有任意频率响应的加窗滤波器由fir2函数设计)。
b=fir1(n,Wn)可得到n阶低通FIR滤波器,滤波器系数包含在b中,这可表示成
b(z)=b
(1)+b
(2)z-1+…·-+b(n十1)z-n
这是一个截止频率为Wn的Hamming(汉明)加窗线性相位滤波器,0≤Wn≤1,Wn=1相应于0.5fs。
当Wn=[W1W2]时,fir1函数可得到带通滤波器,其通带为W1<ω<W2。
b=fir1(n,Wn,'ftype')可设计高通和带阻滤波器,由ftype决定:
·当ftype=high时,设计高通FIR滤波器;
·当ftype=stop时,设计带阻FIR滤波器。
在设计高通和带阻滤波器时,firl函数总是使用阶为偶数的结构,因此当输入的阶次为奇数时,firl函数会自动将阶次加1。
这是因为对奇次阶的滤波器,其在Nyquist频率处的频率响应为零,因此不适合于构成高通和带阻滤波器。
b=fir1(n,Wn,Window)则利用列矢量Window中指定的窗函数进行滤波器设计,Window长度为n+1。
如果不指定Window参数,则fir1函数采用Hamming窗。
b=firl(n,Wn,'ftype',Window)可利用ftype和Window参数,设计各种加窗的滤波器。
由fir1函数设计的FIR滤波器的群延迟为n/2。
例如设计一24阶FIR带通滤波器,通带为0.35<ω<0.65。
其程序如下
b=fir1(48,[0.350.65]);
freqz(b,1,512)
例:
设计FIR低通滤波器,通带允许起伏1dB,要求通带边缘频率fp=100Hz,阻带边缘频率fs=200Hz。
阻带最小衰减大于40dB。
并对信号
u(t)=cos(2π×50t)+0.8cos(2π×250t)+0.1cos(2π×350t)+0.05cos(2π×450t)
进行滤波。
解用理想低通作为逼近滤波器,有
汉宁窗的通带最大衰减为0.11dB,阻带最小衰减为-44dB,选择汉宁窗截断可以满足要求。
用Matlab设计上例要求的FIR滤波器的程序如下:
t=0:
0.02/32:
0.06;
f1=50;
u1=cos(2*pi*f1*t)+0.8*cos(2*pi*5*f1*t)+0.1*cos(2*pi*7*f1*t)+0.05*cos(2*pi*9*f1*t);
b=fir1(32,100/800)%求滤波器的单位脉冲响应
f=0:
1:
800;
h=freqz(b,1,f,1600);%求滤波器的频率响应
y=conv(u1,b);%计算滤波器的输出--卷积和
subplot(411);
plot(t,u1);
axis([0,0.06,-2,2]);
title('带高次谐波的信号');
holdon;
n=0:
32
subplot(412);
stem(n,b);
axis([0,32,-0.05,0.15])
title('滤波器的单位脉冲响应');
holdon;
subplot(413);
plot(f,abs(h));
title('滤波器的频率响应');
holdon;
subplot(414);
plot(t(1:
96),y(1:
96),':
');
%plot(ny,y);
title('通过滤波器后的信号');
holdon;
四、实验内容与思考题:
1.上例中,为什么滤波器输出信号的第一个周期不是正弦波?
2.设计一个只让5次谐波通过的FIR带通滤波器,Wn=[W1W2]=[200300]。
3.将上例输入计算机运行,分析结果。
4.画出三角窗、汉宁窗、汉明窗、布莱克曼窗幅度特性。
实验五切比雪夫最佳一致逼近FIR数字低通滤波器设计
一、实验目的:
(1)了解remez.m文件的各参数及使用方法;
(2)学习切比雪夫最佳一致逼近FIR数字低通滤波器的设计方法。
二、实验设备:
计算机,MATLAB软件。
三、实验原理:
在MATLAB中,remez.m文件可用来设计切比雪夫最佳一致逼近FIR滤波器,同时该文件还可以用来设计希尔伯特变换器和差分器,其调用格式是:
①b=remez(N,F,A);
②b=remez(N,F,A,W);
③b=remez(N,F,A,W,’hilbert’);
④b=remez(N,F,A,W,’differentiator’);
式中N是给定的滤波器的阶次,b是设计的滤波器的系数,其长度为N+1;F是频率向量,A是对应F的各频段上的理想幅频响应,W是各频段上的加权向量。
F的赋值范围是0~1,1对应抽样频率的一半。
需要指出的是,若b的长度为偶数,设计高通和带阻滤波器时有可能出现错误,因此最好保证b的长度为奇数,即N应为偶数。
例:
设计一低通滤波器,要求通带边缘频率Wp=0.6π,阻带边缘频率Ws=0.7π。
实
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数字信号 处理