dsp各种程序.docx
《dsp各种程序.docx》由会员分享，可在线阅读，更多相关《dsp各种程序.docx（27页珍藏版）》请在冰豆网上搜索。
dsp各种程序
DSP题库详解
1已知3阶椭圆IIR数字低通滤波器的性能指标为：
通带截止频率0.4π，通带波纹为0.6dB，最小阻带衰减为32dB。
设计一个6阶全通滤波器对其通带的群延时进行均衡。
绘制低通滤波器和级联滤波器的群延时。
%Progranm1
%Group-delayequalizationofanIIRfilter.
%
[n,d]=ellip（3,0.6,32,0.4）;
[GdH,w]=grpdelay（n,d,512）;
plot（w/pi,GdH）;grid
xlabel（'\omega/\pi'）;ylabel（'Groupdelay,samples'）;
title（'OriginalFilter'）;
F=0:
0.001:
0.4;
g=grpdelay（n,d,F,2）;%Equalizethepassband
Gd=max（g）-g;
%Designtheallpassdelayequalizer
[num,den,tau]=iirgrpdelay（6,F,[00.4],Gd）;%设计六阶的全通滤波器
[b,a]=iirgrpdelay（6,F,[00.4],Gd）;
He1=dfilt.df2（b,a）;
He=dfilt.df2（n,d）;
He_all=dfilt.cascade（He,He1）;
grpdelay（He_all）%DFILT:
DigitalFilterImplementation.
%[GdA,w]=grpdelay（num,den,512）;
%figure
（2）;
%plot（w/pi,GdH+GdA）;grid
%xlabel（'\omega/\pi'）;ylabel（'Groupdelay,samples'）;
%title（'GroupDelayEqualizedFilter'）;
2设计巴特沃兹模拟低通滤波器，其滤波器的阶数和3-dB截止频率由键盘输入，程序能根据输入的参数，绘制滤波器的增益响应。
%Program2
%ProgramtoDesignButterworthAnalogLowpassFilter
%
%Typeinthefilterorderandpassbandedgefrequency设计边缘频率
N=input（'Typeinfilterorder='）;
Wn=input（'3-dBcutoffangularfrequency='）;
%Determinethetransferfunction
[num,den]=butter（N,Wn,'s'）;%计算滤波器的分子和分母
%Computeandplotthefrequencyresponse
[h,w]=freqz（num,den）;
plot（w/pi,20*log10（abs（h）））;
xlabel（'\omega/\pi'）;ylabel（'Magnitude,dB'）;
title（'Themagnituderespponseofthesystem'）
axis（[03-805]）
3%%%%%%%%%%
%Tousetheresidue
%
num=[1,-0.2,0.5,0,0];
den=[1,3.2,1.5,-0.8,1.4];
[r,p,k]=residue（num,den）
4%Program4
%Cheby1analoghighpassfilterdesign
%
Wp=0.75;Ws=0.35;Rs=43;Rp=0.5;
[N,Wn]=cheb1ord（Wp,Ws,Rp,Rs,'s'）;
%Determinethecoefficientsofthetransferfunction
[num,den]=cheby1（N,Rp,Wn,'high','s'）;
%Computeandplotthefrequencyresponse
omega=[0:
1:
30*pi];%w=0:
0.01:
pi;
h=freqs（num,den,omega）;%plot（w/pi,20*log10（abs（h）））
plot（omega/（2*pi）,20*log10（abs（h）））;
xlabel（'\omega/\pi'）;ylabel（'Magnitude,dB'）;
title（'IIRChebyshev1HighpassFilter/digital'）
5.
%Programto绘制30点序列的实部与虚部
%
n=0:
29;
x=0.2*exp（（0.4+0.5*j）.*n）;
subplot（1,2,1）
stem（n,real（x））
title（'Therealpart'）;
xlabel（'n'）;ylabel（'magnitude'）
subplot（1,2,2）
stem（n,imag（x））
title（'Theimaginepart'）
xlabel（'n'）;ylabel（'magnitude'）
6.
%ProgramtoDesignType1ChebyshevanalogLowpassFilter
%
%Readinthefilterorder,passbandedgefrequency
%andpassbandrippleindB
N=input（'Order='）;
Wc=input（'typeinthe3dBcutofffrequency='）;
Rp=input（'Typeinthepeakripplevalueinpassband:
'）;
%Determinethecoefficientsofthetransferfunction
[num,den]=cheby1（N,Rp,Wc,'s'）;
%Computeandplotthefrequencyresponse
[h,w]=freqs（num,den）;
plot（w/pi,20*log10（abs（h）））;
xlabel（'Frequency,Hz'）;ylabel（'magnitude,dB'）;
title（'IIRChebyshev1LowpassFilter/analog'）
gridon
7.
%DeterminationoftheFactoredForm
%ofaRationalz-Transform
%熟悉命令[num,den]=residuez（r,p,k）；[r,p,k]=residuez（num,den）
clearall
r=[10.61.8];
p=[-3.22.42.4];
k=0.2;
[num,den]=residuez（r,p,k）;
disp（'Numeratorpolynomialcoefficients'）;disp（num）
disp（'Denominatorpolynomialcoefficients'）;disp（den）
8,设计IIR数字带通滤波器
%ProgramtoDesigniirdigitalbandpassFilter
clearall
Wp=[0.40.6];%WpandWs都是向量组，带通的性质
Ws=[0.30.7];
Rp=0.6;
Rs=35;
[N,Wn]=buttord（Wp,Ws,Rp,Rs）;
[b,a]=butter（N,Wn）;%巴特卧兹滤波器的分子和分母
[h,w]=freqz（b,a,256）;%求出在频域的画图点，256，即要花的点数
gain=20*log10（abs（h））;
plot（w/pi,gain）;gridon
xlabel（'\omega/\pi'）;ylabel（'magnitude,dB'）;
title（'IIRButterworthBandpassFilter'）;
9%drawtherequestedsequence
clear
n=0:
23;
x=2*0.9.^n;
stem（n,x）
title（'x[n]'）
xlabel（'n'）;ylabel（'magnitude'）
gridon
10%
%ProgramtoDesignEllipticanalogLowpassFilter%椭圆模拟低通滤波器设计
%Readinthefilterorder,passbandedgefrequency,
%passbandrippleindBandminimumstopband
%attenuationindB
N=input（'Order='）;
Wn=input（'typeinthe3dBcutofffrequency='）;
Rp=input（'typeinthepeakripplevalueofpassbandindB='）;
Rs=input（'typeintheMinimumstopbandattenuationindB='）;
[num,den]=ellip（N,Rp,Rs,Wn,'s'）;%Determinethecoefficientsofthetransferfunction
%Computeandplotthefrequencyresponse
[h,w]=freqs（num,den）;
plot（w/2*pi,20*log10（abs（h）））;
xlabel（'Frequency,Hz'）;ylabel（'Gain,dB'）;
title（‘iirelliplowpassfilter’）
girdon
%%%%%
num=[1,-0.2,0.5];
den=[1,3.2,1.5,-0.8,1.4];
y=impz（num,den,30）;%inverse–transform
n=1:
30;
stem（n,y）;gridon
%%%%%%%%%%%
clear
N=5;
Rp=0.8;
Rs=35;
Wn=0.35;%去掉了pi的输入
[num,den]=ellip（N,Rp,Rs,Wn）;
F=0:
0.01:
0.3;
[GdH,w]=grpdelay（num,den,512）;%gd=grpdelay（num,den,w,2）;群时延的另外一个公式
figure
（1）,plot（w/pi,GdH）;grid
xlabel（'\omega/\pi'）;ylabel（'Groupdelay,samples'）;
title（'ThegrpdelayofIIREllipLowpassFilter'）
g=grpdelay（num,den,F,2）;%Equalizethepassband
Gd=max（g）-g;
%Designtheallpassdelayequalizer
[num,den]=iirgrpdelay（10,F,[00.3],Gd）;
He=dfilt.df2（num,den）;
He1=dfilt.df2（num,den）;
He_all=dfilt.cascade（He,He1）;
grpdelay（He_all）%下面是另外一种方法
%[GdA,w]=grpdelay（num,den,512）;
%plot（w/pi,GdH+GdA）;gridon
%xlabel（'\omega/\pi'）;ylabel（'Groupdelay,samples'）;
%title（'GroupDelayEqualizedFilter'）;
13%Programtodrawsignalandsignalwithnoise
clear
N=40;
n=0:
N-1;
s=3*（n.*（0.8.^n））;
d=1.2.*rand（1,N）-0.6;
x=s+d;
figure
（1）
stem（n,s）
xlabel（'timeindexn'）;ylabel（'amplitude'）;title（'originalsignals[n]'）
figure
（2）
stem（n,d）
title（'noised[n]'）
xlabel（'timeindexn'）;ylabel（'amplitude'）
figure（3）
stem（n,x）
title（'signalwithnoisex[n]'）
xlabel（'timeindexn'）;ylabel（'amplitude'）
%
%
%generatemovingaveragefilter
b=ones（4,1）/4;%windowlengthis4
y=filter（b,1,x）;%letsignalgothroughthefilter
figure（4）
stem（n,y）
title（'signalthroughfiltery[n]'）
xlabel（'timeindexn'）;ylabel（'amplitude'）
figure（5）
plot（n,d,'g:
',n,s,'b',n,x,'c--',n,y,'r-.'）;
legend（'noise','s[n]','x[n]','y[n]'）%在图上加方框对其标注
gridon
14%%%%%%
N=10;
x=ones（1,10）;
[h,w]=freqz（x,1,16）;
subplot（2,1,1）
plot（w,abs（h））
xlabel（'\omega/'）;ylabel（'amplitude'）;title（'amplituderesponse'）
subplot（2,1,2）
plot（w,angle（h））
xlabel（'\omega/'）;ylabel（'phase（rad）'）;title（'phaseresponse'）
%notice:
hereomegahasn’tbeen归一化
15%
%已知系统的系统函数为：
%用MATLAB的filter函数求h[n]的前20个样本
clear
N=20;
num=[1,-0.2,0.5];
den=[1,3.2,1.5,-0.8,1.4];
x=zeros（1,N）;
x
（1）=1;%x=[1,zeros（1,19）]
h=filter（num,den,x）;
y=impz（num,den,N）;%用filter和z反变换impz对比
subplot（2,1,1）
stem（[0:
N-1],h）
title（'h[n]'）
subplot（2,1,2）
stem（[0:
N-1],y）
title（'y[n]'）
16%
%利用Hermann公式估计FIR低通滤波器的阶数。
该滤波器的性能指标为：
通带截止频率为1500Hz，阻带截至频率为1800Hz，通带纹波为
=0.015，阻带纹波为
=0.021，抽样频率为5000Hz
%
fedge=[15001800];%unitHz
mval=[10];%desiredmagnitudevaluesineachband
dev=[0.0150.021];%通带和阻带波纹
Ft=5000;
N=remezord（fedge,mval,dev,Ft）
%第二种方法
%clearall
%formatlong
%a1=0.005309;a2=0.07114;a3=-0.4761;a4=0.00266;a5=0.5941;a6=0.4278;
%b1=11.01217;b2=0.51244;
%f=[1500,1800];
%w=2*pi*f/5000;
%wp=w
（1）;ws=w
（2）;
%delta_p=0.015;delta_s=0.021;
%delta=[0.015,0.021];
%ifdelta_p>=delta_s
%m=1;n=2;
%D=（a1*（log10（delta（m）））^2+a2*log10（delta（m））+a3）*log10（delta（n））-（a4*（log10（delta（m）））^2+a%5*log10（delta（m））+a6）;
%F=b1+b2*（log10（delta（m））-log10（delta（n）））;
%else
%m=2;n=1;
%D=（a1*（log10（delta（m）））^2+a2*log10（delta（m））+a3）*log10（delta（n））-（a4*（log10（delta（m）））^2+a%5*log10（delta（m））+a6）;
%F=b1+b2*（log10（delta（m））-log10（delta（n）））;
%end
%N=ceil（（D-F*（（ws-wp）/2/pi）^2）/（（ws-wp）/2/pi））
（anotherversionof13,14,15,16）
13%编写四点滑动平均滤波器。
原始未受干扰的序列为：
，加性噪声信
为随机序列，幅度0.6，受干扰的序列为：
，分别绘制长度为40的原始序列，噪声序列和受干扰序列，以及滑动平均滤波器的输出。
%
R=40;
d=1.2*rand（1,R）-0.6;
n=0:
R-1;
x=3*n.*（0.8）.^n;
s=x+d;
plot（n,d,'bla',n,x,'r',n,s）;
legend（'d[n]','x[n]','s[n]'）
m=4;%input（'pleaseinputthenumberofpointsyouwantfilt='）;
m=（0:
m-1）/m;
y=filter（m,1,x）;
figure
plot（n,s,n,y,'g'）
legend（'s[n]','y[n]'）
14%绘制长度为10点的矩形序列的16点离散傅里叶变换样本的幅度和相位
%
a=input（'pleaseinputthemagnitudeoftheN-pointsequence='）;
N=input（'pleaseinputthelenthofsequence='）;
M=input（'pleaseinputthelenthoftheDFTsequence='）;
u=a*ones（1,N）;
h=fft（u,M）;
m=0:
M-1;
subplot（2,1,1）;
stem（m,abs（h））,title（'magnitude'）;
subplot（2,1,2）;
stem（m,angle（h））,title（'phase'）;
15%
%
N=input（'pleaseinputthelenthofthewantedsequence='）;
n=ones（1,N）;
a=input（'请按降幂输入分子系数='）;
b=input（'请按降幂输入分母系数='）;
y=filter（a,b,n）
16%
%
fedge=input（'Typeinthebandedges='）;
mval=input（'Desiredmagnitudevaluesineachband='）;
dev=input（'Allowabledeviationineachband='）;
FT=input（'Typeinthesamplingfrequency='）;
[N,fpts,mag,wt]=remezord（fedge,mval,dev,FT）;
fprintf（'Filterorderis%d\n',N）
17%编写长度为5的中值滤波器程序，原始未受干扰的序列为：
s[n]=3[n（0.8）^n],
加性噪声信号d【n】为随机序列，幅度为0.6，分别绘制长度为40的受干扰序列
以及中值滤波器的输出
clear
n=1:
40;
s=3*n.*power（0.8,n）;%power表示指数函数
figure
（1）
plot（n,s）;
d=1.2*（rand（1,40）-0.5）;
x=s+d;
figure
（2）
plot（n,x）
xx=medfilt1（x,5）%中值滤波器设计
figure（3）,plot（n,x,'g',n,xx,'b'）
legend（'[n]','y[n]'）
18%已知16点序列x[n]的DFT为：
X[k]=k/160<=k<=15X[k]=0else绘制x[n]的实部和虚部
%
k=0:
15;
X=k/15;
x=ifft（X）;%DFT的反变换IFFT
figure
（1）;
stem（k,real（x））;%实部图
title（'realpart'）;
figure
（2）;
title（'imaginepart'）;
stem（k,imag（x））%虚部图
19.%ProgramofStabilityTest（参考）
&本题详细信息可以参考中文版教材第274页例6.39
clear
num=input（'Typeinthenumeratorvector=（本题中应该是【1-0.20.5】）'）;
den=input（'Typeinthedenominatorvector=（本题中应该是【13.21.50.81.4】）'）;
N=max（length（num）,length（den））;
x=1;y0=0;S=0;zi=zeros（1,N-1）;
forn=1:
1000
[y,zf]=filter（num,den,x,zi）;
ifabs（y）<0.000001,break,end
x=0;
S=S+abs（y）;
y0=y;zi=zf;
end
ifn<1000
disp（'StableTransferFunction'）;
else
disp（'UnstableTransferFunc
升级会员 