关于几个常见数字滤波器的设计.docx
《关于几个常见数字滤波器的设计.docx》由会员分享，可在线阅读，更多相关《关于几个常见数字滤波器的设计.docx（63页珍藏版）》请在冰点文库上搜索。
关于几个常见数字滤波器的设计
实验一信号、系统及系统响应
1、编写常用信号的程序：
（１）：
单位采样序列δ（ｎ）。
程序如下：
function[x,n]=impseq（n0,ns,nf）
n0=0;
ns=-5;
nf=5;
n=[ns:
nf];
x=[（n-n0）==0]
stem（n,x）
grid
title（'单位采样序列'）
xlabel（'n'）
ylabel（'x（n）'）
运行结果如下：
（2）：
单位阶跃序列ｕ（ｎ）。
程序如下：
function[x,n]=stepseq（n0,ns,nf）
n0=0;
ns=-5;
nf=5;
n=[ns:
nf];
x=[（n-n0）>=0]
stem（n,x）
grid
title（'单位阶跃序列'）
xlabel（'n'）
ylabel（'u（n）'）
运行结果如下：
（３）：
单位矩形序列
。
程序如下：
function[x,n]=juxingseq（n0,ns,nf）
n0=4;
ns=0;
nf=9;
n=[ns:
nf];
x=[（（n-ns）>=0）&（（n-nf）<=0）]
stem（n,x）
grid
title（'单位矩形序列'）
xlabel（'n'）
ylabel（'R（n）'）
end
运行结果如下：
（4）正弦序列。
程序如下：
function[x,n]=sin（w,n）
w=0.2*pi;
n=（0:
0.2:
10）;
x=sin（w*n）;
stem（n,x）
title（'正弦序列'）
xlabel（'n'）
ylabel（'x（n）'）
运行结果如下：
（5）：
。
程序如下：
functionx=exp（n）
n=（0:
0.0001:
0.01）;
x=exp（-1000*abs（n））
stem（n,x,'.'）
运行结果如下：
（6）：
。
程序如下：
functionx=rs（n）
n=（0:
0.1:
10）;
A=1;c=2;w=2*pi;
x=A*exp（-c*n）.*sin（w*n）
stem（n,x,'.'）
运行结果如下：
2、认真复习采样理论、离散信号与系统、线性卷积、序列的傅里叶变换及性质等有关内容，阅读本实验原理与方法。
（1）观察并分析采用不同频率时，对函数
的频谱影响。
（ａ）：
以
，对其进行采样得到
。
（ｂ）：
以
，对其进行采样得到
。
（ｃ）：
对
、
采用理想内插函数重建原始信号
。
程序如下：
functionx=exp（n）
n1=（0:
0.0002:
0.01）;
x=exp（-1000*abs（n1））
subplot（2,1,1）
stem（n1,x,'.'）
gridon
n2=（0:
0.001:
0.01）;
x=exp（-1000*abs（n2））
subplot（2,1,2）
stem（n2,x,'.'）
gridon
运行结果如下：
（2）实现两序列的卷积。
利用ｙ＝ｃｏｎｖ（ｘ，ｈ）．例如：
求解
试求：
．
程序如下：
x=[5,9,3,6,-8];
h=[18,7,5,20,11,14,9];
y=conv（x,h）
n=-4:
6;
stem（n,y,'b'）
holdon
gridon
运行结果如下：
（3）系统单位脉冲响应序列产生子程序。
本实验要用到两种FIR系统。
a.ha（n）=R10（n）;
b.hb（n）=δ（n）+2.5δ（n-1）+2.5δ（n-2）+δ（n-3）
hb（n）信号的程序如下：
function[x,n]=impseq（n0,ns,nf）
n0=0;
ns=-5;
nf=5;
n=[ns:
nf];
x0=[（n-n0）==0]
x1=2.5*[（n-1）==0]
x2=2.5*[（n-2）==0]
x3=[（n-3）==0]
x=x0+x1+x2+x3
stem（n,x）
grid
title（'单位采样序列'）
xlabel（'n'）
ylabel（'x（n）'）
运行结果如下：
3.时域离散信号、系统和系统响应分析。
观察信号xb（n）和系统hb（n）的时域和频域特性；利用线性卷积求信号xb（n）通过系统hb（n）的响应y（n），比较所求响应y（n）和hb（n）的时域及频域特性，注意它们之间有无差别，绘图说明，并用所学理论解释所得结果。
b.观察系统ha（n）对信号xc（n）的响应特性。
程序如下：
w=（-3*pi:
0.001:
3*pi）+eps;
Xa=1;
Xb=1+2.5*exp（-j*w）+2.5*exp（-2*j*w）+exp（-3*j*w）;
subplot（2,1,1）
plot（w/pi,abs（Xa）,'b'）
xlabel（'w/pi'）
ylabel（'|Xa|'）
title（'Xa频域特性'）
holdon
gridon
subplot（2,1,2）
plot（w/pi,abs（Xb））
xlabel（'w/pi'）
ylabel（'|Xb|'）
title（'Xb频域特性'）
holdon
gridon
运行结果如下：
程序如下：
function[x,n]=impseq（n0,ns,nf）
n0=0;
ns=0;
nf=3;
n1=[ns:
nf];
x0=[（n1-n0）==0];
xb=[（n1-n0）==0];
x1=2.5*[（n1-1）==0];
x2=2.5*[（n1-2）==0];
x3=[（n1-3）==0];
hb=x0+x1+x2+x3;
y=conv（xb,hb）
n=0:
6;
stem（n,y,'b'）
grid
title（'线性卷积'）
xlabel（'n'）
ylabel（'y（n）'）
运行结果如下：
通过观察响应的y（n）和hb（n）的时域及频域特性，发现两者的时域及频域特性完全相同，说明hb（n）和单位抽样响应相卷积得到的还是hb（n）本身。
用数学表达式为y（n）=hb（n）*δ（ｎ）=hb（n）。
b.程序如下：
function[x,n]=juxingseq（n0,ns,nf）
n0=4;
ns=0;
nf=9;
n=[ns:
nf];
x=[（（n-ns）>=0）&（（n-nf）<=0）]
y=conv（x,x）
n=0:
18;
stem（n,y）
运行结果如下：
系统ha（n）对信号xc（n）的响应特性：
对时域进行采样，频域发生周期延拓。
思考题
在分析理想采样序列特性的实验中，采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同?
它们所对应的模拟频率是否相同?
为什么?
答：
采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量相同，
实验二用FFT作谱分析
1、复习DFT的定义、性质和用DFT作频谱分析的有关内容。
用matlab函数dfs（dft）和idfs（idft）实现dfs（dft）正反变换。
（1）function[Xk]=dfs（xn,N）
function[Xk]=dfs（xn,N）
n=0:
N-1;
k=0:
N-1;
WN=exp（-j*2*pi/N）;
nk=n'*k;
Xk=xn*WN.^nk;
（2）function[xn]=idfs（Xk,N）
n=0:
N-1;
k=0:
N-1;
WN=exp（-j*2*pi/N）;
nk=-n'*k;
Xk=xn*WN.^nk;
2、计算周期序列
的DFS.
N=16;
xn=[4567456745674567];
n=0:
N-1;
Xk=dfs（xn,N）;
stem（n,abs（Xk））
ylabel（'X（k）'）
gridon
3、为了说明高密度频谱和高分辨率频谱之间的区别，考察序列x（n）有限个样本的频谱。
当
时，求x（n）的DFT.
当
时，求x（n）的DFT.
程序如下：
N=11;
xn=cos（0.48*pi*n）+cos（0.52*pi*n）;
n=0:
N-1;
Xk=dfs（xn,N）;
stem（n,abs（Xk））
xlabel（'n'）;
ylabel（'X（k）'）
title（'N=11'）
gridon
运行结果如下：
N=101;
xn=cos（0.48*pi*n）+cos（0.52*pi*n）;
n=0:
N-1;
Xk=dfs（xn,N）;
stem（n,abs（Xk）,'.'）
xlabel（'n'）;
ylabel（'X（k）'）
title（'N=101'）
gridon
4、利用矩阵相乘计算循环卷积。
设
计算其5点、6点、7点、8点、9点循环卷积。
并分析线性卷积与循环卷积的关系。
（1）自定义圆周移位函数cirshift.m
functiony=cirshift（x,m,N）
iflength（x）>N
error（'N必须>=x的长度'）
end
x=[xzeros（1,N-length（x））];
n=[0:
1:
N-1];
n=mod（n-m,N）;
y=x（n+1）;
（2）自定义圆周卷积函数
functiony=circonvt（x1,x2,N）
%checkforlengthofx1
iflength（x1）>N
error（'N必须>=x1的长度'）
end
%checkforlengthofx2
iflength（x2）>N
error（'N必须>=x2的长度'）
end
x1=[x1zeros（1,N-length（x1））];
x2=[x2zeros（1,N-length（x2））];
m=[0:
1:
N-1];
x2=x2（mod（-m,N）+1）;
H=zeros（N,N）;
forn=1:
1:
N
H（n,:
）=cirshift（x2,n-1,N）;
end
y=x1*H'
（3）命令窗口输入：
①5点圆周卷积
x1=[12345];
x2=[6789];
y=circonvt（x1,x2,5）
输出
y=
100958570100
②6点圆周卷积
x1=[12345];
x2=[6789];
y=circonvt（x1,x2,6）
输出：
y=
8264407010094
③7点圆周卷积
x1=[12345];
x2=[6789];
y=circonvt（x1,x2,7）
输出：
y=
511940701009476
④8点圆周卷积
x1=[12345];
x2=[6789];
y=circonvt（x1,x2,8）
输出：
y=
6194070100947645
⑤9点圆周卷积
x1=[12345];
x2=[6789];
y=circonvt（x1,x2,9）
输出：
y=
61940701009476450
6、编制信号产生子程序，产生以下典型信号供谱分析用：
（1）x1（n）=R4（n）
function[x,n]=juxingseq（n0,ns,nf）
n0=1.5;
ns=0;
nf=3;
n=[ns:
nf];
x=[（（n-ns）>=0）&（（n-nf）<=0）];
stem（n,x）;
xlabel（'采样时间n'）;
ylabel（'单位矩形序列R4（n）'）;
title（'单位矩形序列'）
gridon
（2）
X2（n）=
function[x,n]=fenduan（n0,n1,n2,n3）
n0=0;n1=3;n2=4;n3=7;
n=[n0:
n1];
xn=n+1;
n_0=[n2:
n3];
xn_0=8-n_0;
stem（n,xn）
holdon
stem（n_0,xn_0）
xlabel（'采样时间n'）;
ylabel（'x2（n）'）;
gridon
（3）
X3（n）=
function[x,n]=fenduan1（n0,n1,n2,n3）
n0=0;n1=3;n2=4;n3=7;
n=[n0:
n1];
xn=4-n;
n_0=[n2:
n3];
xn_0=n_0-3;
stem（n,xn）
holdon
stem（n_0,xn_0）
xlabel（'采样时间n'）;
ylabel（'x3（n）'）;
gridon
（4）x4（n）=cosn
function[x,n]=cos（w,n）
w=pi/4;
n=（0:
0.4:
8）;
x=sin（w*n）;
stem（n,x）
xlabel（'采样时间n'）;
ylabel（'余弦序列x（n）'）;
title（'余弦序列'）
gridon
（5）x5（n）=sinn
function[x,n]=sin（w,n）
w=pi/8;
n=（0:
0.5:
16）;
x=sin（w*n）;
stem（n,x）
xlabel（'采样时间n'）;
ylabel（'正弦序列x（n）'）;
title（'正弦序列'）
gridon
（6）x6（t）=cos8t+cos16+cos
t=0:
0.005:
0.5;
xt=cos（8*pi*t）+cos（16*pi*t）+cos（20*pi*t）;
plot（t,xt）
gridon
xlabel（'时间t'）;
ylabel（'x6（t）'）;
实验内容
（1）对2中所给出的信号逐个进行谱分析。
下面给出针对各信号的FFT变换
区间N以及对连续信号x6（t）的采样频率fs，供实验时参考。
x1（n）,x2（n）,x3（n）,x4（n）,x5（n）:
N=8,
16=x6（t）:
fs=64（Hz）,N=16,32,64
①xn=[（（n-0）>=0）&（（n-3）<=0）];;
N=8;
fs=64;
D=2*pi*fs/N;
n=0:
N-1;
k=floor（-（N-1）/2:
（N-1）/2）;
X=fftshift（fft（xn,N））;
subplot（2,1,1）;
plot（k*D,abs（X））
title（'幅度频谱'）;
gridon
subplot（2,1,2）;
plot（k*D,angle（X））
title（'相位频谱'）;
gridon
②
N=8;
n=0:
N-1;
xn=fenduan（）;
fs=64;
D=2*pi*fs/N;
k=floor（-（N-1）/2:
（N-1）/2）;
X=fftshift（fft（xn,N））;
subplot（2,1,1）;
plot（k*D,abs（X））
title（'幅度频谱'）;
gridon
subplot（2,1,2）;
plot（k*D,angle（X））
title（'相位频谱'）;
gridon
③
N=8;
n=0:
N-1;
xn=fenduan1（）;
fs=64;
D=2*pi*fs/N;
k=floor（-（N-1）/2:
（N-1）/2）;
X=fftshift（fft（xn,N））;
subplot（2,1,1）;
plot（k*D,abs（X））
title（'幅度频谱'）;
gridon
subplot（2,1,2）;
plot（k*D,angle（X））
title（'相位频谱'）;
gridon
④xn=cos（pi/4*n）;
N=8;
fs=64;
D=2*pi*fs/N;
n=0:
N-1;
k=floor（-（N-1）/2:
（N-1）/2）;
X=fftshift（fft（xn,N））;
subplot（2,1,1）;
plot（k*D,abs（X））
title（'幅度频谱'）;
gridon
subplot（2,1,2）;
plot（k*D,angle（X））
title（'相位频谱'）;
gridon
⑤N=8;
n=0:
N-1;
xn=sin（pi/8*n）;
fs=64;
D=2*pi*fs/N;
k=floor（-（N-1）/2:
（N-1）/2）;
X=fftshift（fft（xn,N））;
subplot（2,1,1）;
plot（k*D,abs（X））
title（'幅度频谱'）;
gridon
subplot（2,1,2）;
plot（k*D,angle（X））
title（'相位频谱'）;
gridon
⑥N0=[163264];
forr=1:
3;
N=N0（r）;
n=0:
N-1;
fs=64;T=1/fs;
xt=cos（8*pi*n*T）+cos（16*pi*n*T）+cos（20*pi*n*T）;
D=2*pi*fs/N;
k=floor（-（N-1）/2:
（N-1）/2）;
X=T*fftshift（fft（xn,N））;
subplot（3,1,r）;
plot（k*D,abs（X））
title（'幅度频谱'）;
axis（[-200,200,1.1*min（abs（X））,1.1*max（abs（X））]）;
gridon
end
（2）令x（n）=x4（n）+x5（n），用FFT计算8点和16点离散傅里叶变换，
X（k）=DFT［x（n）］
N1=8;
N2=16;
N1=0:
N1-1;
n2=0:
N2-1;
xn=cos（pi/4*n）+sin（pi/8*n）;
k1=0:
N1-1;
k2=0:
N2-1;
Xk1=fft（xn,N1）;
Xk2=fft（xn,N2）;
subplot（2,1,1）;
stem（k1,abs（Xk1））
title（'N=8'）;
gridon
subplot（2,1,2）;
stem（k2,abs（Xk2））
title（'N=16'）;
gridon
（3）令x（n）=x4（n）+jx5（n），重复
（2）。
N1=8;
N2=16;
n1=0:
N1-1;
n2=0:
N2-1;
xn1=cos（pi/4*n1）+j*sin（pi/8*n1）;
xn2=cos（pi/4*n2）+j*sin（pi/8*n2）;
k1=0:
N1-1;
k2=0:
N2-1;
Xk1=fft（xn1,N1）;
Xk2=fft（xn2,N2）;
subplot（2,1,1）;
stem（k1,abs（Xk1））
title（'N=8'）;
gridon
subplot（2,1,2）;
stem（k2,abs（Xk2））
title（'N=16'）;
gridon
思考题
在N=8时，x2（n）和x3（n）的幅频特性会相同吗?
为什么?
N=16呢?
答：
N=8，x2（n）和x3（n）8点DFT的模相等，故幅频特性会相同。
N=16时，x2（n）与x3（n）不满足循环移位关系，所以x2（n）和x3（n）16点DFT的模不相等，故幅频特性不相同。
实验三设计IIR数字滤波器
1、用Matlab语言分别设计巴特奥斯低通滤波器和切比雪夫低通滤波器，
其技术指标为：
通带截止频率
，通带最大衰减
；
阻带起始频率
，阻带最小衰减
；
要求：
求出他们的零点、极点、阶数、增益等，并画出图形作比较。
程序如下：
Wp=2*pi*5*10^3;Ws=2*pi*10^4;
rp=3;
rs=30;
wp=1;
ws=Ws/Wp;
[N,wc]=buttord（wp,ws,rp,rs,'s'）;
[z,p,k]=buttap（N）;
[B,A]=zp2tf（z,p,k）;
w=0:
0.05*pi:
2*pi;
[h,w]=freqs（B,A,w）;
plot（w,20*log10（abs（h））,'r'）;
holdon
gridon
[N2,wc2]=cheb1ord（wp,ws,rp,rs,'s'）;
[z2,p2,k2]=cheb1ap（N2,rp）;
[B2,A2]=zp2tf（z2,p2,k2）;
[h2,w]=freqs（B2,A2,w）;
plot（w,20*log10（abs（h2））,'b'）;
holdon
gridon
xlabel（'频率/Hz'）
ylabel（'幅值'）
运行结果如下：
零点、极点、阶数、增益如下：
巴特奥斯低通滤波器参数：
N=5
wc=1.0025
z=[]
p=
-0.3090+0.9511i
-0.3090-0.9511i
-0.8090+0.5878i
-0.8090-0.5878i
-1.0000
k=1
切比雪夫低通滤波器参数：
N2=4
wc2=1
z2=[]
p2=
-0.0852+0.9465i
-0.2056+0.3920i
-0.2056-0.3920i
-0.0852-0.9465i
k2=0.1253
2、切比雪夫低通滤波器，其技术指标为：
通带截止频率
，通带最大衰减
；
阻带起始频率
，阻带最小衰减
；
要求：
求出他们的零点、极点、阶数、增益等，并画出图形作比较。
程序如下：
Wp=2*pi*3*10^6;Ws=2*pi*12*10^6;
rp=0.1;
rs=60;
wp=1;
ws=Ws/Wp;
[N,wc]=buttord（wp,ws,rp,rs,'s'）;
[z,p,k]=buttap（N）;
[B,A]=zp2tf（z,p,k）;
w=0:
0.05*pi:
2*pi;
[h,w]=freqs（B,A,w）;
plot（w,20*log10（abs（h））,'r'）;
holdon
gridon
[N2,wc2]=cheb1ord（wp,ws,rp,rs,'s'）;
[z2,p2,k2]=cheb1ap（N2,rp）;
[B2,A2]=zp2tf（z2,p2,k2）;
[h2,w]=freqs（B2,A2,w）;
plot（w,20*log10（abs（h2））,'b'）;
holdon
gridon
xlabel（'频率/Hz'）
ylabel（'幅值'）
运行结果如下：
巴特奥斯
升级会员 