MATLAB下的数字信号处理实现示例10号江城xin要点.docx
《MATLAB下的数字信号处理实现示例10号江城xin要点.docx》由会员分享，可在线阅读，更多相关《MATLAB下的数字信号处理实现示例10号江城xin要点.docx（20页珍藏版）》请在冰豆网上搜索。
MATLAB下的数字信号处理实现示例10号江城xin要点
MATLAB下的数字信号处理实现示例
附录一信号、系统和系统响应
1、理想采样信号序列
（1）首先产生信号x（n）,0<=n<=50
n=0:
50;%定义序列的长度是50
A=444.128;%设置信号有关的参数
a=50*sqrt（2.0）*pi;
T=0.001;%采样率
w0=50*sqrt（2.0）*pi;
x=A*exp（-a*n*T）.*sin（w0*n*T）;%pi是MATLAB定义的π，信号乘可采用“.*”
closeall%清除已经绘制的x（n）图形
subplot（3,1,1）;stem（x）;%绘制x（n）的图形
title（‘理想采样信号序列’）;
（2）绘制信号x（n）的幅度谱和相位谱
k=-25:
25;
W=（pi/12.5）*k;
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）改变参数为：
T=1,Ω=2.0734,a=0.4,A=1
n=0:
50;%定义序列的长度是50
A=1;%设置信号有关的参数
a=0.4;
T=1;%采样率
w0=2.0734;
x=A*exp（-a*n*T）.*sin（w0*n*T）;%pi是MATLAB定义的π，信号乘可采用“.*”
closeall%清除已经绘制的x（n）图形
subplot（3,1,1）;stem（x）;%绘制x（n）的图形
title（‘理想采样信号序列’）;
k=-25:
25;
W=（pi/12.5）*k;
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;
closeall;
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（‘单位冲击信号的相位谱’）
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;
closeall;
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中。
提供了卷积函数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;
closeall;
subplot（3,1,1）;stem（hb）;title（‘系统hb[n]’）;
m=1:
50;%定义序列的长度是50
A=444.128;%设置信号有关的参数
a=50*sqrt（2.0）*pi;
T=0.001;%采样率
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]’）;
6、卷积定律验证
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（‘输出信号的相位谱’）
%以下将验证的结果显示
XHb=X.*Hb;
figure,Subplot（2,1,1）;stem（abs（XHb））;title（'x（n）的幅度谱与hb（n）幅度谱相乘'）;
Subplot（2,1,2）;stem（abs（Y））;title（'y（n）的幅度谱'）;axis（[0,60,0,8000]）
%连续信号的傅里叶变化，在matlab中实现命令为fft（），求连续信号的福谱图则用abs（fft（x））。
n=0:
15;%定义序列的长度是15
p=8;q=2;
x=exp（-1*（n-p）.^2/q）;
closeall;
subplot（3,1,1）;
stem（abs（fft（x）））
p=8;q=4;
x=exp（-1*（n-p）.^2/q）;
subplot（3,1,2）;
stem（abs（fft（x）））
p=8;q=8;
x=exp（-1*（n-p）.^2/q）;
subplot（3,1,3）;
stem（abs（fft（x）））
n=0:
15;%定义序列的长度是15
a=0.1;f=0.0625;
x=exp（-a*n）.*sin（2*pi*f*n）;
closeall;
subplot（2,1,1）;
stem（x）;
subplot（2,1,2）;
stem（abs（fft（x）））
fori=0:
3
x（i）=i+1;x（i+4）=8-（i+4）;
end
fori=8:
15
x（i）=0;
end
closeall;
subplot（2,1,1）;
stem（x）;
subplot（2,1,2）;
stem（abs（fft（x,16）））%fft（x,16）求关于x的16个点的傅立叶变换
附录三窗函数法设计FIR滤波器
一、在MATLAB中产生窗函数十分简单：
（1）矩形窗（RectangleWindow）
调用格式：
w=boxcar（n），根据长度n产生一个矩形窗w。
（2）三角窗（TriangularWindow）
调用格式：
w=triang（n），根据长度n产生一个三角窗w。
（3）汉宁窗（HanningWindow）
调用格式：
w=hanning（n），根据长度n产生一个汉宁窗w。
（4）海明窗（HammingWindow）
调用格式：
w=hamming（n），根据长度n产生一个海明窗w。
（5）布拉克曼窗（BlackmanWindow）
调用格式：
w=blackman（n），根据长度n产生一个布拉克曼窗w。
（6）恺撒窗（KaiserWindow）
调用格式：
w=kaiser（n,beta），根据长度n和影响窗函数旁瓣的β参数产生一个恺撒窗w。
二、基于窗函数的FIR滤波器设计
利用MATLAB提供的函数firl来实现
调用格式：
firl（n,Wn,’ftype’,Window），n为阶数、Wn是截止频率，取值范围为0[例]设计一个长度为8的线性相位FIR滤波器。
Window=boxcar（8）;%产生长度为8的矩形窗
b=fir1（7,0.4,Window）;
freqz（b,1）
%freqz（）计算滤波器的频率响应，此处freqz（b,1）等同于freqz（b）,只是看作带分母且系数为1
Window=blackman（8）;%加布拉克曼窗
b=fir1（7,0.4,Window）;
freqz（b,1）
[例]设计线性相位带通滤波器，其长度N=15，上下边带截止频率分别为W1=0.3π，w2=0.5π
Window=blackman（16）;
b=fir1（15,[0.30.5],Window）;
freqz（b,1）
设计指标为：
ωp=0.2πRp=0.25dBωa=0.3πAs=50dB
的低通数字FIR滤波器
wp=0.2*pi;ws=0.3*pi;
tr_width=ws-wp;
M=ceil（6.6*pi/tr_width）+1;
N=[0:
1:
M-1];
wc=（ws+wp）/2;
hd=ideal_lp（wc,M）;
w_ham=（boxcar（M））’;
h=hd.*w_ham;
[db,mag,pha,grd,w]=freqz_m（h,[1]）;
delta_w=2*pi/1000;
Rp=-（min（db（1:
1:
wp/delta_w+1）））;
As=-round（max（db（ws/delta_w+1:
1:
501）））;
Closeall;
subplot（2,2,1）;stem（hd）;title（‘理想冲击响应’）
axis（[0M-1–0.10.3]）;ylabel（‘hd[n]’）;
subplot（2,2,2）;stem（w_ham）;title（‘汉明窗’）;
axis（[0M-101.1]）;ylabel（‘w[n]’）;
subplot（2,2,3）;stem（h）;title（‘实际冲击响应’）;
axis（[0M-1–0.10.3]）;ylabel（‘h[n]’）;
subplot（2,2,4）;plot（w/pi,db）;title（‘衰减幅度’）;
axis（[01-10010]）;ylabel（‘Decibles’）;
附录四IIR滤波器的实现
MATLAB中滤波器的分析和实现
1、freqs函数：
模拟滤波器的频率响应
[例]系统传递函数为
的模拟滤波器，在MATLAB中可以用以下程序来实现：
a=[10.41];
b=[0.20.31];
w=logspace（-1,1）;%产生从10^-1到10^1之间的50个等间距点，即50个频率点
freqs（b,a,w）%根据输入的参数绘制幅度谱和相位谱
2、freqz函数：
数字滤波器的频率响应
程序来实现：
a=[10.41];
b=[0.20.31];
%根据输入的参数绘制幅度谱和相位谱，得到0到π之间128个点处的频率响应
freqz（b,a,128）
3、ButterWorth模拟和数字滤波器
（1）butterd函数：
ButterWorth滤波器阶数的选择。
调用格式：
[n,Wn]=butterd（Wp,Ws,Rp,Rs），在给定滤波器性能的情况下（通带临界频率Wp、阻带临界频率Ws、通带内最大衰减Rp和阻带内最小衰减Rs），计算ButterWorth滤波器的阶数n和截止频率Wn。
相同参数条件下的模拟滤波器则调用格式为：
[n,Wn]=butterd（Wp,Ws,Rp,Rs,’s’）
（2）butter函数：
ButterWorth滤波器设计。
调用格式：
[b,a]=butter（n,Wn），根据阶数n和截止频率Wn计算ButterWorth滤波器分子分母系数（b为分子系数的矢量形式，a为分母系数的矢量形式）。
相同参数条件下的模拟滤波器则调用格式为：
[b,a]=butter（n,Wn,’s’）
[例]采样频率为1Hz，通带临界频率fp=0.2Hz，通带内衰减小于1dB（αp=1）；阻带临界频率fs=0.3Hz，阻带内衰减大于25dB（αs=25）。
设计一个数字滤波器满足以上参数。
[n,Wn]=buttord（0.2,0.3,1,25）;
[b,a]=butter（n,Wn）;
freqz（b,a,512,1）;
4、Chebyshev模拟和数字滤波器
（1）cheb1ord函数：
ChebyshevⅠ型Ⅱ滤波器阶数计算。
调用格式：
[n,Wn]=cheb1ord（Wp,Ws,Rp,Rs），在给定滤波器性能的情况下（通带临界频率Wp、阻带临界频率Ws、通带内波纹Rp和阻带内衰减Rs），选择ChebyshevⅠ型滤波器的最小阶n和截止频率Wn。
（2）cheby1函数：
ChebyshevⅠ型滤波器设计。
调用格式：
[b,a]=butter（n,Rp,Wn），根据阶数n、通带内波纹Rp和截止频率Wn计算ButterWorth滤波器分子分母系数（b为分子系数的矢量形式，a为分母系数的矢量形式）。
注：
ChebyshevⅡ型滤波器所用函数和Ⅰ型类似，分别是cheb2ord、cheby2。
[例]实现上例中的滤波器
[n,Wn]=cheb1ord（0.2,0.3,1,25）;
[b,a]=cheby1（n,1,Wn）;
freqz（b,a,512,1）;
（1）脉冲响应不变法设计数字ButterWorth滤波器
调用格式：
[bz,az]=impinvar（b,a,Fs），再给定模拟滤波器参数b，a和取样频率Fs的前提下，计算数字滤波器的参数。
两者的冲激响应不变，即模拟滤波器的冲激响应按Fs取样后等同于数字滤波器的冲激响应。
（2）利用双线性变换法设计数字ButterWorth滤波器
调用格式：
[bz,az]=bilinear[b,a,Fs]，根据给定的分子b、分母系数a和取样频率Fs，根据双线性变换将模拟滤波器变换成离散滤波器，具有分子系数向量bz和分母系数向量az。
模拟域的butter函数说明与数字域的函数说明相同
[b,a]=butter（n,Wn,’s’）可以得到模拟域的Butterworth滤波器
[例]采样频率为1Hz，通带临界频率fp=0.2Hz，通带内衰减小于1dB（αp=1）；阻带临界频率fs=0.3Hz，阻带内衰减大于25dB（αs=25）。
设计一个数字滤波器满足以上参数。
%直接设计数字滤波器
[n,Wn]=buttord（0.2,0.3,1,25）;
[b,a]=butter（n,Wn）;
freqz（b,a,512,1）;
%脉冲响应不变法设计数字滤波器
[n,Wn]=buttord（0.2,0.3,1,25,’s’）;
[b,a]=butter（n,Wn,’s’）;
freqs（b,a）
[bz,az]=impinvar（b,a,1）;
freqz（bz,az,512,1）
%双线性变换法设计ButterWorth数字滤波器
[n,Wn]=buttord（0.2,0.3,1,25,’s’）;
[b,a]=butter（n,Wn,’s’）;
freqs（b,a）
[bz,az]=bilinear（b,a,1）;
freqz（bz,az,512,1）
升级会员 