数字信号处理实验Word格式文档下载.docx
《数字信号处理实验Word格式文档下载.docx》由会员分享,可在线阅读,更多相关《数字信号处理实验Word格式文档下载.docx(13页珍藏版)》请在冰豆网上搜索。
![数字信号处理实验Word格式文档下载.docx](https://file1.bdocx.com/fileroot1/2023-1/1/6d2a3bbe-afab-4a3c-8a58-ba8f16cc9431/6d2a3bbe-afab-4a3c-8a58-ba8f16cc94311.gif)
subplot(4,3,4);
plot(n,y1);
直接滤波的结果'
subplot(4,3,5);
plot(n,abs(y1),'
g'
直接滤波幅频特性曲线'
subplot(4,3,6);
plot(n,angle(y1),'
m'
直接滤波相频特性曲线'
subplot(4,3,7);
plot(n,y2,'
cascade滤波后的波形'
subplot(4,3,8);
plot(n,abs(y2),'
r'
cascade滤波幅频特性曲线'
subplot(4,3,9);
plot(n,angle(y2),'
b'
cascade滤波相频特性曲线'
subplot(4,3,10);
plot(n,y3,'
parallel滤波的波形'
subplot(4,3,11);
plot(n,abs(y3),'
parallel滤波幅频特性曲线'
subplot(4,3,12);
plot(n,angle(y3),'
parallel滤波相频特性曲线'
)
调用的功能函数:
%单位取样序列函数
function[x,n]=impseq(n0,n1,n2)
n=[n1:
n2];
x=[(n-n0)==0];
%IIFDF直接型到级联型的转换函数
function[b0,B,A]=dir2cas(b,a)
%计算增益系数
b0=b
(1);
b=b/b0;
a0=a
(1);
a=a/a0;
b0=b0/a0;
%
M=length(b);
N=length(a);
ifN>
M
b=[b,zeros(1,N-M)];
elseifM>
N
a=[a,zeros(1,M-N)];
N=M;
else
NM=0;
end
K=floor(N/2);
B=zeros(K,3);
A=zeros(K,3);
ifK*2==N
b=[b,0];
a=[a,0];
broots=cplxpair(roots(b));
aroots=cplxpair(roots(a));
fori=1:
2:
2*K
Brow=broots(i:
1:
i+1,:
Brow=real(poly(Brow));
B(fix((i+1)/2),:
)=Brow;
Arow=aroots(i:
Arow=real(poly(Arow));
A(fix((i+1)/2),:
)=Arow;
End
functionhcas=casfilter(b0,B,A,delta)
[K,L]=size(B);
N=length(delta);
w=zeros(K+1,N);
w(1,:
)=delta;
K
w(i+1,:
)=filter(B(i,:
),A(i,:
),w(i,:
));
hcas=b0*w(K+1,:
function[C,B,A]=dir2par(b,a)
[r1,p1,C]=residuez(b,a);
p=cplxpair(p1,10000000*eps);
I=cplxcomp(p1,p);
r=r1(I);
B=zeros(K,2);
fori=1:
N-2
Brow=r(i:
Arow=p(i:
[Brow,Arow]=residuez(Brow,Arow,[]);
)=real(Brow);
)=real(Arow);
end
[Brow,Arow]=residuez(r(N-1),p(N-1),[]);
B(K,:
)=[real(Brow)0];
A(K,:
)=[real(Arow)0];
N-1
functionhpar=parfilter(C,B,A,delta)
)=filter(C,1,delta);
),delta);
hpar=sum(w);
程序运行的结果如下:
n=[-20,20];
[x2,n]=impseq(0,-20,20);
y1=filter(b,a,x2);
y2=casfilter(b0,B,A,x2);
y3=parfilter(C,B,A,x2);
plot(n,x2,'
x2的图像'
直接滤波波形'
cascade滤波波形'
plot(n,y3);
parallel滤波图像'
plot(n,abs(y3));
parallel滤波幅频特性'
plot(n,angle(y3));
parallel滤波相频特性'
2.用Matlab设计一个模拟巴特沃斯低通滤波器,指标
,采用冲激响应不变法和双线性变换设计。
冲激响应不变法:
%冲激响应不变法设计数字巴特沃斯低通滤波器
%数字滤波器指标
wp=0.2*pi;
ws=0.3*pi;
Rp=1;
As=15;
%转换为模拟域指标
T=1;
omegap=wp/T;
omegas=ws/T;
%模拟巴特沃斯滤波器的计算
[cs,ds]=afd_buttap(omegap,omegas,Rp,As);
%冲激不变法
[b,a]=imp_invr(cs,ds,T);
[C,B,A]=dir2par(b,a)
%计算数字滤波器的频率响应
[db,mag,pha,grd,w]=freqz_m(b,a);
subplot(2,2,1);
plot(w/pi,mag);
幅度响应'
gridon;
axis([0,0.8,0,1]);
ylabel('
|幅度|'
xlabel('
以\pi为单位的频率'
subplot(2,2,3);
plot(w/pi,db);
幅度响应(dB)'
axis([0,0.8,-80,0]);
对数幅度/dB'
subplot(2,2,2);
plot(w/pi,pha);
相位响应'
axis([0,0.8,-4,4]);
相位'
subplot(2,2,4);
plot(w/pi,grd);
群延迟'
axis([0,0.8,0,10]);
样本'
function[b,a]=afd_buttap(wp,ws,Rp,As)
N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(wp/ws)));
fprintf('
\n***ButterworthFilterOrder=%2.0f\n'
N)
omega=wp/((10^(Rp/10)-1)^(1/(2*N)));
[b,a]=buttap_o(N,omega);
function[b,a]=imp_invr(c,d,T)
[R,p,k]=residue(c,d);
p=exp(p*T);
[b,a]=residuez(R,p,k);
b=real(b'
a=real(a'
[Brow,Arow]=residues(Brow,Arow,[]);
function[db,mag,pha,grd,w]=freqz_m(b,a)
[H,w]=freqz(b,a,1000,'
whole'
H=(H(1:
501))'
;
w=(w(1:
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
grd=grpdelay(b,a,w);
双线性变换:
%双线性变换设计数字巴特沃斯低通滤波器
Fs=1/T;
omegap=(2/T)*tan(wp/2);
omegas=(2/T)*tan(ws/2);
[cs,ds]=afd_buttap(omegap,omegas,Rp,As)
%双线性变换
[b,a]=bilinear(cs,ds,Fs);
[C,B,A]=sdir2cas(b,a)
axis([0,0.8,-60,0]);
%将模拟滤波器直接形式转换为级联形式的函数
function[C,B,A]=sdir2cas(b,a)
%计算增益系数C
C=b0/a0
Na=length(a)-1;
Nb=length(b)-1;
%计算分母的二阶因子部分
p=cplxpair(roots(a));
K=floor(Na/2);
ifK*2==Na
A=zeros(K,3);
forn=1:
Na
Arow=p(n:
n+1,:
Arow=poly(Arow);
A(fix((n+1)/2),:
elseifNa==1
A=[0real(poly(p))];
A=zeros(K+1,3);
A(K+1,:
)=[0real(poly(p(Na)))];
%计算分子的二阶因子部分
z=cplxpair(roots(b));
K=floor(Nb/2);
ifNb==0
B=[00poly(z)];
elseifK*2==Nb
B=zeros(K,3);
Nb
Brow=z(n:
Brow=poly(Brow);
B(fix((n+1)/2),:
elseifNb==1
B=[0real(poly(z))];
B=zeros(K+1,3);
B(K+1,:
)=[0real(poly(z(Nb)))];
3.用级联型和并联型结构实现以下系统函数,每个二阶节都采用直接II型结构。
级联型:
a=[5,-2.4412*5,2.4412*5,-5];
%输入分子系数向量
b=[1,-1.7728,1.4464,-0.405];
%输入分母系数向量
[z,p,k]=tf2zp(a,b);
%求出各子系统的零极点
[sos,A]=zp2sos(z,p,k);
%求出各二阶节乘系数
disp('
零点:
'
disp(z);
极点:
disp(p);
增益系数:
disp(k);
二阶节:
disp(real(sos));
并联型:
[r,p,k]=residuez(a,b);
留数:
disp(r'
disp(p'
常数:
disp(k'