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

更新时间:2023-09-18 12:51:01 阅读量: 幼儿教育 文档下载

说明:文章内容仅供预览,部分内容可能不全。下载后的文档,内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的,是否完整无缺。

电 工 电 子 实 验 中 心

实 验 指 导 书

数字信号处理 实验教程

二○○九 年 三 月

高等学校电工电子实验系列

数字信号处理实验教程

主编 石海霞 周玉荣

攀枝花学院电气信息工程学院

电工电子实验中心

内容简介

内 容 简 介

数字信号处理是一门理论与实践紧密联系的课程,适当的上机实验有助于深入理解和巩固验证基本理论知识,了解并体会数字信号处理的CAD手段和方法,锻炼初学者用计算机和MATLAB语言及其工具箱函数解决数字信号处理算法的仿真和滤波器设计问题的能力。

本实验指导书结合数字信号处理的基本理论和基本内容设计了八个上机实验,每个实验对应一个主题内容,包括常见离散信号的MATLAB产生和图形显示、离散时间系统的时域分析、离散时间信号的DTFT、离散时间信号的Z变换、离散傅立叶变换DFT、快速傅立叶变换FFT及其应用、基于MATLAB的IIR和FIR数字滤波器设计等。此外,在附录中,还简单介绍了MATLAB的基本用法。每个实验中,均给出了实验方法和步骤,还有部分的MATLAB程序,通过实验可以使学生掌握数字信号处理的基本原理和方法。

目录

目 录

绪 论 ............................................................................................................................ 1 实验一 常见离散信号的MATLAB产生和图形显示.................................................... 2 实验二 离散时间系统的时域分析 ............................................................................. 6 实验三 离散时间信号的DTFT .................................................................................... 9 实验四 离散时间信号的Z变换 ............................................................................... 14 实验五 离散傅立叶变换DFT .................................................................................... 18 实验六 快速傅立叶变换FFT及其应用 ................................................................... 24 实验七 基于MATLAB的IIR数字滤波器设计(设计实验) ...................................... 30 实验八 基于MATLAB的FIR数字滤波器设计(设计实验) ...................................... 33 附 录 .......................................................................................................................... 37 参 考 文 献 .............................................................................................................. 40

绪论

绪 论

随着电子技术迅速地向数字化发展,《数字信号处理》越来越成为广大理工科,特别是IT领域的学生和技术人员的必修内容。

数字信号处理是把信号用数字或符号表示成序列,通过计算机或通用(专用)信号处理设备,用数值计算方法进行各种处理,达到提取有用信息便于应用的目的。数字信号处理的理论和技术一出现就受到人们的极大关注,发展非常迅速。而且随着各种电子技术及计算机技术的飞速发展,数字信号处理的理论和技术还在不断丰富和完善,新的理论和技术层出不穷。目前数字信号处理已广泛地应用在语音、雷达、声纳、地震、图象、通信、控制、生物医学、遥感遥测、地质勘探、航空航天、故障检测、自动化仪表等领域。

数字信号处理是一门理论和实践、原理和应用结合紧密的课程,由于信号处理涉及大量的运算,可以说离开了计算机及相应的软件,就不可能解决任何稍微复杂的实际应用问题。Matlab是1984年美国Math Works公司的产品,MATLAB语言具备高效、可视化及推理能力强等特点,它的推出得到了各个领域专家学者的广泛关注,其强大的扩展功能为各个领域的应用提供了基础,是目前工程界流行最广的科学计算语言。早在20世纪90年代中期,MATLAB就己成为国际公认的信号处理的标准软件和开发平台。从1996年后,美国新出版的信号处理教材就没有一本是不用MATLAB的。

本实验指导书结合数字信号处理的基本理论和基本内容,用科学计算语言MATLAB实现数字信号处理的方法和实践, 通过实验用所学理论来分析解释程序的运行结果,进一步验证、理解和巩固学到的理论知识,从而达到掌握数字信号处理的基本原理和方法的目的。

1

实验一 常见离散信号的MATLAB产生和图形显示

实验一 常见离散信号的MATLAB产生和图形显示

一、实验目的

1. 学会用MATLAB在时域中产生一些基本的离散时间信号。 2. 了解信号的各种运算。

二、实验原理

(一)、序列的产生

由于MATLAB数值计算的特点,用它来分析离散时间信号与系统是很方便的。离散信号是数字信号处理的最基础的内容,由于内存有限,MATLAB无法表示无限序列。在MATLAB中,可以用一个列向量来表示一个有限长度的序列,但是这种表示方法没有包含采样位置的信息,要完全表示x(n),要用x 和n 两个向量,例如

x(n)={2,1,0,2,3,-1,2,3} ↑ 在MATLAB中表示为

n=[-4,-3,-2,-1,0,1,2,3]; x=[2,1,0,2,3,-1,2,3];

当序列从n=0开始,则不需要采样位置信息,这时可以只用x 来表示

1. 单位抽样序列

?(n)??

?0?1n?0n?0

在MATLAB中可以利用zeros()函数实现。

x?zeros(1,N);x(1)?1;

如果?(n)在时间轴上延迟了k个单位,得到?(n?k)即:

?(n?k)??

?0?1n?kn?0

2.单位阶跃序列

2

实验一 常见离散信号的MATLAB产生和图形显示

n?0?1 u(n)?

n?0?0在MATLAB中可以利用ones()函数实现。

x?ones(1,N)

3.正弦序列

x(n)?Asin(2?fn/Fs??)

在MATLAB中

n?0:N?1x?A*sin(2*pi*f*n/Fs?fai)

4.复正弦序列

x(n)?ej?n

在MATLAB中

n?0:N?1x?exp(j*w*n)

5.指数序列

x(n)?a

n在MATLAB中

n?0:N?1x?a.^n

(二)、简单运算 1.信号加

x(n)?x1(n)?x2(n) 在MATLAB中

x=x1+x2

注意:x1和x2应该具有相同的长度,位置对应,才能相加。否则,需要先通过Zeros函数左右补零后再相加。

2.信号延迟

3

实验一 常见离散信号的MATLAB产生和图形显示

给定离散信号x(n),若信号y(n)为 y(n)= x(n-k)

那么y(n)就是信号x(n)在时间轴上右移k个采样周期后得到的新的序列。 在MATLAB中

y(n)= x(n-k) 3.信号乘

x(n)?x1(n)?x2(n)

在MATLAB中

x?x1.?x2

这是信号的点乘运算,所以同样需要信号加需要的x1和x2二者的长度要相等这一前提条件。

4.信号变化幅度

y(n)= k×x(n) 在MATLAB中

y?k?x

5.信号翻转

y(n)= x(-n)

在MATLAB中

y=fliplr(x)

三、实验内容与步骤

1. 产生一个单位样本序列x1(n),起点为ns= -10, 终点为nf=20, 在n0=0时有一单位脉冲并显示它。修改程序,以产生带有延时11个样本的延迟单位样本序列x2(n)= x1(n-11),并显示它。

2.已知 c= -(1/12)+(pi/6)*i;产生一个复数值的指数序列x2(n)=2*exp(c*n),起点为ns= 0, 终点为nf=40;并显示它。

3.产生一个正弦序列x3(n) =1.5*cos(2*pi*f*n); 起点为ns= 0, 终点为nf=40;并显示它。

4.复杂信号的产生:复杂的信号可以通过在简单信号上执行基本的运算来产生 试产生一个振幅调制信号

y(n)?(1?m?cos(2?fLn))?cos(2?fHn)?(1?0.4?cos(2??0.01n))?cos(2??0.1n)

4

实验一 常见离散信号的MATLAB产生和图形显示

n=0:100

四、实验仪器设备

计算机,MATLAB软件

五、实验注意事项

预先阅读附录(MATLAB基础介绍);

六、思考题

1. 讨论复指数序列x2(n)的哪个参数控制该序列的增长或衰减率?哪个参数控制该序列的振幅?

2. 讨论正弦序列x3(n)的哪个参数控制该序列的相位?哪个参数控制该序列的振幅?

3. 讨论算术运算符*和. *之间的区别是什么?

5

实验二 离散时间系统的时域分析

实验二 离散时间系统的时域分析

一、实验目的

1. 运用MATLAB仿真一些简单的离散时间系统,并研究它们的时域特性。 2. 运用MATLAB中的卷积运算计算系统的输出序列,加深对离散系统的差分方程、冲激响应和卷积分析方法的理解。

二、实验原理

离散时间系统其输入、输出关系可用以下差分方程描述:

N?k?0dky[n?k]?M?k?0pkx[n?k]

当输入信号为冲激信号时,系统的输出记为系统单位冲激响应

?[n]?h[n],则系统响应为如下的卷积计算式:

?y[n]?x[n]?h[n]??x[m]h[n?m]

m??? 当h[n]是有限长度的(n:[0,M])时,称系统为FIR系统;反之,称系统为IIR系统。在MATLAB中,可以用函数y=Filter(p,d,x) 求解差分方程,也可以用函数 y=Conv(x,h)计算卷积。 例2.1

clf; n=0:40; a=1;b=2; x1= 0.1*n; x2=sin(2*pi*n); x=a*x1+b*x2; num=[1, 0.5,3];

6

实验二 离散时间系统的时域分析

den=[2 -3 0.1];

ic=[0 0]; %设置零初始条件

y1=filter(num,den,x1,ic); %计算输入为x1(n)时的输出y1(n) y2=filter(num,den,x2,ic); %计算输入为x2(n)时的输出y2(n) y=filter(num,den,x,ic); %计算输入为x (n)时的输出y(n) yt= a*y1+b*y2; %画出输出信号 subplot(2,1,1) stem(n,y);

ylabel(‘振幅’);

title(‘加权输入a*x1+b*x2的输出’); subplot(2,1,2) stem(n,yt); ylabel(‘振幅’);

title(‘加权输出a*y1+b*y2’);

(一)、线性和非线性系统

对线性离散时间系统,若y1(n)和y2(n)分别是输入序列x1(n)和x2(n)的响应,则输入x(n)?ax1(n)?bx2(n)的输出响应为y(n)?ay1(n)?by2(n),即符合叠加性,其中对任意常量a和b以及任意输入x1(n)和x2(n)都成立,否则为非线性系统。 (二)、时不变系统和时变系统

对离散时不变系统,若y1(n)是x1(n)的响应,则输入x(n)=x1(n-n0)的输出响应为y(n)=y1(n-n0),式中n0是任意整数。该输入输出关系,对任意输入序列及其相应的输出成立,若对至少一个输入序列及其相应的输出序列不成立,则系统称之为时变的。 (三)、线性卷积

假设待卷积的两个序列为有限长序列,卷积运算符在MATLAB中可 命令conv实现。例如,可以把系统的冲激响应与给定的有限长输入序列进行卷积,得到有限长冲激响应系统的输出序列。下面的MATLAB程序实现了该方法。 例2.2

clf;

7

实验二 离散时间系统的时域分析

h=[3 2 1 -2 1 0 -4 0 3];%冲激 x=[1 -2 3 -4 3 2 1 ]; %输入序列 y=conv(h,x); n=0:14; stem(n,y);

xlabel(‘时间序号n’);ylabel(‘振幅’); title(‘用卷积得到的输出’);grid;

三、实验内容与步骤

1. 假定一因果系统为

y(n)-0.4y(n-1)+0.75y(n-2)=2.2403x(n)+2.4908x(n-1)+2.2403x(n-2)

用MATLAB程序仿真该系统,输入三个不同的输入序列:

??0.1n),x2(n)?cos(2??0.4n),x?2x1(n)?3x2(n) x1(n)?cos2(计算并并显示相应的输出y1(n), y2(n)和y(n)。

2. 用MATLAB程序仿真步骤1给出的系统,对两个不同的输入序列x(n)和x(n-10),计算并显示相应的输出序列y3(n)和y4(n)。

3.用MATLAB程序仿真计算下列两个有限长序列的卷积和并显示图形。

x1(n)??(n)?3?(n?1)?2?(n?2) x2(n)?u(n)?u(n?3)

四、实验仪器设备

计算机,MATLAB软件

五、实验要求

给出理论计算结果和程序计算结果并讨论。

六、思考题

1. 讨论实验程序1,对由加权输入得到的y(n)与在相同权系数下输出y1(n)和y2(n)相加得到的yt(n)进行比较,这两个序列是否相等?该系统是线性系统吗? 2. 讨论实验程序2,比较输出序列y3(n)和y4(n),这两个序列之间有什么关系?该系统是时不变系统吗?

3. 讨论实验程序3的理论计算结果和程序计算结果是否一致。

8

实验三 离散时间信号的DTFT

实验三 离散时间信号的DTFT

一、实验目的

1. 运用MATLAB计算离散时间系统的频率响应。 2. 运用MATLAB验证离散时间傅立叶变换的性质。

二、实验原理

(一)、计算离散时间系统的DTFT

NNk 已知一个离散时间系统

?ak?0y(n?k)??bk?0kx(n?k),可以用MATLAB函数

frequz非常方便地在给定的L个离散频率点???l处进行计算。由于H(ej?)是ω的连续函数,需要尽可能大地选取L的值(因为严格说,在MATLAB中不使用symbolic工具箱是不能分析模拟信号的,但是当采样时间间隔充分小的时候,可产生平滑的图形),以使得命令plot产生的图形和真实离散时间傅立叶变换的图形尽可能一致。在MATLAB中,freqz计算出序列{b0,b1,?,bM}和{a0,a1,?,aN}的

L点离散傅立叶变换,然后对其离散傅立叶变换值相除得到

j?lH(e),l?1,2,?,L。为了更加方便快速地运算,应将L的值选为2的幂,如256

或者512。

例3.1 运用MATLAB画出以下系统的频率响应。 y(n)-0.6y(n-1)=2x(n)+x(n-1) 程序:

clf;

w=-4*pi:8*pi/511:4*pi; num=[2 1];den=[1 -0.6]; h=freqz(num,den,w); subplot(2,1,1) plot(w/pi,real(h));grid title(‘H(e^{j\\omega}的实部’)) xlabel(‘\\omega/ \\pi’);

9

实验三 离散时间信号的DTFT

ylabel(‘振幅’);

subplot(2,1,1)

plot(w/pi,imag(h));grid title(‘H(e^{j\\omega}的虚部’)) xlabel(‘\\omega/ \\pi’); ylabel(‘振幅’);

(二)、离散时间傅立叶变换DTFT的性质。 1.时移与频移

设 X(ej?)?FT[x(n)], 那么

FT[x(n?n0)]?e?j?n0X(ej?) (2.2.6)

FT[ej?0nx(n)]?X(ej(???0)) (2.2.7)

2.时域卷积定理

如果 y(n)?x(n)?h(n), 那么

Y(ej?

)?X(ej?)?H(ej?)

三、实验内容与步骤

1. 已知因果线性时不变离散时间系统

y(n)-0.4y(n-1)+0.75y(n-2)=2.2403x(n)+2.4908x(n-1)+2.2403x(n-2) 运用MATLAB画出该系统的频率响应。

2.运行下面程序并显示它,验证离散时间傅立叶变换DTFT的时移性。

clf;

w=-pi:2*pi/255: pi;wo=0.4*pi;D=10; num=[1 2 3 4 5 6 7 8 9]; h1=freqz(num,1,w);

h2=freqz([zeros(1,D) num],1,w); subplot(2,2,1)

10

实验三 离散时间信号的DTFT

plot(w/pi,abs(h1));grid title(‘原序列的幅度谱’) subplot(2,2,2)

plot(w/pi,abs(h2));grid title(‘时移后序列的幅度谱’) subplot(2,2,3)

plot(w/pi,angle (h1));grid title(‘原序列的相位谱’) subplot(2,2,4)

plot(w/pi, angle (h2));grid title(‘时移后序列的相位谱’)

3. 运行下面程序并显示它,验证离散时间傅立叶变换DTFT的频移性。

clf;

w=-pi:2*pi/255: pi;wo=0.4*pi;D=10;

num1=[1 3 5 7 9 11 13 15 17];L=length(num1); h1=freqz(num1,1,w);n=0:L-1; num2=exp(wo*i*n).*num1; h2=freqz(num2,1,w); subplot(2,2,1)

plot(w/pi,abs(h1));grid title(‘原序列的幅度谱’) subplot(2,2,2)

plot(w/pi,abs(h2));grid title(‘频移后序列的幅度谱’) subplot(2,2,3)

plot(w/pi,angle (h1));grid title(‘原序列的相位谱’) subplot(2,2,4)

plot(w/pi, angle (h2));grid title(‘频移后序列的相位谱’)

11

实验三 离散时间信号的DTFT

4.运行下面程序并显示它,验证离散时间傅立叶变换时域卷积性质。

clf;

w=-pi:2*pi/255: pi; x1=[1 3 5 7 9 11 13 15 17]; x2=[1 -2 3 -2 1]; y=conv(x1,x2); h1=freqz(x1,1,w); h2=freqz(x2,1,w); hp=hi.*h2; h3=freqz(y,1,w); subplot(2,2,1)

plot(w/pi,abs (hp));grid title(‘幅度谱的乘积’) subplot(2,2,2)

plot(w/pi,abs (h3));grid title(‘卷积后序列的幅度谱’) subplot(2,2,3)

plot(w/pi,angle (hp));grid title(‘相位谱的和’) subplot(2,2,4)

plot(w/pi,angle (h3));grid title(‘卷积后序列的相位谱’)

四、实验仪器设备

计算机,MATLAB软件

五、实验注意事项

课前预先阅读并理解实验程序;

六、思考题

1.讨论实验程序1中的离散时间系统的频率响应是离散的还是连续的,是否是周期的?周期为多少?y与

2.讨论实验程序2中h1和h2的关系是什么?哪个参数控制时移量?

12

实验三 离散时间信号的DTFT

3. 讨论实验程序3中h1和h2的关系是什么?哪个参数控制频移量?

4. 讨论实验程序4中y与x1和x2的关系是什么?h1和h2与x1和x2的关系是什么?h1和hp相等吗?

13

实验四 离散时间信号的Z变换

实验四 离散时间信号的Z变换

一、实验目的

1. 运用MATLAB理解Z变换及其绘制H(z)的零极点图。 2. 运用MATLAB计算逆Z变换。

二、实验原理

(一)、MATLAB在ZT中的应用。

线性时不变离散时间系统的冲激响应h(n)的z变换是其系统函数H(z), 在MATLAB中可以利用性质求解Z变换,例如可以利用线性卷积求的Z变换。若H(z)的收敛域包含单位圆,即系统为稳定系统,即系统在单位圆上z?ej?处计算的是系统的频率响应。

(二)、逆Z变换

Z变换对于分析和表示离散线性时不变系统具有重要作用。但是在MATLAB中不能直接计算Z变换,但是对于一些序列可以进行逆Z变换。

已知序列的Z变换及其收敛域, 求序列称为逆Z变换。 序列的Z变换及共逆Z变换表示如下:

?X(z)?x(n)??n???x(n)z?n,Rx??z?Rx?n?1

12?

通常,直接计算逆Z变换的方法有三种:围线积分法、长除法和部分分式展开法。在实际中,直接计算围线积分比较困难,往往不直接计算围线积分。由于序列的Z变换常为有理函数,因此采用部分分式展开法比较切合实际,它是将留数定律和常用序列的Z变换相结合的一种方法。

设x(n)的Z变换X(z)是有理函数,分母多项式是N阶,分子多项式是M阶,将X(z)展成一些简单的常用的部分分式之和,通过常用序列的Z变换求得各部分的逆变换,再相加即得到原序列x(n)。在MATLAB中提供了函数residuez来实现上述过程,调用格式如下:

14

?j?cX(z)zdz,c?(Rx?,Rx?)

实验四 离散时间信号的Z变换

[R,P,K]= residuez(B,A)

其中B、A分别是有理函数分子多项式的系数和分母多项式的系数,输出R是留数列向量,P是极点列向量。如果分子多项式的阶数大于分母多项式的阶数,则K返回为常数项的系数。 例4.1 计算 X(z)?的Z反变换。

由于分母多项式为:1?1.5z?1?0.5z?2,则 MATLAB实现: clear b=1; a=[1,-1.5,0.5]; [R,P,K]=residuez(b,a) 输出R =

2 -1 P =

1.0000 0.5000 K = []

因此:得到X(z)的部分分式展开为X(z)?得x(n)?(2?0.5)u(n)。

三、实验内容与步骤

1. 运行下面程序,利用线性卷积求Z变换。 设X1(z)?z?2?3z?1n1(1?z?1)(1?0.5z?1),z?1

21?z?1??11?0.5z?1,根据常用序列的z变换可

,X2(z)?2z?4z?3?5z2?1,求X3(z)?X1(z)?X2(z)

由变换定义可知:x1(n)?{1,2,3},n?{?1,0,1} x2(n)?{2,4,3,5},n?{?2,?1,0,1}

15

实验四 离散时间信号的Z变换

通过求x3(n)?x1(n)?x2(n),再求其z变换得到X3(z)?X1(z)?X2(z)。 MATLAB程序: x1=[]1,2,3]; x2=[2,4,3,5]; n1=-1:1; n2=-2:1; x3=conv(x1,x2) nb3=n1(1)+n2(1);

nc3=n1(length(x1))+n2(length(x2)); n3=[nb3:nc3]

2.已知两个线性时不变的因果系统,系统函数分别为

H1(z)?1?z?N ,

H2(z)?1?zN?N?N1?az

分别令N=8,a=0.8,运行下面程序计算并显示示这两个系统的零、极点图及幅频特性。 程序:

b=[1,0,0,0,0,0,0,0,-1]; %H1(z)和H2(z)的分子多项式系数向量 a0=1; %H1(z)分母多项式系数向量

a1=[1,0,0,0,0,0,0,0,-(0.8)^8]; % H2(z)的分母多项式系数向量 [H, w]=freqz(b,a0); [H1, w1]=freqz(b,a1); subplot(2,2,1);zplane(b,a0);

xlabel(‘实部’);ylabel(‘虚部’); title(‘H1 (z)系统的零极点图’); subplot(2,2,2);zplane(b,a1);

16

附录A MATLAB简介

附 录

MATLAB是一种用于科学计算且功能很强的高级程序设计语言。它易于学会,主要用来解决复杂的工程问题。

MATLAB由内建于解释器中的函数或以M文件存在的函数组成,他是包含有实现某种算法的程序语句序列。MATLAB在电脑屏幕上以三种类型的窗口工作。它们是命令窗口、图形窗口和编辑器窗口。命令窗口的标题为Command,图形窗口的标题为Figure No.1,而编辑器窗口的标题或是一个已经打开的现存M文件的名字,或者为Unititled,代表创建一个新的M文件。命令窗口显示符>>,表明它准备执行MATLAB命令。大部分输入命令的结果在命令窗口显示,同时这个窗口也可用来运行小程序和保存M文件。由画图命令产生的所有图形在图形窗口出现。

1. 数字和数据表示

MATLAB使用常规的十进制符号表示数字。负数可用一个位于最前面的负号来表示。MATLAB中可表示的数字的近似范围从10-308到10308。很大或很小的数字可以以指数形式表示。下面是一些有效数字表示的典型例子: 1234.56789 l 23456.789E·2 1.23456789e3 —123456.789 注意,在指数前不应有空格。

MATLAB中的数据以矩阵表示且不需要预先指定维数,元素可以是实数或复数。矩阵元素可以用两种不同的方式输入。既可以在用分号分开的单行键入,又可以在不同的行键人。例如,3x 4矩阵A ?1? A?2???9341056117??8 ?12??可以以

A=[1 3 5 7;2 4 6 8;9 10 11 12];或以

A=[1 3 5 7

2 4 6 8 9 10 11 12];

的形式输入。

37

附录A MATLAB简介

2. 算术运算

如下所示.MATLAB中有两种不同类型的算术运算用来对已存储的数据进 行操作,其中x和y表示两个不同的矩阵。若x和y具有相同的维数,x和y的相加由表达式x+y实现。相加运算+也可用来把一个标量加到一个矩阵中。同样地,从x减去y由表达式x—y实现。减法运算—也可用来将一个标量从一个矩阵中减去。

若x的列数与Y的行数相同,矩阵相乘x﹡y产生x和y的线性代数乘积。乘法运算﹡也可用来将一个标量乘以一个矩阵。若x和y有同样的维数,x.﹡y是一个数组相乘,它形成X和y的逐元素的乘积。

若y是一个方阵.而x是一个具有与y的列数相同的矩阵,则矩阵右除x/y相当于x﹡lnv(y),其中inv(y)表示y的逆。若它们中的一个是标量,右除运算x/y也能实现。若y是一个方阵而x是一个带有与y的行数相同的矩阵,则矩阵左除X\Y相当于lnv(y)﹡x。

若x和y有相同的维数,数组右除用表达式x./y实现,得到第一个矩阵,其第(r,s)个元素由x(r,s)/y(r,s)给出。

若在一个语句中使用乘法运算,在求取这个表达式时遵循通常的优先规则,然而,括号可以改变运算的优先顺序。

3. 关系运算符

在MATLAB中,关系运算<,<=,>,>=,==和-=分别表示比较运算小于、 小于等于、大于、大于等于、等于和不等于。使用这些运算符对两个同样大小的矩阵进行比较时,MATLAB将把对应的元素逐个进行比较,其结果也将以同样大小的矩阵出现。且当关系为真时,矩阵中的该元素被设为1,而当关系为假时,该元素被设为0。在复数值的矩阵情况下,运算符<,<=,>,>=将只用来比较矩阵元素的实数部分.但运算符==和-=可用来比较实数和虚数两部分。 4.逻辑运算符

在MATLAB中,三个逻辑运算—,|和- 分别实现与(AND)、或(OR)及非(NOT)运算。当应用到矩阵时,它们将逐个元素进行运算,用0表示FALSE而用1表示TRUE。 5.流程控制

MATLAB的流程控制命令是break,else,elseif,end,error,for,if,return和while。这些命令允许有条件地执行某些程序命令。命令for用于以给定的次数重复一组程

38

附录A MATLAB简介

序指令。命令if用来有条件地执行一组程序指令,而命令while可以用来以一个无法确定的次数重复程序指令。跟在命令for,while和if后面的指令必须用命令end结束。命令break用于结束一个循环的执行。命令else和elseif需与命令if一起使用,它可在一个循环的执行中使其有条件的停止。命令error用于显示错误消息并中断函数。 6.输出数据格式

MATLAB中的所有算术运算都以双精度进行。然而,在命令窗口可用不同的格式来显示操作的结果。若所有的结果是严格的证书,则它们以无小数点的形式显示。若一个或更多数据元素不是整数,则结果可按不同的精度显示。format short显示五个有效小数位,这是默认的格式。format short e显示带有两个正或负十进制指数的五个有效小数位。format long显示15个有效小数位,而format long e显示15个有效小数位加两个正或负十进制指数。

7.图形

MATLAB能将计算结果以高级的图形方式显示出来。在许多情况下,我们仅涉及二维图形,只在一些特殊情况下用到三维图形。对于二维图形,图形可以以不同的形式表示——坐标轴即可以用线性尺度,又可以用对数尺度,或者一个轴用线性尺度而另一个轴用对数尺度。我们可以在图中加入栅格线,在两个轴上加标记,在图顶部加标题等等。

39

参考文献

参 考 文 献

[1] 程佩青.数字信号处理(第三版)[M].北京:清华大学出版社. 2007.2.

[2] 丁玉美,高西全.数字信号处理[M].西安电子科技大学出版社. 2001.1

[3] (美)米特拉(Mitra, S.k.)著.孙洪,余翔宇译.数字信号处理实验指导书(MATLAB版)[M]. 北京:电子工业出版社,2005.1

[4] 聂祥飞,王海宝,谭泽富.MATLAB程序设计及其在信号处理中的应用 [M]. 成都:西南

交通大学出版社,2005.4

[5] 从玉良,王宏志.数字信号处理原理及其MATLAB实现[M].北京:电子工业出版社,2005.7 [6] 赵红怡,张常年.数字信号处理及其MATLAB实现 [M]. 北京:化学工业出版社,2002.1

40

参考文献

41

实验四 离散时间信号的Z变换

xlabel(‘实部’);ylabel(‘虚部’); title(‘H2(z)系统的零极点图’); subplot(2,2,3);plot(w/pi,abs(H)); title(‘H1 (z)系统幅频响应曲线’); subplot(2,2,4);plot(w/pi,abs(H1)); title(‘H2(z)系统幅频响应曲线’);

3.编写MATLAB程序 计算系统

X(z)?1(1?0.3z?1?1?1?1)(1?0.4z)(1?0.5z)(1?0.8z),z?0.8的z反变换。

四、实验仪器设备

计算机,MATLAB软件

五、实验注意事项

课前预先阅读并理解实验程序;

六、思考题

1.讨论实验程序1中的n3代表什么含义?根据x3和n3的结果写出

X3(z)?X1(z)?X2(z)的结果。

2.讨论实验程序2中极点与系统稳定性的关系?根据程序运行结果判断该系统的稳定性。 3. 根据实验程序3的运行结果写出z反变换x(n)。

17

实验五 离散傅立叶变换DFT

实验五 离散傅立叶变换DFT

一、实验目的

1. 运用MATLAB计算有限长序列的DFT和IDFT。 2. 运用MATLAB验证离散傅立叶变换的性质。 3 .运用MATLAB计算有限长序列的圆周卷积。

二、实验原理

(一)、离散傅立叶变换DFT的定义

一个有限长度的序列x(n)(0≤n

(0???2?)上对X(ej?)均匀采样得到

?X(k)?X(ej?)

??2?k/N??n???x(n)e?j2?kn/N 0?k?N?1

可以看到X(k)也是频域上的有限长序列,长度为N。序列X(k)称为序列x(n)的N点DFT。N称为DFT变换区间长度。 通常表示

WN?e?j2?/N

可将定义式表示为

?X(k)?

?x(n)Wn???kn 0?k?N?1

X(k)的离散傅里叶逆变换(IDFT)为

x(n)?1N?

1.圆周移位

?X(k)Wn????kn 0?n?N?1

(二)、DFT的性质

定义序列x(n)的m单位的圆周移位y(n)为:

y(n)?~x(n?m)RN(n)?x((n?m))NRN(n)18

实验五 离散傅立叶变换DFT

(x((n?m))N即对x(n)以N为周期进行周期延拓的序列~RN(n)表x(n)的m点移位,示对此延拓移位后再取主值序列)

2. 圆周卷积

??X1(k) 0?k?N?1 设 x1(n)??NDFT??X2(k) 0?k?N?1 x2(n)??NDFT??X1(k)X2(k) 0?k?N?1 则 x1(n) x2(n)??NDFT这里 x1(n) x2(n) 表示x1(n)与 x2(n)的N点循环卷积。

N?1x1(n) x2(n)??xm?02(m)[x1((n?m))NRN(n)],n?0,1,?,N?1

3. 共轭对称性

x(n)?xep(n)?xop(n),0?n?N?1

19

实验五 离散傅立叶变换DFT

1?*x(n)?[x(n)?x(N?n)]?ep2??,0?n?N?1 1*?xop(n)?[x(n)?x(N?n)]2?DFT??X(k) x(n)??Nxep(n)????NDFT12[X(k)?X(k)]?Re[X(k)]?Xr(k)

*实际应用中,利用上述对称性质可以减少DFT的运算量,提高运算效率。

三、实验内容与步骤

1. 构造离散傅立叶正、反变换函数的MATLAB程序,其中dft(xn,N)为离散傅立叶正变换,idft(xn,N)为离散傅立叶反变换。 function[Xk]=dft(xn,N) n=[0:1:N-1]; k=n;

WN=exp(-j*2*pi/N); nk=n’*k;

WNnk=WN.^nk; Xk=xn*WNnk;

function[xn]=idft(xn,N) n=[0:1:N-1]; k=n;

WN=exp(-j*2*pi/N); nk=n’*k;

WNnk=WN.^(-nk); xn =(Xk*WNnk)/N;

如果x(n)?sin(n?/8)?sin(n?/4)是一个N=16的有限长序列,利用离散傅立叶变换函数求其16点DFT,并显示其DFT结果。

2. 利用MATLAB程序求有限长序列x(n)=8(0.4)n, 0≤n<20的圆周移位

xm(n)?x[(n?10)]20R20(n),并显示其图形。

程序:

20

实验五 离散傅立叶变换DFT

N=20; m=10; n=0:1:N-1; x=8* (0.4).^n;

n1=mod((n+m),N); %求余

xm=x(n1+1); % 求余后加1是因为MATLAB向量下标从1开始 subplot(2,1,1) stem(n,x); title(‘原序列’); xlabel(‘n’); ylabel(‘x(n)’); subplot(2,1,2) stem(n,xm);

title(‘圆周移位序列’); xlabel(‘n’); ylabel(‘xm(n)’);

3.构造一个圆周卷积函数,专门用来计算两个序列(假设均从0点开始)的循环卷积。

Function[y]=circonv(x1,x2) xn2=[x2(1),fliplr(x2)]; xn2(length(xn2))=[]; C=xn2;R=x2; M=toeplitz(C,R); y=x1* M; 如果

x1(n)??(n)?2?(n?1)?3?(n?2)?4?(n?3),

x2(n)?4?(n)?3?(n?1)?2?(n?2)??(n?3),利用上面的圆周卷积函数求

其4点圆周卷积,并显示其结果。再利用时域圆周卷积定理,调用前面的离散傅立叶正、反变换函数计算它们的4点圆周卷积,并显示其结果。

21

实验五 离散傅立叶变换DFT

4. 调用离散傅立叶正变换函数运行下面MATLAB程序,验证DFT的共轭对

称性。

x=[1 2 4 2 6 32 6 4 2 zeros(1,247)];

x1=[x(1) x(256:-1:2)];

xep=0.5*(x+x1); xop=0.5*(x-x1); XF=dft(x,256); XEPF=dft(xep,256); XOPF=dft(xop,256); clf k=0:255; subplot(2,2,1);

plot(k/128,real(XF));grid; ylabel(‘振幅’);

title(‘原序列DFT的实部’);

subplot(2,2,2);

plot(k/128,imag(XF));grid; ylabel(‘振幅’);

title(‘原序列DFT的虚部’); subplot(2,2,3);

plot(k/128,real(XEPF));grid; ylabel(‘振幅’);

title(‘其共轭对称序列DFT’);

subplot(2,2,4);

plot(k/128,imag(XOPF));grid; ylabel(‘振幅’);

title(‘其共轭反对称序列DFT’);

四、实验仪器设备

计算机,MATLAB软件

五、实验注意事项

课前预先阅读并理解实验程序;

22

实验五 离散傅立叶变换DFT

六、思考题

1.讨论实验程序1中的.^代表什么含义?

2.讨论实验程序2中圆周移位与线性移位的关系。指令mod的含义是什么? 中极点与系统稳定性的关系?根据程序运行结果判断该系统的稳定性。 3 .计算实验程序3的理论计算结果,对比两种程序运行结果,结果是否一致? 4 .讨论实验程序4说明实序列的共轭对称部分的离散傅立叶变换与原序列的傅立叶变换实部的关系,共轭反对称部分的离散傅立叶变换与原序列的傅立叶变换虚部的关系。

23

实验六 快速傅立叶变换FFT及其应用

实验六 快速傅立叶变换FFT及其应用

一、实验目的

1. 利用MATLAB的快速傅立叶变换来计算信号的离散傅立叶变换。 2. 利用MATLAB程序,理解进一步离散傅立叶变换的物理意义。 3. 利用MATLAB程序,理解快速卷积算法。

二、实验原理

在MATLAB中,使用函数fft可以很容易地计算有限长序列x(n)的离散傅立叶变换X[k]。此函数有两种形式,fft(x)计算序列x(n) 的离散傅立叶变换X(k),这里X(k)的长度与x(n)的长度相等。fft(x,L)计算序列x(n) 的L点离散傅立叶变换,其中L≥N。若L>N,在计算离散傅立叶变换之前,对x(n)尾部的L-N个值进行补零。同样,离散傅立叶变换序列X(k)的离散傅立叶逆变换x(n)用函数ifft计算,它也有两种形式。 (一)、基本序列的离散傅立叶变换计算

N点离散傅立叶变换的一种物理解释就是,X[k]是x(n)以N为周期的周期延拓序列的

~~离散傅立叶级数系数X(k)的主值区间序列,即X(k)?X(k)RN(k)。例如序列

cos(n)RN(n),当N=16时,cos(n)RN(n)正好是cos(n)的一个周期,所以

888cos(????8n)RN(n)的周期延拓序列就是这种单一频率的正弦序列。而当N=8时,

cos(n)RN(n)正好是cos(n)的半个周期,cos(n)RN(n)的周期延拓就不再是单一频888???率的正弦序列,而是含有丰富的谐波成分,其离散傅立叶级数的系数与N=16时的差别很大,因此对信号进行谱分析时,一定要截取整个周期,否则得到错误的频谱。

(二)、验证N点DFT的物理意义

假如x(n)非周期、有限长,则傅立叶变换存在,那么对X(e?k=2πk/N, k=0,1,…,,N-1取样,则可得X(k)。

24

j?)在N个等间隔频率

实验六 快速傅立叶变换FFT及其应用

X(k)?X(?)??2?k/N???x(n)en????j2?kn/N 0?k?N?1

序列x(n)的N点DFT的物理意义是对X(ω)在[0,2π]上进行N点的等间隔采样。

(三)、利用FFT计算序列的线性卷积

直接计算线性卷积计算量大,并且计算机无法判断y(n)的长度,需要计算多少的y(n)值,若输入为无限长,就更无法计算,其运算量随长度成级数增长。由于可以利用FFT对DFT进行有效的计算,我们希望能够利用DFT来计算线性卷积。

设 x(n) 和 h(n) 是长度分别为M和N的有限长序列, 令 L=M+N-1,定义两个长度L的有限长序列: x'(n)??

h'(n)???h(n),?0,0?n?N?1N?n?L?1?x(n),?0,0?n?M?1M?n?L?1 (3.4.8)

(3.4.9)

通过对x(n) 和 h(n)补充零样本值得到上面两个序列。那么:

yl(n)?x(n)?h(n)?yc(n)?x'(n) h'(n) (3.4.10) 上面的过程如下图所示:

计算线性卷积也可以直接调用函数con来计算,因为MATLAB中的计时比较粗糙,所以只有M和N较大的时候,才能比较两种方法的执行时间快慢。

25

实验六 快速傅立叶变换FFT及其应用

三、实验内容与步骤

1. 对复正弦序列x(n)?ej?8nTLAB程序求当N=16和N=8时的离散傅RN(n),利用MA

立叶变换,并显示其图形。 程序:

N=16;N1=8; n=0:N-1;k=0: N1-1; x=exp(j*pi*n/8); X1=fft(x,N); X2=fft(x,N1); subplot(2,1,1); stem(n,abs(X1)); axis([0,20,0,20]); ylabel(‘∣X1(k)∣’); title(‘16点的DFT’);

subplt(2,1, 2);

stem(n, abs(X2)); axis([0, 20, 0, 20]); ylabel(‘∣X2(k)∣’); title(‘8点的DFT’);

2.已知x(n)?R4(n),X(?)?N=16时的DFT。

1?e?j4??j?1?e, 绘制相应的幅频和相频曲线,并计算N=8和

程序:

N1=8;N2=16;

n=0:N1-1;k1=0:N1-1;k2=0:N2-1; w=2*pi*(0:2047)/2048;

26

实验六 快速傅立叶变换FFT及其应用

Xw=(1-exp(-j*4*w))./(1-exp(-j*));%对x(n)的频谱函数采样2048个点可以近似地看

作是连续的频谱

xn=[(n>=0)&(n<4)];%产生x(n) X1k=fft(xn,N1); X2k=fft(xn,N2);

subplot(3,2,1);plot(w/pi,abs(Xw));xlabel(‘w/π’); title(‘x(n)的幅频曲线’);

subplot(3,2,2);plot(w/pi,angle(Xw));axis([0,2,-pi,pi]); line([0,2],[0,0]);xlabel(‘w/π’); title(‘x(n)的相频曲线’);

subplot(3,2,3);stem(k1,abs(X1k),’.’);

xlabel(‘k(w=2πk/N1)’);ylabel(‘∣X1(k)∣’);hold on plot(N1/2*w/pi,abs(Xw));%图形上叠加连续频谱的幅度曲线 subplot(3,2,4);stem(k1,angle(X1k),’.’); axis([0,N1,-pi,pi]);line([0,N1],[0,0]);

xlabel(‘k(w=2πk/N1)’);ylabel(‘Arg[X1(k)]’);hold on

plot(N1/2*w/pi,angle(Xw));%图形上叠加连续频谱的相位曲线 subplot(3,2,5);stem(k2,abs(X2k),’.’);

xlabel(‘k(w=2πk/N2)’);ylabel(‘∣X2(k)∣’);hold on plot(N2/2*w/pi,abs(Xw));%图形上叠加连续频谱的幅度曲线 subplot(3,2,6);stem(k2,angle(X2k),’.’);

27

实验六 快速傅立叶变换FFT及其应用

axis([0,N2,-pi,pi]);line([0,N2],[0,0]);

xlabel(‘k(w=2πk/N2)’);ylabel(‘Ang[X2(k)]’);hold on

plot(N2/2*w/pi,angle(Xw));%图形上叠加连续频谱的相位曲线

3.分别利用快速卷积法以及conv函数计算下面两个序列的线性卷积。

h=[3 2 1 -2 1 0 -4 0 3];

x=[1 -2 3 -4 3 2 1 ] ↑

程序1:快速卷积 clf;

h=[3 2 1 -2 1 0 -4 0 3];%冲激 x=[1 -2 3 -4 3 2 1 ]; %输入序列 L=pow2(nextpow2(length(x)+length(h)-1)); Xk=fft(x,L); Hk=fft(h,L); Yk=Xk.*Hk; y=ifft(Yk,L);

nh=0:8;nx=0:6;ny=0:L-1;

subplot(3,1,1);stem(nx,x);title(‘x(n)’);

subplot(3,1,2);stem(nh,h);title(‘h(n)’);

subplot(3,1,3);stem(nx,x); xlabel(‘时间序号n’);ylabel(‘振幅’);title(‘卷积y(n)’); 程序2: clf;

h=[3 2 1 -2 1 0 -4 0 3];%冲激 x=[1 -2 3 -4 3 2 1 ]; %输入序列 y=conv(h,x); n=0:14; stem(n,y);

28

实验六 快速傅立叶变换FFT及其应用

xlabel(‘时间序号n’);ylabel(‘振幅’); title(‘用卷积函数conv得到的输出’);grid;

四、实验仪器设备

计算机,MATLAB软件

五、实验注意事项

课前预先阅读并理解实验程序;

六、思考题

1.分析实验程序1的图形中,两种N值下DFT是否有差别及产生差别的原因。 2.根据实验程序2的结果图分析由N点DFT的物理意义。

3 .实验程序3两个程序的输出是否一致?计算该卷积的理论结果,与两个程序的输出是否一致?若改变L的长度(改为8时),两个程序的输出是否一致?FFT的变换长度L必须满足什么条件,输出y(n)才等于x(n)和h(n)的线性卷积,这种快速卷积的方法是利用什么运算与线性卷积运算的关系得到的?

29

实验七 基于MATLAB的IIR数字滤波器设计(设计实验)

实验七 基于MATLAB的IIR数字滤波器设计(设计实验)

一、实验目的

1. 进一步熟悉IIR数字滤波器的理论知识。

2. 熟悉与IIR数字滤波器设计有关的MATLAB函数。

3 . 学会通过MATLAB,利用脉冲响应不变法和双线性变换法设计IIR数字滤波器,加深对数字滤波器的常用指标和设计过程的理解。

二、实验原理

(一)、低通滤波器的常用指标:

1??P?G(eG(ej?j?)?1??P,for???P

)??S,for?S????

通带边缘频率:?p,阻带边缘频率:?s , 通带起伏:?p,通带峰值起伏: ?最小阻带衰减:?S??20log(二)、IIR数字滤波器设计

目前,设计IIR数字滤波器的通用方法是先设计相应的低通滤波器,然后再通过双线性变换法和频率变换得到所需要的数字滤波器。模拟滤波器从功能上分有低通、高通、带通及带阻四种,从类型上分有巴特沃斯滤波器、切比雪夫滤波器、椭圆滤波器以及贝塞尔滤波器等。

1、利用模拟滤波器设计IIR数字低通滤波器的步骤。

(1)确定数字低通滤波器的技术指标:通带截止频率ωp、通带衰减αp、阻带截止频率ωs、阻带衰减αs。

(2)将数字低通滤波器的技术指标转换成模拟低通滤波器的技术指标。

脉冲响应不变法:

双线性变换法:

??1tan(?)T2210p??20log10(1??p)[dB],阻带起伏:?s

(?s)[dB]。

???T(3)按照模拟低通滤波器的技术指标设计模拟低通滤波器。

30

实验七 基于MATLAB的IIR数字滤波器设计(设计实验)

(4)将模拟滤波器Ha(s),从s平面转换到z平面,得到数字低通滤波器系统函数H(z)。

2、下面给出与IIR数字滤波器设计有关的MATLAB文件。

(1)buttord.m

用来确定数字低通或模拟低通滤波器的阶次,其调用格式分别是 a. [N,Wn]=buttord(Wp,Ws,Rp,Rs) b. [N,Wn]=buttord(Wp,Ws,Rp,Rs,’s’)

格式a对应数字滤波器,式中Wp,Ws分别是通带和阻带的截止频率,实际上它们是归一化频率,其值在0-1之间,1对应π(即对π的归一化)。Rp,Rs分别是通带和阻带衰减,单位为dB。N是求出的相应低通滤波器的阶次,Wn是求出的3dB频率。

格式b对应模拟滤波器,式中各个变量的含义和格式a相同,但Wp,Ws及Wn是模拟角频率,单位为rad/s。

(2)buttap.m

用来设计模拟低通原型(归一化)滤波器Ha(p),其调用的格式为 [z , p, k]=buttap(N)

N是欲设计的低通原型(归一化)滤波器的阶次,z, p和k分别是设计出Ha(p)的极点、零点及增益。

(3)lp2lp.m

将模拟低通原型(归一化)滤波器Ha(p)转换为实际的低通滤波器Ha(s)。(去归一化),其调用格式为:

[B,A]=lp2lp(b,a,Wn)

b,a分别是模拟低通原型滤波器Ha(p)的分子、分母多项式的系数向量,其中B,A是去归一化后Ha(s) 的分子、分母多项式的系数向量, Wn为截止频率。

(4)bilinear.m

实现双线性变换,即由模拟滤波器Ha(s)得到数字滤波器H(z)。其调用格式是: [Bz,Az]=bilinear(B,A,Fs)

B,A是去归一化后Ha(s) 的分子、分母多项式的系数向量,Bz,Az是H (z) 的分子、分母多项式的系数向量, Fs是抽样频率。

(4)impinvar.m

由脉冲响应不变法将模拟滤波器Ha(s)转换为数字滤波器H(z)。其调用格式

31

实验七 基于MATLAB的IIR数字滤波器设计(设计实验)

是: [Bz,Az]= impinvar(B,A,Fs)

B,A是去归一化后Ha(s) 的分子、分母多项式的系数向量,Bz,Az是H (z) 的分子、分母多项式的系数向量, Fs是抽样频率。

(5) butter.m

用来直接设计巴特沃斯数字滤波器(双线性变换法),实际上它把buttord.m,buttap.m,lp2lp.m及bilinear.m等文件都包含进去,从而使设计过程更简捷,其调用格式为: a. [B,A]=butter(N, Wn)

b. [B,A]=butter(N,Wn,‘s’)

格式a是设计低通数字滤波器,格式b是设计低通模拟滤波器。B,A是H (z) 的分子、分母多项式的系数向量,Wn是截止频率。

三、实验内容与步骤

1. 设计MATLAB程序,采用脉冲响应不变法设计一个巴特沃斯低通数字滤波器,其通带上限临界频率为400Hz,阻带临界频率为600Hz,抽样频率是1000Hz,在通带内的最大衰减为0.3dB, 阻带内的最小衰减为60dB,并绘出幅频特性曲线。 2. 设计MATLAB程序,采用双线性变换法设计一个巴特沃斯低通数字滤波器,要求在通带[0,0.2π]内衰减不大于3dB, 在阻带[0.6π,π]内衰减不小于40dB,并绘出幅频特性曲线。

四、实验仪器设备

计算机,MATLAB软件

五、实验要求

根据要求独立编程设计,并根据程序运行结果写出滤波器的系统函数。

32

实验八 基于MATLAB的FIR数字滤波器设计(设计实验)

实验八 基于MATLAB的FIR数字滤波器设计(设计实验)

一、实验目的

1. 进一步熟悉FIR数字滤波器的理论知识。

2. 熟悉与FIR数字滤波器设计有关的MATLAB函数。 3. 学会通过MATLAB,利用窗函数法设计FIR数字滤波器。

二、实验原理

设计FIR滤波器实际上是要在满足线性相位的条件下,实现幅度响应的逼近。而一个FIR滤波器若是符合线性相位,则必须满足一定的条件,即:

一个FIR滤波器若是线性相位的,则其单位冲激响应必然满足 h(n)??h(N?1?n) n=0,1,…,N-1 h(n)是关于(N-1)/2对称(奇对称或偶对称) 即,

(1) h(n)是偶对称序列

N?1? ? ? ?

?2 ?,0?n?N?1?h?n??h?N?1?n? (1) h(n) 是奇对称(反对称)序列

N?1??? ? 2 ???h?n???h?N?1?n???

设滤波器要求的理想频率响应为Hd(ejw) , 那么FIR滤波器的设计问题在于

N?1——寻找一系统函数H(z)??h(n)zn?0?n ,使其频率响应H(ejw)?H(z)|z?ejw逼近

Hd(ejw)。若要求FIR滤波器具有线性相位特性,则h(n)必须满足上节所述的对称条件。逼近的方法有三种:

窗口设计法(时域逼近);频率采样法Frequency-sampling(频域逼近);

33

最优化设计Optimum Equiripple(等波纹逼近)。

实验八 基于MATLAB的FIR数字滤波器设计(设计实验)

窗函数法又称傅立叶级数法,是设计FIR数字滤波器的最简单的方法。FIR数字滤波器的设计问题就是要使所设计的FIR数字滤波器的频率响应H(w)去逼近所要求的理想滤波器的响应Hd(w)。从单位采样响应序列来,就是使所设计滤波器的h(n)逼近理想单位采样响应序列hd(n),这可以用hd(n)和一个窗函数w(n)的乘积来得到。 (一)、设计原理。

窗函数设计FIR数字滤波器的步骤如下: (1)给定要求的频率响应函数Hd(w); (2)计算hd(n);

(3)根据过渡带宽及阻带最小衰减的要求,选定窗的性状以及窗的大小N; (4)根据所选择的合适的窗函数w(n)来修正hd(n),得到所设计的FIR数字滤波器的单位采样响应序列h (n)= hd(n) w(n),n=0,1,…,N-1 。 (二)、函数的应用

MATLAB中用fir1函数来设计具有标准频率响应的FIR滤波器。其调用方式: b=fir1(n,wn)——设计n阶低通FIR滤波器,返回的向量b为滤波器的系数(即h(n)的值),它的阶数为n+1;截止频率为wn(对π归一化后的值)。 b=fir1(n,wn,’hign’)——设计n阶高通FIR滤波器 b=fir1(n,wn,’low’)——设计n阶低通FIR滤波器 b=fir1(n,wn,’bandpass’)——设计n阶带通FIR滤波器 b=fir1(n,wn,’stop’)——设计n阶带阻FIR滤波器

b=fir1(n,wn,win)——输入参数win用来指定使用的窗函数的类型,其长度为n+1,缺省情况下,默认为汉明窗。

三、实验内容与步骤

1. 用矩形窗、三角窗、汉宁窗、汉明窗分别设计低通数字滤波器。信号采样频率为1000Hz,数字滤波器的截止频率为 100Hz,滤波器的阶数为80。

程序:

passrad=2*pi*100/1000;

w1=boxcar(81); %矩形窗 w2=triang(81); %三角窗 w3=hanning(81); %汉宁窗

34

实验八 基于MATLAB的FIR数字滤波器设计(设计实验)

w4=hamming(81); %矩形窗 n=1:1:81;

hd=sin(passrad*(n-41))./(pi*(n-41)); hd(41)= passrad/pi; h1=hd. *rot90(w1); h2=hd. *rot90(w2); h3=hd. *rot90(w3); h4=hd. *rot90(w4); [MAG1,RAD]=freqz(h1); [MAG2,RAD]=freqz(h2); [MAG3,RAD]=freqz(h3); [MAG4,RAD]=freqz(h4); subplot(2,2,1);

plot(RAD,20*log10(abs(MAG1))); grid on; subplot(2,2,2);

plot(RAD,20*log10(abs(MAG2))); grid on; subplot(2,2,3);

plot(RAD,20*log10(abs(MAG3))); grid on; subplot(2,2,4);

plot(RAD,20*log10(abs(MAG4))); grid on;

运行上面程序,并显示其图形。

2.编写MATLAB程序,利用窗函数法设计线性相位FIR低通数字滤波器,实现对模拟信号采样后进行数字低通滤波,对模拟信号的滤波要求如下:

通带截止频率:fp=2kHz;阻带截止频率:fs=3kHz; 阻带最小衰减:αs=40dB;

采样频率:Fs=10kHz

选择合适的窗函数及其长度,求出h(n),并画出幅频特性衰减曲线。

35

实验八 基于MATLAB的FIR数字滤波器设计(设计实验)

3.编写MATLAB程序,利用窗函数法设计一线性相位FIR数字低通滤波器,截止频率wc=0.2π,过渡带宽 △w<0.4π, 阻带衰减αs>40dB。

四、实验仪器设备

计算机,MATLAB软件

五、实验注意事项

根据要求独立编程设计。

六、思考题

1. 根据实验程序1的幅频特性曲线,分析不同窗函数下的幅频特性曲线有什么差别。

2. 根据实验内容2和3所编写程序的运行结果,讨论它们的幅频曲线是否符合设计要求?它们的h(n)有什么特点?

36

本文来源:https://www.bwwdw.com/article/atah.html

Top