数字信号处理实验指导书(3)

2018-11-22 18:54

如果所希望的滤波器的理想频率响应函数为Hd(ejw),则其对应的单位脉冲响应为

hd(n)?12?????Hd(ejw)ejwndw

窗函数设计法的基本原理是用有限长单位脉冲响应h(n)逼近hd(n)。由于hd(n)往往是无限长序列,且是非因果的,所以用窗函数w(n)将hd(n)截断,并进行加权处理,得到:

h(n)?hd(n)?w(n)

h(n)就作为实际设计的FIR数字滤波器的单位脉冲响应序列,其频率响应函数H(ejw)为

H(e)??h(n)e?jwn

jwn?0N?1式中,N为所选窗函数w(n)的长度。

这种对理想单位取样响应的加窗处理对滤波器的频率响应会产生以下三点影响:

(1)使理想特性不连续的边沿加宽,形成一过渡带,过渡带的宽度取决于窗函数频谱的主瓣宽度。 (2)在过渡带两旁产生肩峰和余振,它们取决于窗函数频谱的旁瓣;旁瓣越多,余振也越多;旁瓣相对值越大,肩峰则越强。

(3)增加截断长度N,只能缩小窗函数频谱的主瓣宽度而不能改变旁瓣的相对值;旁瓣与主瓣的相对关系只决定于窗函数的形状。因此增加N,只能相对应减小过渡带宽。而不能改变肩峰值。肩峰值的大小直接决定通带内的平稳和阻带的衰减,对滤波器性能有很大关系。例如矩形窗的情况下,肩峰达8.95%,致使阻带最小衰减只有21分贝,这在工程上往往是不够的。怎样才能改善阻带的衰减特性呢?只能从改善窗函数的形状上找出路,所以希望的窗函数频谱中应该减少旁瓣,使能量集中在主瓣,这样可以减少肩峰和余振,提高阻带衰减。而且要求主瓣宽度尽量窄,以获得较陡的过渡带,然而这两个要求总不能同时兼得,往往需要用增加主瓣宽度带换取较大的阻带衰急,于是提出了海明窗、汉宁窗、布莱克曼窗、凯塞窗、切比雪夫窗等窗函数。

所以,用窗函数法设计的滤波器性能取决于窗函数w(n)的类型及窗口长度N的取值。设计过程中,要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度。设待求滤波器的过渡带用?w表示,它近似等于窗函数主瓣宽度。因过渡带?w近似与窗口长度成反比,N?A/?w,A决定于窗口形式。例如,矩形窗A=4π,海明窗A=8π等。按照过渡带及阻带衰减情况,选择窗函数形式。原则是在保证阻带衰减满足要求的情况下, 尽量选择主瓣窄的窗函数。

这样选定窗函数类型和窗口长度N后,求出单位脉冲响应h(n)?hd(n)?w(n),再求出H(ejw)。H(ejw)是否满足要求,要进行验算。一般在h(n)的尾部加零使长度满足2的整数次幂,以便用FFT计算H(ejw)。如果要观察细节,补零点数增多即可。如果H(ejw)不满足要求,则要重新选择窗函数类型和长度N,再次验算,直至满足要求。

如果要求线性相位特性,则h(n)还必须满足:

h(n)??h(N?1?n)

根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类。例如,要设计线性相位低通特性,可选择h(n)?h(N?1?n)一类, 而不能选

h(n)??h(N?1?n)一类。 4、实验步骤及内容

(1)根据下列技术指标,设计一个线性相位的FIR数字低通滤波器。通带截止频率wp?0.2?,通带允许波动Ap?0.25dB;阻带截止频率wS?0.4?,阻带衰减AS?50dB。 ?

(2)调用fir1()函数得到所设计的低通滤波器的单位脉冲响应,调用fft()函数进行频响验证。打印输出各部分结果。

(3)编程验证窗长和窗形状对实际滤波器性能的影响。如要求用窗函数法设计一个线性相位FIR数字低通滤波器,用理想低通滤波器作为逼近滤波器,截止频率wc?

11

?4rad,用四种窗函数(矩形窗,汉宁窗(升余弦

窗),哈明窗(改进的升余弦窗),布莱克曼窗)设计该滤波器,选择窗函数的长度N?15,33两种情况。 (4)输入为20Hz正弦和200Hz的正弦的叠加波形,要求用窗函数法设计一数字低通滤波器滤除200Hz的正弦,使输出中只保留20Hz的正弦波。并绘制出滤波前和滤波后的波形。 5、实验用MATLAB函数介绍

fir1(); fft(); freqz(); boxcar(); hamming(); hanning(); blackman(); sin();figure(); plot(); stem(); abs();angle(); title(); xlabel(); ylabel(); text(); hold on; axis(); grid on; subplot();等

6、思考题

如果要求用窗函数法设计带通滤波器, 且给定上、 下边带截止频率为ω1和ω2,试求理想带通的单位脉冲响应hd(n)。

7、实验报告要求

(1)简述实验目的及实验原理。

(2)编程实现各实验内容,列出实验清单及说明。 (3)总结窗函数法设计FIR数字滤波器的步骤。 (4)简要回答思考题。

附 录

一、常用的MATLAB函数介绍

1.plot

功能:线型绘图函数。 格式:plot(v);plot(x,y) 说明:

plot(v),v是长度为N的数值向量。作用是在坐标系中顺序地用直线连接顶点,生成一条折(曲)线。当向量元素充分多时,即可生成一条光滑的曲线。在实验中,若FFT点数足够多时,用plot打印的幅频特性就很接近|X(ejw)|连续曲线。

plot(x,y),x和y都是长度为n的向量。作用在坐标系中生成顺序连接顶点的折(曲)线。这种调用可被用来生成参数方程的图形。 2.stem

功能:绘制离散序列图。

格式:stem(y); stem(x,y) ;stem(…,’线端符号’) ;stem(…,’线型’); stem(…,‘线型’,‘线端符号’) 说明:

stem(y)和stem(x,y)分别与plot(v)和plot(x,y)的绘图规则相同,只是stem绘制的是离散序列图。实验中用于绘制时域序列x(n)的波形图和序列的离散傅里叶变换X(k)的幅度图。

3.subplot

功能:多坐标设置与定位当前坐标系。 格式:subplot(m,n,k) 说明:

subplot(m,n,k) 将图形窗口分成m行n列的m×n块子区域,按行从上到下,从左到右的顺序,在第i块子区定义一个坐标系,使其成为当前坐标系,随后的绘图函数将在该坐标系输出图形。另外,同一个图形窗口的坐标系可以重叠,这样可以产生前面的坐标系遮住后面坐标系的各种图形效果。 4.figure

功能:创建新的图形窗口(用于输出图形的窗口)。

12

格式:figure;figure(h) 说明:

figure 函数创建一个新的图形窗口,并成为当前图形窗口,所创建的图形窗口的序号(句柄值)是按同一MATLAB程序中创建的顺序号。使用figure(h)函数,该方法常用在程序设计中,用于控制将各种波形图输出到相应的图形窗口中。打印输出或存储时,一个图形窗口打印一张图纸或存储一个图形文件。 5.axis

功能:设置图形的坐标范围。

格式:axis([xmin xmax ymin ymax]) 6.grid on

功能:画出图形的分格线。 7.title

功能: 书写图名。 格式:title(‘s’) 8.xlabel

功能: 横坐标名称。 格式:xlabel(‘s’) 9.ylabel

功能: 纵坐标名称。 格式:ylabel(‘s’) 10.text

功能:在图面(x,y)位置处书写字符注释。 格式:text(x,y,‘s’) 11.hold on

功能:通过该命令,保持当前图形不变化,准备在当前的图形窗口上绘制新的曲线。 12.abs

功能:求绝对值(模值)。 格式:y=abs(x) 说明:

y=abs(x)用于计算x的绝对值,当x为复数时,得到的是复数的模值。当x为字符串时,abs(x)得到字符串的各个字符的ASCII码,例如,x=’123’,则,abs(x)得到:49 50 51。 13.angle

功能:求相角。 格式:ψ=angle(h) 说明:

ψ=angle(h)用于求复矢量或复矩阵的相角(以弧度为单位),相角介于??和??之间。 14.conv

功能:求卷积。 格式:c=conv(a,b)

说明:conv(a,b)用于求矢量a和b的卷积,即

c(n)??a(k?1)b(n?k),n?1,2,?

k?0N?1式中N为矢量a和b的最大长度。例如,当a=[1 2 3],b=[4 5 6]时,则

c=conv(a,b) c=

4 13 28 27 18

此函数可直接用于求两个有限长序列的卷积。设x(n)和h(n)的长度分别为M和N,则计算二者卷积的MATLAB语句如下:

y=conv(x,h) y的长度为M+N-1。 15.filter

13

功能:利用IIR滤波器或FIR滤波器对数据进行滤波。

格式:y=filter(b,a,x); [y,zf]=filter(b,a,x); y=filter(b,a,x,zi) 说明:

filter利用数字滤波器对数据进行滤波,其实现采用直接Ⅱ型结构,因而适用于IIR和FIR两种滤波器。滤波器的系统函数为

b0?b1z?1???bMz?M H(z)??1?N1?a1z???aNz即滤波器系数a=[a0 a1 a2 …aN],b=[b0 b1 … bM]。

这里的标准形式为a0=1,如果输入矢量a时,a0≠1,则MATLAB将自动进行归一化系数的操作;如果a0=0,则给出出错信息。

y=filter(b,a,x)利用给定系数矢量a和b对x中的数据进行滤波,结果放入y矢量中。 y=filter(b,a,x,zi)可在zi中指定x的初始状态。

[y,zf]=filter(b,a,x)除得到矢量y外,还得到x的最终状态矢量zf。 16.fftfilt

功能:利用重叠相加法实现短序列b(n)和长序列x(n)的块之间的卷积。 格式:y=fftfilt(b,x) 17.freqz

功能:数字滤波器的频率响应。 格式:[H,w]=freqz(b,a,N);[H,f]=freqz(b,a,N,Fs);

H=freqz(b,a,w); H=freqz(b,a,f,Fs); freqz(b,a) 说明:

freqz用于计算数字滤波器H(z)的频率响应函数H(ejw)。

[H,w]=freqz(b,a,N)可得到数字滤波器的N点频率响应值,这N个点均匀地分布在[0,?]上,并将这N个频点的频率记录在w中,相应的频响值记录在H中。要求N为大于零的整数,最好为2的整数次幂,以便采用FFT计算,以提高速度。

[H,f]=freqz(b,a,N,Fs)用于对H(e)在[0,Fs/2]上等间隔采样N点,采样点频率及相应频响值分别记录在f和H中。由用户指定Fs (以Hz为单位)值。

H=freqz(b,a,w)用于对H(ejw)在[0, 2?]上进行采样,采样频率点由矢量w指定。 H=freqz(b,a,f,Fs) 用于对H(ejw)在[0, Fs]上采样,采样频率点由矢量f给定。

freqz(b,a)用于在当前窗口中绘制出幅频和相频特性曲线。 18.impz

功能:计算H(z)相应的单位脉冲响应h(n)。 格式:[h,t]=impz(b,a); [h,t]=impz(b,a,z); 19.fft

功能:一维快速傅里叶变换(FFT) 格式:y=fft(x); y=fft(x,N) 说明:

fft函数用于计算矢量或矩阵的离散傅里叶变换。

y=fft(x)利用FFT算法计算矢量x的离散傅里叶变换,当x为矩阵时,y为矩阵x每一列的FFT。当x的长度为2的整数次幂时,fft采用基2FFT算法,否则采用稍慢的混合基算法。

y=fft(x,N)采用N点FFT。当x长度小于N时,fft函数自动在x尾部补零,以构成N点数据;当x的长度大于N时,fft截取x的前面N点数据进行FFT。 20.ifft

功能:一维逆快速傅里叶变换(IFFT) 格式:y=ifft(x); y=ifft(x,N)

jw二、MATLAB在信号处理中的应用举例 1、线性卷积与循环卷积的计算

对于无限长序列不能用MATLAB直接计算线性卷积,在MATLAB内部只提供了一个conv函数计算两个有限长序列的线性卷积。对于循环卷积MATLAB内部没有提供现成的函数,我们可以按照定义式直接编程计算。

例1: 已知两序列:

14

?0.8n0?n?11x(n)??其它?00?n?5?1h(n)??其它?0

求它们的线性卷积yl(n)?h(n)?x(n)和N点的循环卷积yc(n)?h(n)?x(n),并研究两者之间的关系。

MATLAB实现程序: 循环卷积的函数

function yc=circonv(x1,x2,N)

%realize circular convolution use direct method %y=circonv(x1,x2,N) %y:output sequences %x1,x2:input sequences %N:circulation length if length(x1)>N

error('N must not be less than length of x1'); end

if length(x2)>N

error('N must not be less than length of x2'); end

%以上语句判断两个序列的长度是否小于N x1=[x1,zeros(1,N-length(x1))];

%填充序列x1(n)使其长度为N1+N2-1(序列%h(n)的长度为N1,序列x(n)的长度为N2) x2=[x2,zeros(1,N-length(x2))]; %填充序列x2(n)使其长度为N1+N2-1 n=[0:1:N-1];

x2=x2(mod(-n,N)+1); %生成序列x2((-n))N H=zeros(N,N); for n=1:1:N

H(n,:)=cirshiftd(x2,n-1,N); %该矩阵的k行为x2((k-1-n))N end

yc=x1*H'; %计算循环卷积

function y=cirshiftd(x,m,N)

%directly realize circular shift for sequence x %y=cirshiftd(x,m,N);

%x:input sequence whose length is less than N %m:how much to shift %N:circular length

%y:output shifted sequence if length(x)>N

error('length of x must be less than N'); end

x=[x,zeros(1,N-length(x))]; n=[0:1:N-1];

y=x(mod(n-m,N)+1); 研究两者之间的关系 clear all; n=[0:1:11]; m=[0:1:5]; N1=length(n); N2=length(m);

15


数字信号处理实验指导书(3).doc 将本文的Word文档下载到电脑 下载失败或者文档不完整,请联系客服人员解决!

下一篇:2016年吉林省高低压电器装配工技能考试试题

相关阅读
本类排行
× 注册会员免费下载(下载后可以自由复制和排版)

马上注册会员

注:下载文档有可能“只有目录或者内容不全”等情况,请下载之前注意辨别,如果您已付费且无法下载或内容有问题,请联系我们协助你处理。
微信: QQ: