《dsp各种程序.docx》由会员分享，可在线阅读，更多相关《dsp各种程序.docx（26页珍藏版）》请在三一办公上搜索。
1、dsp各种程序DSP题库详解 1已知3阶椭圆IIR数字低通滤波器的性能指标为：通带截止频率0.4，通带波纹为0.6dB，最小阻带衰减为32dB。设计一个6阶全通滤波器对其通带的群延时进行均衡。绘制低通滤波器和级联滤波器的群延时。 %Progranm 1 % Group-delay equalization of an IIR filter. % n,d = ellip(3,0.6,32,0.4); GdH,w = grpdelay(n,d,512); plot(w/pi,GdH); grid xlabel(omega/pi); ylabel(Group delay, samples); tit
2、le(Original Filter); F = 0:0.001:0.4; g = grpdelay(n,d,F,2); % Equalize the passband Gd = max(g)-g; % Design the allpass delay equalizer num,den,tau = iirgrpdelay(6, F, 0 0.4, Gd); %设计六阶的全通滤波器 b,a=iirgrpdelay(6,F,0 0.4,Gd); He1=dfilt.df2(b,a); He=dfilt.df2(n,d); He_all=dfilt.cascade(He,He1); grpdela
3、y(He_all) %DFILT:Digital Filter Implementation. GdA,w = grpdelay(num,den,512); figure(2); plot(w/pi,GdH+GdA); grid xlabel(omega/pi);ylabel(Group delay, samples); title(Group Delay Equalized Filter); 2设计巴特沃兹模拟低通滤波器，其滤波器的阶数和3-dB截止频率由键盘输入，程序能根据输入的参数，绘制滤波器的增益响应。 % Program 2 % Program to Design Butterwor
4、th Analog Lowpass Filter % % Type in the filter order and passband edge frequency设计边缘频率 N = input(Type in filter order = ); Wn = input(3-dB cutoff angular frequency = ); % Determine the transfer function num,den = butter(N,Wn,s); %计算滤波器的分子和分母 % Compute and plot the frequency response h,w=freqz(num,d
5、en); plot (w/pi,20*log10(abs(h); xlabel(omega/pi); ylabel(Magnitude, dB); title(The magnitude respponse of the system) axis(0 3 -80 5) 3% % To use the residue % num=1,-0.2,0.5,0,0; den=1,3.2,1.5,-0.8,1.4; r,p,k=residuez(num,den) 4% Program 4 % Cheby1 analog highpass filter design % Wp=0.75;Ws=0.35;R
6、s=43;Rp=0.5; N,Wn=cheb1ord(Wp,Ws,Rp,Rs,s); % Determine the coefficients of the transfer function num,den=cheby1(N,Rp,Wn,high,s); % Compute and plot the frequency response 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
7、/pi);ylabel(Magnitude,dB); title(IIR Chebyshev1 Highpass Filter/digital) 5. % Program to 绘制30点序列的实部与虚部 % n=0:29; x=0.2*exp(0.4+0.5*j).*n); subplot(1,2,1) stem(n,real(x) title(The real part ); xlabel(n);ylabel(magnitude) subplot(1,2,2) stem(n,imag(x) title(The imagine part) xlabel(n);ylabel(magnitude
8、) 6. % Program to Design Type 1 Chebyshev analog Lowpass Filter % % Read in the filter order, passband edge frequency % and passband ripple in dB N = input(Order = ); Wc = input(type in the 3 dB cutoff frequency=); Rp = input( Type in the peak ripple value in passband: ); % Determine the coefficient
9、s of the transfer function num,den = cheby1(N,Rp,Wc,s); % Compute and plot the frequency response h,w=freqs(num,den); plot (w/pi,20*log10(abs(h); xlabel(Frequency, Hz); ylabel(magnitude, dB); title(IIR Chebyshev1 Lowpass Filter/analog) grid on 7. % Determination of the Factored Form % of a Rational
10、z-Transform %熟悉命令num,den=residuez(r,p,k)；r,p,k=residuez(num,den) clear all r = 1 0.6 1.8; p =-3.2 2.4 2.4; k =0.2; num, den = residuez(r,p,k); disp(Numerator polynomial coefficients); disp(num) disp(Denominator polynomial coefficients); disp(den) 8, 设计IIR数字带通滤波器 % Program to Design iir digital bandp
11、ass Filter clear all Wp = 0.4 0.6; %Wp and Ws 都是向量组，带通的性质 Ws = 0.3 0.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);grid on xlabel(omega/pi); ylabel(magnitude, dB); title(II
12、R Butterworth Bandpass Filter); 9%draw the requested sequence clear n=0:23; x=2*0.9.n; stem(n,x) title(xn) xlabel(n);ylabel(magnitude) grid on 10% % Program to Design Elliptic analog Lowpass Filter %椭圆模拟低通滤波器设计 % Read in the filter order, passband edge frequency, % passband ripple in dB and minimum
13、stopband % attenuation in dB N = input(Order = ); Wn = input(type in the 3 dB cutoff frequency= ); Rp = input(type in the peak ripple value of passband in dB = ); Rs = input(type in the Minimum stopband attenuation in dB = ); num,den = ellip(N,Rp,Rs,Wn,s); %Determine the coefficients of the transfer
14、 function % Compute and plot the frequency response h,w=freqs(num,den); plot (w/2*pi,20*log10(abs(h); xlabel(Frequency, Hz); ylabel(Gain, dB); title(iir ellip lowpass filter) grid on % num=1,-0.2,0.5; den=1,3.2,1.5,-0.8,1.4; y= impz(num,den,30); n=1:30; stem(n,y);grid on %inverse transform % clear N
15、=5; Rp=0.8; Rs=35; Wn=0.35; num,den = ellip(N,Rp,Rs,Wn); F = 0:0.01:0.3; %去掉了pi的输入 GdH,w = grpdelay(num,den,512);%gd=grpdelay(num,den,w,2);群时延的另外一个公式 figure(1),plot(w/pi,GdH);grid xlabel(omega/pi); ylabel(Group delay, samples); title(The grpdelay of IIR Ellip Lowpass Filter) g = grpdelay(num,den,F,2
16、); % Equalize the passband Gd = max(g)-g; % Design the allpass delay equalizer num,den = iirgrpdelay(10, F, 0 0.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); grid on xlabel(omega/pi)
17、;ylabel(Group delay, samples); title(Group Delay Equalized Filter); 13%Program to draw signal and signal with noise 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(time index n); ylabel(amplitude);title(original signal sn) figure(2) stem(n,d) title(noise
18、 dn) xlabel(time index n); ylabel(amplitude) figure(3) stem(n,x) title(signal with noise xn) xlabel(time index n); ylabel(amplitude) % % %generate moving average filter b=ones(4,1)/4; %window length is 4 y=filter(b,1,x); figure(4) stem(n,y) title(signal through filter yn) xlabel(time index n); ylabe
19、l(amplitude) figure(5) plot(n,d,g:,n,s,b,n,x,c-,n,y,r-.); legend(noise,sn,xn,yn) grid on 14% %let signal go through the filter %在图上加方框对其标注 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(amplitude response) subplot(2,1,2) plot(w,angle(h) xl
20、abel(omega/); ylabel(phase(rad);title(phase response) %notice: here omega hasnt been归一化 15% %已知系统的系统函数为： %用MATLAB的filter函数求hn的前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) s
21、tem(0:N-1,h) title(hn) subplot(2,1,2) stem(0:N-1,y) title(yn) 16% %利用Hermann公式估计FIR低通滤波器的阶数。该滤波器的性能指标为：通带截止频率为1500Hz，阻带截至频率为1800Hz，通带纹波为dp=0.015，阻带纹波为ds=0.021，抽样频率为5000Hz % fedge=1500 1800; mval=1 0; % unit Hz %desired magnitude values in each band dev=0.015 0.021; %通带和阻带波纹 Ft=5000; N=remezord(fedge
22、,mval,dev,Ft) %第二种方法 % clear all % format long % 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; % if delta_p=delta_s % m=1;n=2; %D=(a1*(log10(delta(m)2+
23、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=
24、ceil(D-F*(ws-wp)/2/pi)2)/(ws-wp)/2/pi) nsn=3n(0.8)，13%编写四点滑动平均滤波器。原始未受干扰的序列为：加性噪声信dn为随机序列，幅度0.6，受干扰的序列为：xn=sn+dn，分别绘制长度为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(dn,xn,sn) m=4;%input(please input the number of points
25、 you want filt=); m=(0:m-1)/m; y=filter(m,1,x); figure plot(n,s,n,y,g) legend(sn,yn) 14%绘制长度为10点的矩形序列的16点离散傅里叶变换样本的幅度和相位 % a=input(please input the magnitude of the N-point sequence=); N=input(please input the lenth of sequence=); M=input(please input the lenth of the DFT sequence=); u=a*ones(1,N);
26、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(please input the lenth of the wanted sequence= ); n=ones(1,N); a=input(请按降幂输入分子系数=); b=input(请按降幂输入分母系数=); y=filter(a,b,n) 16% % fedge = input(Type in the bandedges = ); m
27、val = input(Desired magnitude values in each band = ); dev = input(Allowable deviation in each band = ); FT = input(Type in the sampling frequency = ); N, fpts, mag, wt = remezord(fedge, mval, dev, FT); fprintf(Filter order is %d n,N) 17%编写长度为5的中值滤波器程序，原始未受干扰的序列为：sn=3n(0.8)n, 加性噪声信号d为随机序列，幅度为0.6，分别绘
28、制长度为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.6); x=s+d; figure(2) plot(n,x) xx=medfilt1(x,5) %中值滤波器设计 figure(3),plot(n,x,g,n,xx,b) legend(n,yn) 18%已知16点序列xn的DFT为： Xk=k/16 0=k=15 Xk=0 else 绘制xn的实部和虚部 % k=0:15; X=k/15; x=ifft(X); %DFT
29、的反变换IFFT figure(1); stem(k,real(x); %实部图 title(real part); figure(2); title(imagine part); stem(k,imag(x) %虚部图 19. %Program of Stability Test &本题详细信息可以参考中文版教材第274页 例6.39 clear num = (1 -0.2 0.5); den = (1 3.2 1.5 0.8 1.4); N = max(length(num),length(den); x = 1; y0 = 0; S = 0;zi = zeros(1,N-1); for
30、n = 1:1000 y,zf = filter(num,den,x,zi); if abs(y) 0.000001, break, end x = 0; S = S + abs(y); y0 = y;zi = zf; end if n =1 disp(The system is not stable!) else disp(The system is stable!) end 20 % Program 10_2 % Design of Equiripple Linear-Phase FIR Filters等波纹线性相位滤波器设计 %详细信息见书P449 例10.15 % format lon
31、g %作用何在呢？ %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,mag,wt = remezord(fedge,mval,dev,Ft); %b = remez(N,fpts,mag,wt); fp=1500;fs=1800;delta_p=0.015;delta_s
32、=0.021;Ft=5000 N,f,mag,Wt=remezord(fp,fs,1,0,delta_p,delta_s,Ft) b=remez(N,f,mag,Wt) %fir filter coefficients h,w = freqz(b,1,256); %展示256个点的频域图形 plot(w/pi,20*log10(abs(h);grid on xlabel(omega/pi); ylabel(Gain, dB); 21% clear x1=2.2,3,-1.5,4.2,-1.8; x2=0.8,-1,1.6,0.8; N = length(x1) ; M = length(x2)
33、 ; L = N + M - 1 ; re=conv(x1,x2);%卷积 n=1:1:L; stem(n,re,r) 22% %clearsyms wn=0:1:15;x=cos(pi*n./2); X=fft(x,16);stem(n, abs(X) %调用函数fft实现x的DFT变换figure(1)w=(0:0.1:pi)/piY=freqz(x,1,w)stem(w,abs(Y) %法二 clear all n=0:15; k=0:1023; x=cos(n*pi/2); X_dft=fft(x,16); X_dtft=fft(x,1024); plot(n/16,X_dft,ro,
34、k/1024,X_dtft) title(DTFT/DFT) grid on figure subplot(2,1,1) stem(n,X_dft) title(DFT) subplot(2,1,2) plot(k/1024,X_dtft) title(DTFT) xlabel(omega/pi);ylabel(magnitude) 23%已知FIR滤波器的系统函数为： %H(Z)=2.4+3.2+1.5+0.8+1.4+3.6+5.2 %用MATLAB将系统函数分解为二次多项式之积，并写出各二次多项式的表达式。 % clear P=2.4,3.2,1.5,0.8,1.4,3.6,5.2; r
35、=roots(P);%调用函数计算 syms z s1=simple(z-r(1)*(z-r(2); d1=simple(s1./z2) s2=simple(z-r(3)*(z-r(4); d2=simple(s2./z2) s3=simple(z-r(5)*(z-r(6); d3=simple(s3./z2) Q=2.4*d1*d2*d3 %方法二： clear all b=2.4 3.2 1.5 0.8 1.4 3.6 5.2; sos,G=tf2sos(b,1) 24% %已知FIR数字低通滤波器的性能指标为通带截止频率0.35，阻带截止频率0.45，通带和阻带波纹=0.01，设计满足该
36、滤波器的Kaisers窗函数，绘制出Kaisers窗函数的增益响应。 % Lowpass Filter Design Using the Kaiser Window % fpts = 0.35 0.45; %输入wp和ws,不带Pi的值哦 mag = 1 0; %输入各频段期望的幅度值 dev = 0.01 0.01; %输入通带和阻带波纹 N,Wn,beta,ftype = kaiserord(fpts, mag, dev); w = kaiser(N+1, beta); b=w/sum(w); %归一化吗？ %b = fir1(N, Wn, w) h,W = freqz(b,1,512);
37、 %调用函数得出Kaiser的增益响应 plot(W/pi,20*log10(abs(h); grid on xlabel(omega/pi); ylabel(Gain, dB); 25%program_25 clear k=500; % number of frequency samples is 500 num=1 -.2 .5 2 -.6; %Numerator coefficients den=1 3.2 1.5 -.8 1.4; %Denominator coefficients w=0:pi/(k-1):pi; h=freqz(num,den,w); %Compute the fr
38、equency response subplot(1,2,1) plot(w/pi,abs(h) title(Magnitude Spectrum) xlabel(omega/pi);ylabel(Magnitude) subplot(1,2,2) plot(w/pi,unwrap(angle(h) %unwrapped phase function title(Phase Spectrum) xlabel(omega/pi);ylabel(Phase,radians) sos,g=tf2sos(num,den) 26% % x1=2.2 3 -1.5 4.2 -1.8; x2=.8 -1 1
39、.6 .8; l = length(x1)+length(x2)-1; DFT_X1=fft(x1,l); DFT_X2=fft(x2,l); % Compute the DFTs xn=ifft(DFT_X1.*DFT_X2) %Determine the IDFT of the product % Plot the sequence generated by DFT-based convolution n=0:l-1; stem(n,xn) title(Result of DFT-based linear convolution) xlabel(Time index n);ylabel(A
40、mplitude); 27% %已知IIR滤波器的系统函数为: %用matlab将系统函数表示为级联型结构形式,并写出各级联子系统的表达式. % b=2 5 1 -3 4 6; %Numerator a=1 3 -5 2 -4 3; %Denominator syms z; %定义符号变量 num=2+5*z-1+z(-2)-3*z(-3)+4*z(-4)+6*z(-5); den=1+3*z-1-5*z-2+2*z-3-4*z-4+3*z-5; numerator=factor(num); denominator=factor(den); %方法二： clear all num=2 5 1
41、-3 4 6; den=1 3 -5 2 -4 3; 2+5z-1+z-2-3z-3+4z-4+6z-5H(z)=1+3z-1-5z-2+2z-3-4z-4+3z-5sos,G=tf2sos(num,den) %级联型结构形式,并写出各级联子系统的表达式 28%用kaisers窗函数设计FIR数字高通滤波器，去滤波器的性能指标为：通带截止频率0.55，阻带截止频率0.45，通带和阻带波纹=0.04。绘制出该滤波器的增益响应。 %program_28 fpts=.45 .55; %the bandedges mag=0 1; %the desired magnitude values dev=.
42、04 .04; %the ripples in each band N,Wn,beta,ftype = kaiserord(fpts,mag,dev) kw = kaiser(N+1,beta); b = fir1(N,Wn,high, kw); h,omega = freqz(b,1); plot(omega/pi,20*log10(abs(h);grid; xlabel(omega/pi); ylabel(Gain, dB); another version of 25,26,28 25 % num=1 -0.2 0.5 2 -0.6; den=1 3.2 1.5 -0.8 1.4; w=
43、0:pi/255:2*pi; H=freqz(num,den,w); subplot(2,1,1); plot(w,abs(H); subplot(2,1,2); plot(w/pi,angle(H); %求出频率响应函数 %画出相频曲线 % x1=2.2 3 -1.5 4.2 -1.8; x2=0.8 -1 1.6 0.8; x2=x2 0; Fx1=fft(x1,5); Fx2=fft(x2,5); % 把数据对齐 H=Fx1.*Fx2; x3=ifft(H,5); n=0:4; stem(n,x3) %总的系统函数 %傅里叶反变换 % n,wn,beta,type=kaiserord(0.45,0.55,0 1,0.04 0.04); b=fir1(n,wn,high,kaiser(n+1,beta),noscale);
