《数字信号处理》实验指导.docx
《《数字信号处理》实验指导.docx》由会员分享，可在线阅读，更多相关《《数字信号处理》实验指导.docx（20页珍藏版）》请在冰豆网上搜索。
《数字信号处理》实验指导
《数字信号处理》实验指导
康文静
一、Matlab与信号处理有关的基本函数
1.单位采样序列——zeros函数
单位采样序列，或称为单位脉冲序列（编程序常用Delta表示）：
如果
在时间轴上延迟了k个单位，得到
，即：
Matlab函数zeros（1,N）可产生一个包含N个零的行向量。
下图为产生delta函数的几种编程方法。
2.单位阶跃序列（Unitstepsequence）——函数ones
单位阶跃序列：
在MATLAB中可以利用ones（）函数实现，
函数ones（1,N）可产生一个包含N个1的行向量。
下图为产生单位阶跃序列的几种编程方法。
3.矩形序列——rectanglesequence
矩形序列：
编程：
>>rect=[zeros（1,3）,ones（1,4）,zeros（1,3）]
rect=
0001111000
4.实指数序列——stem函数，幂的表示
实指数序列：
用matlab实现
，并画出相应图形。
>>n=[0:
8];
>>x=（0.6）.^n;
>>stem（n,x）;
5.正弦序列——函数plot，sin，cos；常量pi
用matlab实现如下序列，并画出相应图形。
>>n=[0:
0.1:
10];
>>x=2*sin（0.6*pi*n）+cos（0.1*pi*n+pi/4）;
>>plot（n,x）;%对应第一幅图
>>stem（n,x）;%对应第二幅图
6.复指数序列——exp函数
复指数序列：
在MATLAB中实现：
用matlab实现
。
>>n=1:
6;
>>x=exp（j*0.5*n）
x=0.8776+0.4794i0.5403+0.8415i0.0707+0.9975i-0.4161+0.9093i-0.8011+0.5985i-0.9900+0.1411i
7.序列的翻转——fliplr函数
，用matlab实现
，并画出相应图形。
>>n=[0:
0.1:
10];
>>x=2*sin（0.6*pi*n）+cos（0.1*pi*n+pi/4）;
>>y=fliplr（x）;
>>n=-fliplr（n）;
>>plot（n,y）;
>>stem（n,y）;
8.计算序列的线性卷积——函数conv
已知，
程序如下：
>>x=[1,1,2,3];
>>h=[2,1,8,6,3]
>>y=conv（x,h）
y=2313222839249
9.序列的傅里叶变换（DTFT）-freqz函数
freqz函数计算序列的离散时间傅里叶变换在给定的离散频率点上的抽样值。
假设：
freqz的调用方式如下：
（1）[H,w]=freqz（b,a,N）,其中b和a分别是分子和分母多项式的系数向量。
此函数在单位圆上半部[0,pi]上等间隔地计算N个频率响应。
返回该系统的N点频率响应矢量w和N点复数频率响应矢量H，如果N没有说明，缺省值为512。
（2）H=freqz（b,a,w）,返回矢量w指定的那些频率点上的频率响应，频率范围为[0,pi]。
（3）[H,w]=freqz（b,a,N，‘whole’）,在整个单位圆上等间隔地计算N个频率响应，即频率范围为[0,2pi]。
编程题目：
因果系统
，利用freqz函数试画出
的幅频相应
和相频相应
。
利用函数abs/angle/real/image等计算DTFT的幅度、相位以及实部和虚部。
程序如下：
>>b=[0.5];a=[1,-0.85];
>>[H,w]=freqz（b,a,200,'whole'）;
>>magH=abs（H（1:
101））;%计算幅度
>>phaH=angle（H（1:
101））;%计算相位
>>w=w（1:
101）;
>>%画图
>>subplot（2,1,1）;
>>plot（w/pi,magH）;
>>grid;%加网格
>>xlabel（'frequencyUnit:
pi'）;
>>ylabel（'Magnitude'）;
>>title（'MagnitudeResponse'）;
>>subplot（2,1,2）;
>>plot（w/pi,phaH/pi）;
>>grid;
>>xlabel（'frequencyUnit:
pi'）;
>>ylabel（'PhaseUnit:
pi'）;
>>title（'PhaseResponse'）;
10.Z变换——函数tf2zp/zp2tf/zplane
假设：
利用函数[z,p,k]=tf2zp（b,a）可确定分子分母多项式按z的降幂排列的有理z变换式的零点z/极点p/增益常数k。
函数[b,a]=zp2tf（z,p,k）用来实现相反的过程。
函数zplane可以用来画出z变换的零极点图，有如下两种调用方式：
zplane（zeros，poles）
zplane（b，a）
编程题目：
已知离散系统差分方程为：
，
求其系统函数，画出零/极点示意图，并判断系统的稳定性。
由差分方程可得
程序如下：
>>b=[1,-0.5,0];
>>a=[1,0.75,0.125];
>>[z,p,k]=tf2zp（b,a）;
>>zplane（b,a）
>>title（'Pole-ZeroPlot'）
由上图可知，系统的极点全部在单位圆内，因此系统是稳定的。
11.求解差分方程——fliter函数
filter函数求解给定输入x（n）的差分方程的解，该函数调用形式为：
y=filter（b,a,x）
编程题目：
已知离散系统差分方程为：
计算并画出冲击响应h（n）（n=-10,-9,…,0,1,…,50）;由h（n）确定系统是否稳定。
程序如下：
>>b=[1];a=[1,-1,0.5];
>>%求单位脉冲响应
>>n=[-10:
50];x=[（n）==0];
>>h=filter（b,a,x）
>>stem（n,h）
>>xlabel（'n'）;
>>ylabel（'h（n）'）;
>>title（'ImpulseResponse'）;
>>axis（[-10,50,-1,1.5]）
>>sum（abs（h））%
12.计算信号的DFT和IDFT——函数fft和ifft
调用形式如下：
●fft（x）函数用来计算M点DFT，其中M是序列x的长度；
●fft（x,N）计算N点DFT，N是用户指定的长度。
若序列x的长度M>N，则将序列截短为N点序列，再作N点DFT；
若序列x的长度M如果N为2的正整数幂，则可以得到基2fft算法，如果N既不是2的正整数幂，也不是质数，则函数将N分解为质数，得到较慢的组合数FFT算法（混合基）；如果N为质数，则fft函数采用的是原来的DFT算法。
●ifft（x）用来计算M点IDFT，其中M是序列x的长度；
●ifft（x,N）计算N点IDFT，N是用户指定的长度。
若序列x的长度M>N，则将序列截短为N点序列，再作N点IDFT；
若序列x的长度M编程题目：
已知一长度为16的有限长序列
，计算其16点和256点DFT。
程序如下：
>>n=0:
15;
>>x=sin（0.25*pi*n）;%序列x（n）
>>dft_16=fft（x）;%16点DFT
>>dft_256=fft（x,256）;%256点DFT
>>L=0:
255;
>>plot（L/256,abs（dft_256））;
>>xlabel（'normalizedfrequency'）;
>>ylabel（'Magnitude'）;
>>title（'256点DFT'）;
>>holdon;
>>plot（n/16,abs（dft_16）,'*'）;
>>title（'256点DFT和16点DFT（*）'）;
编程题目：
已知
，为16点的有限长序列，计算其16点和256点DFT。
程序如下：
>>n=0:
15;
>>x=sin（0.25*pi*n）;%序列x（n）
>>dft_16=fft（x）;%16点DFT
>>dft_256=fft（x,256）;%256点DFT
>>L=0:
255;
>>subplot（2,1,1）;stem（L/256,abs（dft_256））;
>>xlabel（'normalizedfrequency'）;
>>ylabel（'Magnitude'）;
>>title（'256点DFT'）;
>>subplot（2,1,2）;stem（n/16,abs（dft_16）,'*'）;
>>xlabel（'normalizedfrequency'）;
>>ylabel（'Magnitude'）;
>>title（'16点DFT（*）'）;
二、基于Matlab的数字滤波器设计
1.IIR滤波器设计函数
（1）巴特沃斯模拟滤波器
[N，Wn]=buttord（Wp,Ws,Rp,Rs,‘s’）
此函数可获得巴特沃斯滤波器的参数N和Wn，分别是在给定通带边界频率Wp、阻带边界频率Ws、通带最大衰减Rp和阻带最小衰减Rs的条件下，所需要的Butterworth滤波器的最小阶数和对应的3dB带宽，这里Wp和Ws的单位是rad/s。
[z,p,k]=buttap（N）
求出N阶巴特沃斯低通滤波器的零点、极点和增益。
其中p是长度为N的列向量，分布在单位圆的左半平面，而零点z是空矩阵。
[B,A]=butter（N,Wn,‘s’）
获得巴特沃斯低通滤波器系统函数的分子多项式B和分母多项式A的系数。
例：
设计一个模拟巴特沃斯低通滤波器满足以下指标：
通带截止频率为fp=1Hz，通带最大衰减1dB（αp=1）；阻带截止频率fs=4Hz，阻带最小衰减30dB（αs=30）。
>>%Butterworth低通滤波器设计
>>Wp=2*pi*1;Ws=2*pi*4;
>>Rp=3;Rs=30;
>>%设计滤波器
>>[N,Wn]=buttord（Wp,Ws,Rp,Rs,'s'）
N=3
Wn=7.9490
>>[z,p,k]=buttap（N）
z=[]
p=
-0.5000+0.8660i
-0.5000-0.8660i
-1.0000
k=1.0000
>>[B,A]=butter（N,Wn,'s'）
B=000502.2695
A=1.000015.8980126.3731502.2695
因此要满足所要求的指标需要3阶巴特沃斯滤波器，其系统函数为：
>>画图Butterworth低通滤波器
>>f1=linspace（0,Wp,5）
>>f2=linspace（Wp,Ws,15）
>>f3=linspace（Ws,2*pi*10,30）
>>h1=20*log10（abs（freqs（B,A,f1）））
>>h2=20*log10（abs（freqs（B,A,f2）））
>>h3=20*log10（abs（freqs（B,A,f3）））
>>plot（[f1f2f3]/（2*pi）,[h1,h2,h3]）
>>grid
>>xlabel（'FrequencyinHz'）
>>ylabel（'GainindB'）
此外，chebwin,cheb1ord,cheb2ord,cheby1,cheby2,ellip,ellipord等函数可分别设计切比雪夫和椭圆滤波器。
（2）脉冲响应不变法的Matlab实现
[bz,az]=impinvar（b,a,Fs）
此函数实现用脉冲响应不变法将模拟滤波器转换为数字滤波器。
其中b和a分别是模拟滤波器的系统函数H（s）的分子多项式和分母多项式的系数，Fs是脉冲响应不变法的采样频率，单位为Hz。
如果Fs没有说明，其缺省值为1Hz。
运算结果bz、az分别表示数字滤波器的系统函数H（z）的分子多项式和分母多项式的系数。
例：
利用impinvar函数将模拟系统函数Ha（s）转变为数字系统函数H（z），采样周期为Ts=0.5。
其中
。
>>b=[0,3];a=[1,4,3];T=0.5;Fs=1/T;
>>[bz,az]=impinvar（b,a,Fs）
bz=
00.2876
az=
1.0000-0.82970.1353
因此，数字滤波器为
（3）双线性变换法的Matlab实现
[zd,pd,kd]=bilinear（z,p,k,Fs）
其中z，p，k分别表示模拟滤波器系统函数的零点、极点和增益，Fs是采样频率，zd,pd,kd是对应的数字滤波器系统函数的零点、极点和增益。
[numd，dend]=bilinear（num，den，Fs）
其中num，den和Fs分别表示模拟滤波器系统函数分子、分母多项式系数，Fs是采样频率；numd，dend是对应的数字滤波器系统函数的分子、分母多项式系数。
例：
利用双线性变换法设计低通巴特沃斯数字滤波器，其技术指标为：
Wp=0.2*pi,Rp=1dB;Ws=0.4*pi,Rs=15dB,T=1s.
>>Wp=0.2*pi;Rp=1;Ws=0.4*pi;Rs=15;T=1,
>>Fs=1/T
>>omegap=（2/T）*tan（Wp/2）;%通带截止频率预畸
>>omegas=（2/T）*tan（Ws/2）;%阻带截止频率预畸
>>[N,Wn]=buttord（omegap,omegas,Rp,Rs,'s'）
>>[B,A]=butter（N,Wn,'s'）
>>[b,a]=bilinear（B,A,Fs）
>>%画图
>>[h,w]=freqz（b,a,256）;
>>h1=20*log10（abs（h））;%增益（dB）
>>plot（w/pi,h1）;
>>grid
>>xlabel（'DigitalFrequencyinpiunits'）
>>ylabel（'GainindB'）
>>axis（[01-8010]）
2.FIR滤波器设计函数（窗函数法）
（1）常用的窗函数
（a）矩形窗
调用格式：
w=boxcar（n），根据长度n产生一个矩形窗w。
（b）三角窗
调用格式：
w=triang（n），根据长度n产生一个三角窗w。
（c）汉宁窗
调用格式：
w=hanning（n），根据长度n产生一个汉宁窗w。
（d）哈明窗
调用格式：
w=hamming（n），根据长度n产生一个哈明窗w。
（e）布莱克曼窗
调用格式：
w=blackman（n），根据长度n产生一个布莱克曼窗w。
（f）凯塞窗
调用格式：
w=kaiser（n，beta），根据长度n产生一个布莱克曼窗w，beta是控制窗形状的参数。
上述窗函数返回的变量w是一个长度为n的列向量，表示窗函数在这n点的取值。
（2）窗函数法设计FIR滤波器
（a）函数fir1
调用格式：
b=fir1（n,wn,’ftype’,Window）
其中，n为阶数，因此h（n）长度为n+1。
wn是截止频率，其取值在0~1之间，这是以采样频率为基准频率的标称值，故1对应采样频率；若Wn是标量，则用来设计低通滤波器；如果要设计带通滤波器，输入wn=[w1,w2]为一矢量。
ftype是滤波器的类型，低通时省略该参数，高通为high，带阻为stop。
Window表示选用的窗函数类型，以列向量形式表示，向量Window的长度必须为n+1，若缺省，则fir1默认使用哈明窗。
B对应设计好的滤波器的单位取样响应h（n）。
（a）函数fir2
调用格式：
b=fir2（n,f,m,Window）
该函数采用窗函数法设计具有任意频率响应的FIR数字滤波器。
其中f是频率向量，其值在0~1（标称值）之间，1对应采样频率，其中第一个点必须是0，最后一个点必须是1，而且频率点必须是递增的。
m对应于频率点f处期望的幅频响应，f和m的长度必须相等。
例：
设计长度为8的线性相位FIR滤波器。
其理想幅频特性满足
>>Window=boxcar（8）;%长度为8的矩形窗Window
>>b=fir1（7,0.4,Window）;%用矩形窗作为脉冲响应的窗函数
>>freqz（b,1）；%矩形窗对应的频率响应
freqz（B,A,...）withnooutputargumentsplotsthemagnitudeandunwrappedphaseofthefilterinthecurrentfigurewindow.
>>window=blackman（8）;%长度为8的布拉克曼窗Window
>>b=fir1（7,0.4,Window）;
>>freqz（b,1）
例：
设计线性相位带通滤波器，阶数N=15，上下边带截至频率分别为0.3*pi，0.5*pi。
>>window=blackman（16）;
>>b=fir1（15,[0.30.5],window）;
>>freqz（b,1）
例：
用矩形窗和哈明窗设计线性相位FIR低通滤波器，wc=0.25*pi，N=10。
>>N=10;
>>M=128;
>>b1=fir1（N,0.25,boxcar（N+1））;
>>b2=fir1（N,0.25,hamming（N+1））;
>>h1=freqz（b1,1,M）;
>>h2=freqz（b2,1,M）;
>>f=0:
0.5/M:
0.5-0.5/M;
>>plot（f,abs（h1）,'-.',f,abs（h2））;%画出幅频响应
>>legend（'矩形窗','哈明窗'）
>>grid
二、建立一个M-file的方法
（1）建立一个M-file
（2）编写M-file，实现阶跃序列程序。
（3）保存
升级会员 