数字信号处理-MATLAB参考资料.doc

数字信号处理 MATLAB参考资料

———————————————————————————————— 作者：
———————————————————————————————— 日期：

23

个人收集整理 勿做商业用途
MATLAB相关内容一 信号与系统
一、 目的
1、 熟悉连续信号理想采样前后的频谱变化，加深对时域采样定理的理解.
2、 熟悉时域
二、 内容
1、理想采样信号序列
（1）首先产生信号x(n），0<=n〈=50
n=0:50; % 定义序列的长度是50
A=444。128; ％ 设置信号有关的参数
a=50*sqrt （2.0）＊pi；
T=0.001; % T为采样率
w0=50*sqrt(2.0)*pi； ％ ω符号在MatLab 中不能输入,用w 代替
x=A*exp（—a*n*T)。*sin(w0*n＊T）; %pi 是MATLAB 定义的π，信号乘可采用“。* "
close all ％ 清除已经绘制的x(n） 图形
subplot（3，1，1）;stem(x)； % 绘制x(n) 的图形
title（'理想采样信号序列')；
（2）绘制信号x(n）的幅度谱和相位谱
k=-25：25;
X=x*（exp(-j*pi/12.5)).^（n’＊k);
magX=abs（X)； % 绘制x(n） 的幅度谱
subplot(3，1,2)；stem（magX）;title（’理想采样信号序列的幅度谱’）;
angX=angle(X）； % 绘制x（n） 的相位谱
subplot(3，1,3）;stem（angX） ； title ('理想采样信号序列的相位谱’) ；

(3）改变参数为： A=1 , a=0。4, w0=2.0734, T=1
A=1; a=0.4; w0=2.0734； T=1； ％ 设置信号有关的参数和采样率T
x=A＊exp(-a*n*T).＊sin(w0＊n＊T)； ％pi 是MATLAB 定义的π，信号乘可采用“。* "
close all % 清除已经绘制的x（n) 图形
subplot(3，1，1);stem（x); % 绘制x（n) 的图形
title('理想采样信号序列’）；
k=-25:25;
X=x＊(exp（—j*pi/12。5)）。^(n’＊k）;
magX=abs(X)； % 绘制x(n) 的幅度谱
subplot(3,1,2)；stem(magX);title（’理想采样信号序列的幅度谱’)；
angX=angle（X）; % 绘制x(n) 的相位谱
subplot（3,1,3)；stem（angX) ; title ('理想采样信号序列的相位谱’） ;

2、单位脉冲序列
在MatLab 中，这一函数可以用zeros 函数实现：
n=1:50; ％ 定义序列的长度是50
x=zeros（1，50); ％ 注意:MATLAB 中数组下标从1 开始
x（1）=1；
close all；
subplot（3，1,1);stem(x);title(’单位冲击信号序列')；
k=—25：25；
X=x＊（exp(-j＊pi/12.5)).^(n'*k）；
magX=abs（X）; % 绘制x(n) 的幅度谱
subplot(3,1，2）;stem(magX）；title(’单位冲击信号的幅度谱'）;
angX=angle（X); ％ 绘制x(n） 的相位谱
subplot(3,1，3）;stem(angX) ； title （'单位冲击信号的相位谱’) ;

3、矩形序列
n=1:50；
x=sign（sign（10—n）+1);
close all;
subplot（3,1，1)； stem(x）;title(’矩形信号序列’）；
k=—25:25;
X=x*(exp（—j*pi/25））。^（n’＊k)；
magX=abs（X); % 绘制x(n) 的幅度谱
subplot(3,1，2）;stem（magX)；title(’矩形序列的幅度谱'）；
angX=angle(X）； % 绘制x（n） 的相位谱
subplot（3,1，3）；stem(angX） ; title （’矩形序列相位谱') ;

4、特定冲击串
n=1：50； % 定义序列的长度是50
x=zeros（1,50); % 注意：MATLAB 中数组下标从1 开始
x(1）=1； x（2)=2.5； x（3）=2.5; x(4)=1；
close all；
subplot(3，1，1);stem（x)；title(’特定冲击串序列'）；
k=—25:25；
X=x＊(exp（—j＊pi/12.5)）.^(n'*k）；
magX=abs（X)； ％ 绘制x（n) 的幅度谱
subplot(3,1,2）;stem（magX)；title('特定冲击串序列的幅度谱'）;
angX=angle（X）; ％ 绘制x（n) 的相位谱
subplot（3，1，3);stem（angX） ； title ('特定冲击串序列相位谱’） ；

MATLAB相关内容二 卷积计算及定理
一、 目的
利用卷积方法观察分析信号、系统的频谱特性
二、内容
1、卷积计算
在MATLAB 中,提供了卷积函数conv，即y=conv(x,h），调用十分方便。
n=1:50； % 定义序列的长度是50
hb=zeros(1，50)； ％ 注意：MATLAB 中数组下标从1 开始
hb(1)=1; hb（2)=2。5； hb(3）=2。5； hb(4）=1；
close all;
subplot（3，1,1);stem(hb);title('系统hb［n］');
m=1：50; ％ 定义序列的长度
T=0.001; ％ 定义序列的采样率
A=444。128; %设置信号有关的参数
a=50*sqrt(2.0）＊pi；
w0=50＊sqrt(2。0)*pi；
x=A＊exp（—a＊m*T)。*sin(w0*m＊T)； ％pi 是MATLAB 定义的π，信号乘可采用“.＊ ”
subplot(3,1，2);stem（x)；title（’输入信号x[n］’）；
y=conv（x，hb)；
subplot（3，1，3)；stem(y）；title（’输出信号y［n]');

2、卷积定律验证
(1）
n=1:50; % 定义序列的长度是50
hb=zeros(1，50)； % 注意:MATLAB 中数组下标从1 开始
hb(1)=1； hb（2)=2.5； hb(3)=2。5; hb(4）=1；
m=1:50； ％ 定义序列的长度
T=0。001; % 定义序列的采样率
A=444。128； ％设置信号有关的参数
a=50*sqrt（2。0）*pi;
w0=50＊sqrt（2.0）*pi；
x=A＊exp(-a＊m*T).＊sin(w0*m*T); ％pi 是MATLAB 定义的π，信号乘可采用“.＊ "
y=conv（x,hb)；
k=—25：25;
X=x*（exp（—j*pi/12。5)）。^（n'*k);
magX=abs（X)； ％ 绘制x(n） 的幅度谱
subplot(3，2，1);stem(magX）；title（'入信号的幅度谱'）;
angX=angle（X); ％ 绘制x（n) 的相位谱
subplot（3，2，2)；stem（angX） ； title （’输入信号的相位谱'）
Hb=hb＊（exp（—j＊pi/12.5)）。^（n’＊k）;
magHb=abs（Hb); % 绘制hb（n） 的幅度谱
subplot（3,2,3)；stem（magHb）;title（'系统响应的幅度谱'）;
angHb=angle（Hb）; ％ 绘制hb(n） 的相位谱
subplot（3,2,4）;stem（angHb) ； title （'系统响应的相位谱’） ;
n=1：99； k=1:99；
Y=y＊（exp(-j*pi/12.5))。^（n'＊k);
magY=abs（Y）； ％ 绘制y（n) 的幅度谱
subplot(3,2，5)；stem(magY）;title('输出信号的幅度谱'）;
angY=angle（Y）; % 绘制y（n) 的相位谱
subplot(3,2,6);stem（angY） ; title （'输出信号的相位谱’)

(2)
n=1:50; ％ 定义序列的长度是50
hb=zeros（1，50)； ％ 注意:MATLAB 中数组下标从1 开始
hb（1）=1; hb（2）=2。5； hb(3）=2。5； hb（4）=1；
m=1:50; ％ 定义序列的长度
T=0.001; ％ 定义序列的采样率
A=444。128; ％设置信号有关的参数
a=50＊sqrt（2.0）*pi；
w0=50*sqrt（2.0）*pi；
x=A＊exp(—a＊m＊T）.＊sin（w0＊m＊T)； %pi 是MATLAB 定义的π，信号乘可采用“。* ”
y=conv(x，hb）；
k=—25:25；
X=x＊(exp（—j*pi/12.5)).^（n’＊k);
Hb=hb＊（exp(—j*pi/12.5））。^（n'＊k);
n=1：99； k=1:99;
Y=y＊(exp(—j＊pi/12.5)）。^(n’＊k）；
XHb=X。*Hb;
Subplot（2，1，1）；stem(abs（XHb））;title（’x(n)的幅度谱与hb(n)幅度谱相乘’）;
Subplot(2,1，2）；stem(abs(Y))；title(’y（n）的幅度谱')；

MATLAB相关内容三 用FFT作频谱分析
一、 目的
掌握用FFT方法对连续信号和离散信号进行频谱分析
二、 内容
1、 高斯序列 p=8、q=2时高斯序列的频谱特性
n=0:15;
p=8；
q=2；
x =exp（-1＊（n—p)。^2/q)；
close all；
subplot(3，1,1);
stem（fft（x)) ； ％ 利用fft 函数实现命傅氏变换
subplot（3,1，2);
stem（abs（fft（x)));
subplot(3，1,3）；
stem(angle（fft(x))）;

2、衰减正弦序列
n=0：15； ％定义序列的长度是15
a=0.1;
f=0.0625；
x=exp（-a＊n)。*sin(2*pi＊f*n)；
close all;
subplot（2，1,1）；
stem(x）;
title（’衰减正弦序列'）；
subplot(2，1,2);
stem(abs(fft(x）））；
title('x信号的频谱')；

3、三角序列
for i=1：4
x(i）=i;
end
for i=5:8
x（i）=9—i；
end
close all
subplot(2，1,1);stem（x）；
subplot(2，1，2);stem（abs(fft（x))）;

4、反三角
for i=1：4
x（i）=5—i；
end
for i=5：8
x（i)=i-4;
end
close all
subplot（2，1，1）；stem（x)；
subplot（2,1,2);stem（abs（fft(x）））;

MATLAB相关内容四 IIR滤波器的实现
（IIR滤波器是无限脉冲响应数字滤波器的简称)
一、 目的
熟悉用双线性变换法设计IIR数字滤波器和方法。

二、内容
（一）基本概念
数字滤波器是指输入、输出均为数字信号，通过一定运算关系，改变输入信号所含频率成分的相对比例或者滤除某些频率成分的器件.数字滤波器与模拟滤波器概念相同，只是信号的形式和实现滤波的方法不同。一般数字滤波器从现实的网络结构或者从单位脉冲响应分类，可以分成无限脉冲脉冲响应(IIR）滤波器和有限脉冲响应（FIR）滤波器。
IIR滤波器设计的主要方法是先设计低通模拟滤波器，进行频率变换，将其转换为相应的高通、带通模拟滤波器，再将模拟滤波器转换为相应的数字滤波器。
N阶模拟滤波器系统函数的一般形式为

N阶数字模拟滤波器系统函数的一般形式为

(二）求滤波器频率响应的函数
1、freqs函数:求模拟滤波器的频率响应
freqs(b,a,w）计算由向量w(rad/s)指定的频率点上（或W直接为采样点数）频率响应，其中b和a分别为模拟滤波器系统函数H（s）的分子与分母。freqs函数将自动绘出幅频和相频曲线.
实验4—1：系统传输函数为的模拟滤波器，绘出其幅频和相频谱.
a=[1，0.4，1］;
b=［0。2，0.3，1];
w=logspace（—1,1); ％产生从10—1到101之间地50个等间距点,即50个频率点
freqs（b，a,w)；
运行得幅频特性和相频特性如下：

2、freqz函数:求数字滤波器的频率响应
freqz（b,a，w)计算由向量w指定的频率点上模拟滤波器H(z)的频率响应，其中b和a分别为数字滤波器系统函数H(z）的分子与分母。freqz函数自动将频率点均匀设置在频率范围[0,Π］上，并将自动绘出幅频和相频曲线.
实验4—2：系统传输函数为的数字滤波器,绘出其幅频和相频谱。

a=[1,0.4,1］；
b=[0.2,0.3,1］;
freqz（b，a,128）；
运行得幅频特性和相频特性如下：

三、ButterWorth（巴特沃斯）模拟和数字滤波器的设计
1、直接设计
（1）buttord函数：ButterWorth滤波器阶数的选择(即ButterWorth滤波器的最小阶数的计算）
ButterWorth模拟滤波器： 调用格式：［n,Wn]=buttord(Wp，Ws,Rp，Rs)
ButterWorth数字滤波器： 调用格式:［n，Wn］=buttord(Wp,Ws,Rp，Rs,’s’),其中’s’为变元
在给定滤波器性能的情况下（通带临界频率Wp、阻带临界频率Ws、通带由最大衰减Rp
和阻带内最小衰减Rs）,计算ButterWorth滤波器的阶数n和截止频率Wn.
（2）butter函数：ButterWorth滤波器设计
ButterWorth模拟滤波器： 调用格式:［b，a]=butter（n,Wn)
ButterWorth数字滤波器： 调用格式：［b，a]=butter（n,Wn，’s’），其中’s’为变元
根据滤波器的阶数n和截止频率Wn计算ButterWorth滤波器分子和分母系数（b为分子系数的矢量形式，a为分母系数的矢量形式）。
实验4—3： 采用直接设计方法设计一个ButterWorth数字滤波器满足以下参数:采样频率为1Hz,通带临界频率Wp=0.2Hz、通带内衰减小于1dB（Rp=1）;阻带临界频率Ws=0。3Hz、阻带内衰减大于25dB（Rs=25).

［n,Wn］=buttord（0.2，0。3,1，25）；
［b,a]=butter(n，Wn);
freqz（b，a，128);
运行得幅频特性和相频特性如下:

2、脉冲响应不变法设计ButterWorth(巴特沃斯）数字滤波器
（设计思路:先设计模拟滤波器，再模拟滤波器转换成数字滤波器）
调用格式:［bz,az]=impinvar（b，a,Fs),在给定模拟滤波器参数b、a和采样频率Fs的前提下,把模拟滤波器的参数变为近似等价的数字滤波器参数，从而设计数字滤波器.
实验4-4:采用脉冲响应不变法设计上述ButterWorth(巴特沃斯)数字滤波器
［n,Wn］=buttord（0。2，0。3,1,25,'s’)；
［b,a］=butter(n,Wn，'s'）；
[bz，az]=impinvar（b,a，1）;
freqz(bz,az，128）；
运行得幅频特性和相频特性如下:

3、双线性不变法设计ButterWorth（巴特沃斯)数字滤波器
（设计思路：先设计模拟滤波器，再模拟滤波器转换成数字滤波器）
调用格式:[bz，az］=bilinear（b,a,Fs），在给定模拟滤波器参数b,a和采样频率Fs的前提下,把模拟滤波器的参数变为近似等价的数字滤波器参数,从而设计数字滤波器.
实验4-5:采用双线性不变法设计上述ButterWorth（巴特沃斯)数字滤波器

[n,Wn］=buttord(0。2，0。3，1,25，’s’);
[b，a]=butter（n，Wn，’s’);
[bz，az］=bilinear(b，a,1)；
freqz（bz,az，128）;
运行得幅频特性和相频特性如下：

四、Chebyshev(切比雪夫)Ⅰ型模拟和数字滤波器的设计
(1）cheb1ord函数：Chebyshev滤波器阶数的选择（即Chebyshev滤波器最小阶数的计算）
Chebyshev数字滤波器: 调用格式：[n，Wn］= cheb1ord （Wp，Ws,Rp，Rs）， 在给定滤波器性能的情况下（通带临界频率Wp、阻带临界频率Ws、通带内波纹Rp和阻带内衰减Rs），计算Chebyshev数字滤波器的阶数n和截止频率Wn.
（2）cheby1函数：Chebyshev滤波器设计
Chebyshev数字滤波器: 调用格式:［b，a]=cheby1(n， Rp ，Wn）， 根据滤波器的阶数n、通带内波纹Rp和截止频率Wn计算Chebyshev数字滤波器分子和分母系数(b为分子系数的矢量形式，a为分母系数的矢量形式）.

实验4-6： 设计一个Chebyshev数字滤波器满足以下参数:采样频率为1Hz，通带临界频率Wp=0。2Hz、通带内衰减小于1dB（Rp=1)；阻带临界频率Ws=0.3Hz、阻带内衰减大于25dB(Rs=25).(参数同前）
[n，Wn]= cheb1ord （0。2,0。3，1，25）；
［b，a]= cheby1(n，1, Wn)；
freqz(b，a,128)；
运行得幅频特性和相频特性如下:

MATLAB相关内容五 窗函数法设计FIR滤波器
（FIR滤波器是有限脉冲响应滤波器的简称）

FIR滤波器最大的优点是容易设计成线性相位特性，而且不存在稳定性问题。线性相位特性滤波器在图像处理和数字通信等领域非常有用。
一、 授课目的
熟悉用窗函数法设计FIR数字滤波器
二、授课内容
(一)在MATLAB中产生窗函数的方法
在MATLAB中产生窗函数十分简单
1、矩形窗（Rectangle Window）
调用格式：w=boxcar(n），根据长度n产生一个矩形窗w.
2、三角窗(Triangular Window）
调用格式：w=triang（n)，根据长度n产生一个三角窗w。
3、汉宁窗（HanningWindow）
调用格式：w=hanning(n)，根据长度n产生一个汉宁窗w。
4、海明窗（Hamming Window）
调用格式：w=hamming(n)，根据长度n产生一个海明窗w。
5、布拉克曼窗（Blackman Window）
调用格式:w=blackman(n)，根据长度n产生一个布拉克曼窗w。
6、恺撒窗(Kaiser Window）
调用格式：w=Kaiser(n, β）,根据长度n和影响窗函数旁瓣的β参数产生一个恺撒窗w。
（二）基本窗函数的FIR滤波器的设计—-—-—标准频率响应形状
利用MATLAB提供的函数fir1来实现。
调用格式：fir1（n,Wn，’ftype'，window），n为阶数、Wn是截止频率（如果输入是形如［W1 W2］的矢量时，本函数将设计带通滤波器,其通带为W1 〈ω>W2）、ftype是滤波器的类型(低通--—省略、高通-——ftype=high、带阻——-ftype=stop）、window是窗函数,默认为海明窗。

实验5-1：设计一个长度为8、截止频率为0。4Hz的线性相位FIR滤波器（利用矩形窗)
window=boxcar(8);
b=fir1(7,0.4,window);
freqz（b，1）;
运行得幅频相频特性如下:

实验5—2: 设计一个长度为8、截止频率为0。4Hz的线性相位FIR滤波器（利用布拉克曼窗）
window=blackman（8);
b=fir1（7,0.4,window)；
freqz(b，1）；
运行得幅频相频特性如下:

实验5—3: 设计线性相位FIR带通滤波器，其长度N=15，上下边带截止频率分别为W1=0.3Π，W2=0。5Π
window=blackman(16);
b=fir1(15,[0.3 0。5],window）；
freqz（b，1)；
运行得幅频相频特性如下：