南理工数字信号处理matlab题库53页.doc

1、DSP题库详解1已知3阶椭圆IIR数字低通滤波器的性能指标为:通带截止频率04门,通带波纹为0.6dB, 最小阻带衰减为32dBo设计一个6阶全通滤波器对其通带的群延时进行均衡。绘制低通滤 波器和级联滤波器的群延时。%Progranm 1% Group-delay equalization of an IIR filter.%n,d = ellip(3,0.6,32,0.4);GdH,w = grpdelay(n,dz512);plot(w/pi,GdH); gridxlabel(omega/pi); ylabel(Group delay, samples);title(Original Fi

2、lter*);F = 0:0.001:0.4;g = grpdelay(nF,2);% Equalize the passbandGd = max(g)-g;% Design the al I pass delay equalizernumzden,tau = iirgrpdelay(6, F, 0 0.4, Gd);%设计六阶的全通滤波器bza=iirgrpdelay(6,F,0 0.4,Gd);Hel=dfilt.df2(bza);He=dfilt.df2(n,d);He_all=dfilt.cascade(He,Hel);grpdelay(He_all)%DFILT:Digital Fi

3、lter Implementation.%GdA,w = grpdelay(num,den,512);%figure(2);%plot(w/pi,GdH+GdA); grid%xlabel(omega/pi,);ylabel(Group delay, samples);%title(fGroup Delay Equalized Filter1);2设计巴特沃兹模拟低通滤波器，其滤波器的阶数和3dB截止频率由键盘输入，程序能根 据输入的参数，绘制滤波器的增益响应。% Program 2% Program to Design Butterworth Analog Lowpass Filter% T

4、ype in the filter order and passband edge frequency 设计边缘频率N = input(Type in filter order =);Wn = input(3-dB cutoff angular frequency =);% Determine the transfer functionnum,den = butter(N,Wn,s*);%计算滤波器的分子和分母% Compute and plot the frequency responseh,w=freqz( nu mzden);plot (w/pi,20*logl0(abs(h);xlab

5、el(、omega/pi); ylabel(Magnitude, dB);title(The magnitude respponse of the system1)axis(0 3 -80 5)3% To use the residuenum=lz-0.2,0.5,0,0；den=l,3.2z1.5,-0.8,1.4;rzP/k=residue (nu m,de n)4% Program 4% Chebyl analog highpass filter design%Wp=0.75;Ws=0.35;Rs=43;Rp=0.5;N,Wn=cheblord(Wp,WszRp,Rs/s,);% Det

6、ermine the coefficients of the transfer function numzden=chebyl(NzRpzWn/high7s);% Compute and plot the frequency responseomega=0:l:30*pi;%w=0:0.01:pi;h=freqs(num,denzomega);%plot(w/pi,20*logl0(abs(h)plot(omega/(2*pi),20*logl0(abs(h);xlabel(、omega/pr);ylabel(Mag nitude,dB);title(IIR Chebyshevl Highpa

7、ss Filter/digital)5.% Program to绘制30点序列的实部与虚部%n=0:29;x=0.2*exp(0.4+0.5*j).*n);subplot(lz2,l)stem(nzreal(x)title(The real part )；xlabel(rT);ylabel(mag nitude1)subplot( 1,2,2)stem(n,imag(x)title(The imagine part)xlabelCnJjylabelCmag nitude1)6.% Program to Design Type 1 Chebyshev analog Lowpass Filter%

8、 % Read in the filter order, passband edge frequency% and passband ripple in dBN = input(Order =);Wc = input(type in the 3 dB cutoff frequency=)；Rp = input( Type in the peak ripple value in passband:);% Determine the coefficients of the transfer functionnum,den = chebyl(N，Rp/Wc,s)；% Compute and plot

9、 the frequency responseh,w=freqs( nu m,de n);plot (w/pi,20*logl0(abs(h);xlabel(fFrequency，Hz1); ylabel(fmagnitude, dB);titleCIIR Chebyshevl Lowpass Filter/analog1)grid on7.% Determination of the Factored Form% of a Rational z-Transform%熟悉命令numzden=residuez(r,p,k)； r,p,k=residuez(numzden)clear allr =

10、 1 0.6 1.8;p =-3.2 2.4 2.4;k =0.2;num, den = residuez(r,pzk);dispCNumerator polynomial coefficients); disp(num)disp(Denominator polynomial coefficients); disp(den)8, 设计IIR数字带通滤波器% Program to Design iir digital bandpass Filter clear allWp = 0.4 0.6；的性质%Wp and Ws都是向量组，带通Ws = 0.3 0.7;Rp = 0.6;Rs = 35;N

11、Wn = buttord(Wp, Ws, Rp, Rs);bza = butter(N,Wn);h,w = freqz(b,az256);数%巴特卧兹滤波器的分子和分母%求出在频域的画图点，256,即要花的点gain = 20*logl0(abs(h);plot (w/pi,gain);grid onxlabel(、omega/pi)； ylabel(magnitude, dB);title(IIR Butterworth Bandpass Filter);9%draw the requested sequeneeclearn=0:23;x=2*0.9.An;stem(n,x)titleCx

12、tn1)xlabel(n);ylabel(mag nitude)grid on10%椭圆模拟低通滤波器设% Program to Design Elliptic analog Lowpass Filter 计% Read in the filter order, passband edge frequency,% passband ripple in dB and minimum stopband % attenuation in dBN = input(fOrder =);Wn = input(type in the 3 dB cutoff frequency= *);Rp = input(

13、ftype in the peak ripple value of passband in dB =);Rs = inputCtype in the Minimum stopband attenuation in dB = !);num,den = ellip(NzRp,RszWn,s)； %Determine the coefficients of the transfer function% Compute and plot the frequency responseh,w=freqs( nu m,de n);plot (w/2*pi,20*logl0(abs(h);xlabel(fFr

14、equencyz Hz1); ylabel(!Gain, dB!);titlefiir ellip lowpass filter?gird on%num=lz-0.2,0.5;den=lz3.2,1.5,-0.8,1.4;y= impz(num,de n,30);n=l:30;stem(n,y);grid on%clearN=5;Rp=0.8；Rs=35;Wn=0.35;%in verse 一tran sform%去掉了 pi的输入num,den = ellip(N,RpzRszW n);F = 0:0.01:0.3;GdH,w = grpdelay(numzden,512);%gd=grpd

15、elay(numzdenzw,2);群时延的另外一个公式 figure(l),plot(w/pi,GdH);gridxlabel(omega/pi); ylabel(Group delay, samples);title(The grpdelay of HR Ellip Lowpass Filter)g = grpdelay(numzdenF,2);% Equalizethe passbandGd = max(g)-g;% Desig n the al I pass delay equalizernum,den = iirgrpdelay(10, F, 0 0.3, Gd);He=dfilt.

16、df2(num,den);Hel=dfilt.df2(numzden);He_all=dfilt.cascade(He,Hel);grpdelay(He_all)%下面是另外一种方法%GdA,w = grpdelay(numzden,512);%plot(w/pi,GdH+GdA); grid on%xlabel(omega/pi,);ylabel(Group delay, samples);%title(Group Delay Equalized Filter);13%Program to draw signal and signal with noise clearN=40;n=0:N-l

17、s=3*(n.*(0.8.An);d=1.2.*rand(l,N)-0.6;x=s+d;figure(l)stem(nzs)xlabel(!time index nJ; ylabelCamplitude);titleCoriginal signal sn) figure(2)stem(nzd)title(noise dn)xlabeftime index rT); ylabel(amplitude)figure(3)stem(n,x)title(signal with noise xfn1)xlabel(time index rT); ylabeKmplitude1)%window leng

18、th is 4% let signal go through the filter%在图上加方框对其标注%generate moving average filterb=ones(4,l)/4;y=filter(b,l,x);figure(4)stem(nzy)title(signal through filter yn)xlabel(ftime index rT); ylabel(amplitude) figure(5)plot(n,d,g：ms,b 九 x,cYnMr-);Iegend(,noise7sn7xn7yn)grid on14%N=10;x=ones(l,10);h,w=freq

19、z(x,l,16);subplot(2,lzl)plot(w,abs(h)xlabel(omega/*); ylabelCamplitudeJjtitleCamplitude response)subplot(2,lz2)plot(wzangle(h)xlabel(omega/*); ylabel(,phase(rad);title(phase response*)%notice: here omega has nt bee n 归一化15%已知系统的系统函数为：%用MATLAB的filter函数求hn的前20个样本clearN=20;num=lz-0.2z0.5;den=lz3.2z1.5,

20、0.8,1.4;x=zeros(l,N);x(l)=l;%x=l,zeros(l,19)h=filter( nu mzde n,x);y=impz(numzden,N);%用 filter 和 z 反变换 impz 对比subplot(2zl,l)stem(0:N-l,h)titleChtn*)subplot(2zlz2)stem(0:N-lzy)titleCytn*)16% unit Hz%desired magnitude%通带和阻带波纹%利用Hermann公式估计FIR低通滤波器的阶数。该滤波器的性能指标为：通带截止频率 为1500Hz,阻带截至频率为1800Hz,通带纹波为=0.01

21、5,阻带纹波为=0.021,抽样频率 为 5000Hz fedge=1500 1800;mval=l 0;values in each banddev=0.015 0.021;Ft=5000;N=remezord(fedgezmval,dev,Ft)%第二种方法% clear all% format long % al=0.005309;a2=0.07114;a3=-0.4761;a4=0.00266;a5=0.5941;a6=0.4278;% bl=11.01217;b2=0.51244;% f=1500,1800;% w=2*pi*f/5000;% wp=w(l);ws=w(2);% de

22、lta_p=0.015;delta_s=0.021;% delta二0.015,0.021;% if delta_p=delta_s%m=l;n=2;%D=(al*(loglO(delta(m)A2+a2*loglO(delta(m)+a3)*loglO(delta(n)-(a4*(loglO(delt a(m)A2+a%5*loglO(delta(m)+a6);% F=bl+b2*(logl0(delta(m)-logl0(delta(n);% else%m=2; n=l;%D=(al*(loglO(delta(m)A2+a2*loglO(delta(m)+a3)*loglO(delta(n

23、)-(a4*(loglO(delt a(m)A2+a%5*loglO(delta(m)+a6);% F=bl+b2*(logl0(delta(m)-logl0(delta(n);% end% N=ceil(D-F*(ws-wp)/2/pi)A2)/(ws-wp)/2/pi)(another version of 13,14,15,16)13%编写四点滑动平均滤波器。原始未受干扰的序列为：，加性噪声信为随机序列，幅 度0.6,受干扰的序列为：* MERGEFORMAT,分别绘制长度为40的原始序列，噪声序列 和受干扰序列，以及滑动平均滤波器的输出。%R=40；d=1.2*rand(lzR)-0.

24、6;n=0:R-l;x=3*n.*(0.8).An;s=x+d;plot(rd,bla：n,x,r：rs);legend(,dn,；xn,z,snj,)m=4;%i nput(please in put the nu mber of points you want filt=);m=(0:m-l)/m;y=filter(m,l,x);figurePlot(n,sz nyg)legendCsn7yn)14%绘制长度为10点的矩形序列的16点离散傅里叶变换样本的幅度和相位% a=in put(fplease in put the mag nitude of the N-poi nt seque n

25、ce=J;N=in put(please in put the lenth of sequence=J;M=input(!please in put the lenth of the DFT sequence=); u=a*ones(lzN);h=fft(uzM);m=0:M-l;subplot(2,l,l);stem(mzabs(h),title(,magnitude,);subplot(2,lz2);stem(m,a ngle(h),title(phase);15%N=in put(please in put the lenth of the wanted sequence=); n=on

26、es(lzN);a=inputCW按降幕输入分子系数=);b=inputCW按降幕输入分母系数=);y=filter(a,b, n)16%fedge = input(Type in the bandedges =);mval = input(Desired magnitude values in each band = f);dev = input(Allowable deviation in each band =);FT = input(fType in the sampling frequency = f);N, fpts, mag, wt = remezord(fedge, mval,

27、 dev, FT); fprintf(Filter order is %d n,N)17%编写长度为5的中值滤波器程序，原始未受干扰的序列为：sn=3n(0.8)An,加性噪声信号d【n】为随机序列，幅度为0.6,分别绘制长度为40的受干扰序列 以及中值滤波器的输出clearn=1:40;s=3*n.*power(0.8,n);% power 表示指数函数figure(l)Plot(n,s);d=1.2* (rand(l,40)-0.5);x=s+d;figure(2)plot(nx)xx=medfiltl(xz5)%中值滤波器设计figure(3),plot(n,x,g几 xx,b)lege

28、ndCn7yn-)18% 已知 16 点序列刈n的 DFT 为：Xk=k/16 0=k=15 Xk=0 else 绘制 xn的实部 和虚部%k=0:15;X=k/15;x=ifft(X);%DFT 的反变换 IFFTfigure(l);stem(k,real(x); %实部图title(real part);figure(2);title(imagine part);stem(k,imag(x) % 虚部图19. % Prog ram of Stability Test (参考)&木题详细信息可以参考中文版教材第274页例6.39clearnum = i叩ut(*Type in the num

29、erator vector =(本题中应该是【1 -0.2 0.5);den = input(Type in the denominator vector =(本题中应该是1 3.2 1.5 0.8 1.4);N = max(le ngth( nu m),le ngth(de n);x = 1; y0 = 0; S = 0;zi = zeros(l,N-l);for n = 1:1000y,zf = filter( nu mzde npzi);if abs(y) =ldisp(The system is not stable!1) elsedisp(fThe system is stable!

30、1)end20 % Program 10_2% Design of Equiripple Linear-Phase FIR Filters 等波纹线性相位滤波器设计%详细信息见书P449例10.15%format long%作用何在呢?%fedge = input(Band edges in Hz =);%mval = input(Desired magnitude values in each band =);%dev = input(Desired ripple in each band =);%Ft= input(Sampling frequency in Hz =);%N,fpts,m

31、agzwt = remezord(fedge,mval,dev,Ft);%b = remez(Nzfptszmag,wt);fp=1500;fs=1800;delta_p=0.015;delta_s=0.021;Ft=5000N,fzmagz Wt=remezord(fpzfszlzOzdelta_p,delta_s/Ft)b=remez(N,f,magzWt)%fir filter coefficientsh,w = freqz(bzl,256);%展示256个点的频域图形plot(w/piz20*logl0(abs(h);grid onxlabel(omega/pi); ylabel(Ga

32、inz dB);21%clearxl=2.2,3,-1.5,4.2z-1.8;x2=0.8z-l,1.6,0.8;N = length(xl);M = Iength(x2);L=N + Mre=conv(xl,x2);% 卷积n=l:l:L;stem(n/e 刖)22%clearsyms wn=0:l:15;x=cos(pi*n./2);X=fft(x,16);stem(n, abs(X)%调用函数fft实现x的DFT变换figure(l)w=(0:0.1:pi)/piY=freqz(x,l,w)stem(w,abs(Y)%法二clear alln=0:15;k=0:1023;x=cos(n*

33、pi/2);X_dft=fft(x,16);X_dtft=fft(x,1024);plot(n/16,X_dft/rozk/1024zX_dtft)titleCDTFT/DFT*)grid onfiguresubplot(2zl,l)stem(nzX_dft)titleCDFT*)subplot(2,l,2)plot(k/1024zX_dtft)titleCDTFT*)xlabel(、omega八pr);ylabel(mag nitude)23%已知FIR滤波器的系统函数为：%H(Z)=2.4+3.2z-l +1.5Z-2 +0.8z-3 +1.4z4 +3.6z-5 +52z6%用MATLA

34、B将系统函数分解为二次多项式之积，并写出各二次多项式的表达式。%clearP=2.4z3.2z1.5/0.8/1.4z3.6,5.2;r=roots(P);%调用函数计算syms zsl=simple(z-r(l)*(z-r(2);dl=simple(sl./zA2)s2=simple(z-r(3)*(z-r(4);d2=simple(s2./zA2)s3=simple(z-r(5)*(z-r(6);d3=simple(s3./zA2)Q=2.4*dl*d2*d3%方法二：clear allb=2.4 3.2 1.5 0.8 1.4 3.6 5.2;sos,G=tf2sos(b,l)24%己知

35、FIR数字低通滤波器的性能指标为通带截止频率0.35 ,阻带截止频率045 ,通带 和阻带波纹=0.01,设计满足该滤波器的Kaisers窗函数，绘制出Kaisers窗函数的增益响 应。% Lowpass Filter Design Using the Kaiser Window%fpts = 0.3 0.4;%输入wp和ws,不带Pi的值哦mag = 1 0;%输入各频段期望的幅度值dev = 0.01 0.01;%输入通带和阻带波纹N/Wnzbetazftype = kaiserord(fpts, mag, dev); w = kaiser(N+l, beta);b=w/sum(w);%归

36、一化吗？%b = firl(N, Wn, w)hzW = freqz(b,l/512);%调用函数得出Kaiser的增益响应plot(W/pi,20*logl0(abs(h);grid onxlabel(、omega/pi)；ylabelCGain, dB);25%program_25clearnum=l -.2 .5 2 -.6;%Numerator coefficientsden=l 3-2 1.5 -.8 1-4;%Denominator coefficientsw=0:pi/(k-l):pi;h=freqz( nu m,denzw);%Compute the frequency res

37、ponsesubplot(l,2,l)plot(w/pi,abs(h)title(Magnitude Spectrum1)xlabel(、omega/pi);ylabel(Mag nitude)subplot( 1,2,2)plot(w/pi,unwrap(angle(h)%unwrapped phase functiontitle(Phase Spectrum)xlabel(、omega八pi);ylabelCPhasezradia ns)sos,g=tf2sos( nu mzde n)26%xl=2.2 3-1.5 4.2 -1.8;x2=.8 -1 1.6 .8;I = len gth(

38、xl)+length(x2)-l;DFT_Xl=fft(xlzl);DFT_X2=fft(x2,l); % Compute the DFTsxn=ifft(DFT_Xl.*DFT_X2)%Determine the IDFT of the product% Plot the seque nee gen erated by DFT-based con volutio nn=0:l-l;stem(n,x n)title(Result of DFT-based linear convolutiorT)xlabel(fTime index nJjylabelCAmplitude);27%已知IIR滤波

39、器的系统函数为：%用matlab将系统函数表示为级联型结构形式，并写出各级联子系统的表达式.%b=2 5 1 3 4 6;%Numeratora=l 3 5 2 -4 3;%Denomin atorsyms 石变量%定义符号num=2+5*zA-l+zA(-2)-3*zA(-3)+4*zA(-4)+6*zA(-5);den=l+3*zA-l-5*zA-2+2*zA-3-4*zA-4+3*zA-5;nu merator=factor(num);denomin ator=factor(de n);%方法二：clear allnum=2 5 1 -34 6;den=l 3-5 2-4 3;sos,G

40、tf2sos(numzden) %级联型结构形式，并写出各级联子系统的表达式28%用kaisers窗函数设计FIR数字高通滤波器，去滤波器的性能指标为：通带截止频率 0.55n,阻带截止频率0.45n,通带和阻带波纹6=0.04o绘制出该滤波器的增益响应。%program_28fpts=.45 .55;%the bandedgesmag=0 1;%the desired magnitude valuesdev=.O4 .04;%the ripples in each bandN/Wn zbetazftype = kaiserord(fpts,mag,dev) kw = kaiser(N+lz

41、beta);b = firl(N/WrVhigh： kw);h,omega = freqz(b,l);plot(omega/piz20*logl0(abs(h);grid;xlabel(omega/pi); ylabel(Gainz dB);another version of 25,26,2825%num=l -0.2 0.5 2 -0.6;den=l 3.2 1.5 -0.8 1.4; w=0:pi/255:2*pi;H=freqz(num,de nw);subplot(2,lzl)； plot(w,abs(H);subplot(2,l,2); plot(w/pi,angle(H);%xl

42、2.2 3 15 4.2 -1.8；x2=0.81.6 0.8;x2=x2 0;Fxl=fft(xl,5);Fx2=fft(x2,5);H=Fxl.*Fx2;x3=ifft(H,5);n=0:4;%求岀频率响应函数%画出相频曲线%把数据对齐%总的系统函数%傅里叶反变换stem(n,x3)%其出单位冲击响应%求出系统函数%画出增益函数n,wn,betaztype=kaiserord(0.45z0.55,0 l,0.04 0.04); b=firl (rwn/highkaiserf n+lbeta), no scale);h,omega=freqz(bzl/256);plot(omega/pi,-20*log(h)29%绘制六点滑动平均滤波器的幅频特性和相频特性%w = 0:pi/255:pi;% 何解呢num=on es(l,6)/6;h = freqz(num, 1, w);subplot(2,l,l)plot(w/pi,abs(h);gridtitle(Magnitude Spectrum1)xlabel(、omega/pi); ylabel(Magnitude)subplot(2,lz2)plot(w/pi,a ngle(h);gridtitle(Phase Spectrum)xlabel(、omega