dsp含40题程序.doc

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); gridxlabel(omega/pi); ylabel(Group delay, samples);title(Original Filter);F = 0:0.001:0.4;g = grpdelay(n,d,F,2); % Equalize the passbandGd = max(g)-g;% Design the allpass delay equalizernum,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);grpdelay(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 Butterworth 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 functionnum,den = butter(N,Wn,s);%计算滤波器的分子和分母% Compute and plot the frequency responseh,w=freqz(num,den);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=residue(num,den)4% Program 4% Cheby1 analog highpass filter design%Wp=0.75;Ws=0.35;Rs=43;Rp=0.5;N,Wn=cheb1ord(Wp,Ws,Rp,Rs,s);% Determine the coefficients of the transfer functionnum,den=cheby1(N,Rp,Wn,high,s);% Compute and plot the frequency responseomega=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(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)6. % Program to Design Type 1 Chebyshev analog Lowpass Filter% 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 = cheby1(N,Rp,Wc,s);% Compute and plot the frequency responseh,w=freqs(num,den);plot (w/pi,20*log10(abs(h);xlabel(Frequency, Hz); ylabel(magnitude, dB);title(IIR Chebyshev1 Lowpass Filter/analog)grid on7.% Determination of the Factored Form% of a Rational z-Transform%熟悉命令num,den=residuez(r,p,k)；r,p,k=residuez(num,den) clear allr = 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 bandpass Filterclear allWp = 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 onxlabel(omega/pi); ylabel(magnitude, dB);title(IIR Butterworth Bandpass Filter);9%draw the requested sequenceclearn=0:23;x=2*0.9.n;stem(n,x)title(xn)xlabel(n);ylabel(magnitude)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(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 function% Compute and plot the frequency responseh,w=freqs(num,den);plot (w/2*pi,20*log10(abs(h);xlabel(Frequency, Hz); ylabel(Gain, dB);title(iir ellip lowpass filter)gird on%num=1,-0.2,0.5; den=1,3.2,1.5,-0.8,1.4; y= impz(num,den,30); %inverse transformn=1:30;stem(n,y);grid on%clearN=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);gridxlabel(omega/pi); ylabel(Group delay, samples);title(The grpdelay of IIR Ellip Lowpass Filter)g = grpdelay(num,den,F,2); % Equalize the passbandGd = max(g)-g;% Design the allpass delay equalizernum,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);ylabel(Group delay, samples); %title(Group Delay Equalized Filter); 13%Program to draw signal and signal with noiseclearN=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 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 filterb=ones(4,1)/4; %window length is 4y=filter(b,1,x); %let signal go through the filterfigure(4)stem(n,y)title(signal through filter yn)xlabel(time index n); ylabel(amplitude)figure(5)plot(n,d,g:,n,s,b,n,x,c-,n,y,r-.);legend(noise,sn,xn,yn)%在图上加方框对其标注grid on14%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)xlabel(omega/); ylabel(phase(rad);title(phase response)%notice: here omega hasnt been归一化15% %已知系统的系统函数为：%用MATLAB的filter函数求hn的前20个样本clearN=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(hn)subplot(2,1,2)stem(0:N-1,y)title(yn)16%利用Hermann公式估计FIR低通滤波器的阶数。该滤波器的性能指标为：通带截止频率为1500Hz，阻带截至频率为1800Hz，通带纹波为=0.015，阻带纹波为=0.021，抽样频率为5000Hz%fedge=1500 1800;% unit Hzmval=1 0;%desired magnitude values in each banddev=0.015 0.021; %通带和阻带波纹Ft=5000;N=remezord(fedge,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+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)（another version of 13,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(dn,xn,sn)m=4;%input(please input the number of points you want filt=);m=(0:m-1)/m;y=filter(m,1,x);figureplot(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);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 = );mval = 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【n】为随机序列，幅度为0.6，分别绘制长度为40的受干扰序列以及中值滤波器的输出clearn=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,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的反变换IFFTfigure(1);stem(k,real(x); %实部图title(real part);figure(2);title(imagine part);stem(k,imag(x) %虚部图19. %Program of Stability Test（参考）&本题详细信息可以参考中文版教材第274页 例6.39clearnum = input(Type in the numerator vector =(本题中应该是【1 -0.2 0.5】) );den = input(Type in the denominator vector =(本题中应该是【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 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;endif n =1 disp(The system is not stable!)else disp(The system is stable!)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,mag,wt = remezord(fedge,mval,dev,Ft);%b = remez(N,fpts,mag,wt);fp=1500;fs=1800;delta_p=0.015;delta_s=0.021;Ft=5000N,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 onxlabel(omega/pi); ylabel(Gain, dB);21%clearx1=2.2,3,-1.5,4.2,-1.8; x2=0.8,-1,1.6,0.8; N = length(x1) ; M = length(x2) ; 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 alln=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,k/1024,X_dtft)title(DTFT/DFT)grid onfiguresubplot(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.2z-1+1.5z-2+0.8z-3+1.4z-4+3.6z-5+5.2z-6%用MATLAB将系统函数分解为二次多项式之积，并写出各二次多项式的表达式。%clear P=2.4,3.2,1.5,0.8,1.4,3.6,5.2; r=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 allb=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，设计满足该滤波器的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,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);%调用函数得出Kaiser的增益响应plot(W/pi,20*log10(abs(h);grid onxlabel(omega/pi); ylabel(Gain, dB);25%program_25 cleark=500; % number of frequency samples is 500num=1 -.2 .5 2 -.6; %Numerator coefficientsden=1 3.2 1.5 -.8 1.4; %Denominator coefficientsw=0:pi/(k-1):pi;h=freqz(num,den,w); %Compute the frequency responsesubplot(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 functiontitle(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.6 .8; l = length(x1)+length(x2)-1; DFT_X1=fft(x1,l);DFT_X2=fft(x2,l); % Compute the DFTsxn=ifft(DFT_X1.*DFT_X2) %Determine the IDFT of the product% Plot the sequence generated by DFT-based convolutionn=0:l-1;stem(n,xn)title(Result of DFT-based linear convolution)xlabel(Time index n);ylabel(Amplitude);27%已知IIR滤波器的系统函数为:%用matlab将系统函数表示为级联型结构形式,并写出各级联子系统的表达式.%b=2 5 1 -3 4 6; %Numeratora=1 3 -5 2 -4 3; %Denominatorsyms 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 allnum=2 5 1 -3 4 6;den=1 3 -5 2 -4 3;sos,G=tf2sos(num,den) %级联型结构形式,并写出各级联子系统的表达式28%用kaisers窗函数设计FIR数字高通滤波器，去滤波器的性能指标为：通带截止频率0.55，阻带截止频率0.45，通带和阻带波纹=0.04。绘制出该滤波器的增益响应。%program_28fpts=.45 .55; %the bandedges mag=0 1; %the desired magnitude valuesdev=.04 .04; %the ripples in each bandN,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,2825%num=1 -0.2 0.5 2 -0.6;den=1 3.2 1.5 -0.8 1.4; w=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);%其出单位冲击响应h,omega=freqz(b,1,256);%求出系统函数plot(omega/pi,-20*log(h)%画出增益函数29%绘制六点滑动平均滤波器的幅频特性和相频特性% w = 0:pi/255:pi; %何解呢num=ones(1,6)/6;h = freqz(num, 1, w);subplot(2,1,1)plot(w/pi,abs(h);gridtitle(Magnitude Spectrum)xlabel(omega/pi); ylabel(Magnitude)subplot(2,1,2)plot(w/pi,angle(h);gridtitle(Phase Spectrum)xlabel(omega/pi); ylabel(Phase, radians)%方法二：clear allb=ones(1,6)/6;freqz(b,1)30%原始序列为：sn=3*n*(0.8n),加性噪声dn为随机序列，幅度为0.6，受干扰的序列为xn=sn+dn,使用重叠相加法实现5点滑动平均滤波器对xn的处理。绘制未受干扰的序列sn和滤波器输出的有噪声序列clearR = 50;d = 1.2*(rand(1,R)-0.5);m = 0:1:R-1;s = 3*m.*(0.8.m);x =s+d;%重叠相加法b = ones(5,1)/5;%五阶的滑动平均滤波器y= filter(b,1,x);%x 是待输入的信号，要滤波的信号plot(m,s,r-,m,y,b-)legend(sn,yn);%或者是用 %y=fftfilt(b,x);%fftfilt的方法用来Y = FFTFILT(B,X) filters X with the FIR filter B using the overlap/add，method, using internal parameters (FFT size and block length) which guarantee efficient execution.%plot(n,s,b,n,y,r:)%legend(sn,yn)xlabel (Time index n);ylabel(Amplitude)31已知IIR滤波器的系统函数为%num = input(Numerator coefficient vector = );den = input(Denominator coefficient vector = );r1,p1,k1 = residuez(num,den);%系统并联I型分解r2,p2,k2 = residue(num,den);%系统并联II型分解disp(Parallel Form I)disp(Residues are);disp(r1);disp(Poles are at);disp(p1);disp(Constant value);disp(k1);disp(Parallel Form II)disp(Residues are);disp(r2);disp(Poles are at);disp(p2);disp(Constant value);disp(k2);32%用海明窗设计多频带FIR滤波器，该滤波器满足如下条件：在频率范围0到0.32内幅度为0.6，在频率范围0.35到0.65内幅度为0.2，在频率范围0.68到内幅度为0.8，绘制改滤波器的幅频特性。%fpts = 0 0.32 0.35 0.65 0.68 1;mval = 0.6 0.6 0.2 0.2 0.8 0.8;b = fir2(100,fpts,mval);h,w=freqz(b,1);plot(w/pi,abs(h);grid onxlabel(omega/pi); ylabel(Magnitude);33%已知滤波器的差分方程和输入信号为： ，绘制出该系统的输入序列和输出序列的包络。%包络,管庆的意思，就是用plot画出来，囧clear; n=0:63;x=0,0;%因为matlab矩阵的下标只能是正整数，为了下面的n-2要预填充两个0x1=(cos(0.1*n)+cos(0.4*n);x=x,x1;subplot(2,1,1);stem(n,x1);grid%画离散谱线的title(输入序列);y=-6.76195*x(n+3)+13.456335*x(n+2)-6.76169*x(n+1);%y是长度为64的序列subplot(2,1,2);%Index exceeds matrix dimensions.plot(y);gridtitle(输出包络);34%已知系统的系统函数为：%绘制该系统的零极点分布图num=1,-0.2,0.5,2,-0.6;%分子系数降幂排列den=1,3.2,1.5,-0.8,1.4;%分母系数降幂排列z,p,k=tf2zp(num,den);zplane(z,p);%画零极点的函数，表示毫无压力35%num=3,-4,2,-5,3,1;den=1,3,-5,2,-4,3;num=num/num(1);%先归一化，少用个乘法器k,alpha = tf2latc(num,den);%格型(latice)的参数计算%方法二 Program 8_4% Cascaded Lattice Realization of an% Allpass Transfer Function%format longden = input(Denominator coefficient vector = );k = poly2rc(den);knew = fliplr(k);disp(The lattice section multipliers