南理工數(shù)字信號處理matlab題庫

1、DSP題庫詳解 1已知3階橢圓IIR數(shù)字低通濾波器的性能指標為：通帶截止頻率0.4，通帶波紋為0.6dB，最小阻帶衰減為32dB。設計一個6階全通濾波器對其通帶的群延時進行均衡。繪制低通濾波器和級聯(lián)濾波器的群延時。%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

2、');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);              

3、60;    %設計六階的全通濾波器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)                   

4、;                         %DFILT:Digital Filter Implementation.%GdA,w = grpdelay(num,den,512);%figure(2);%plot(w/pi,GdH+GdA); grid%xlabel('omega/pi');ylabel('Group dela

5、y, samples');%title('Group Delay Equalized Filter'); 2設計巴特沃茲模擬低通濾波器，其濾波器的階數(shù)和3-dB截止頻率由鍵盤輸入，程序能根據(jù)輸入的參數(shù)，繪制濾波器的增益響應。% 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

6、 = input('3-dB cutoff angular frequency = ');% Determine the transfer functionnum,den = butter(N,Wn,'s');                              

7、;            %計算濾波器的分子和分母%  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

8、5) 3%          % To use the residuenum=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 coefficien

9、ts of the transfer functionnum,den=cheby1(N,Rp,Wn,'high','s');% Compute and plot the frequency responseomega=0:1:30*pi;                           

10、60;                   %w=0:0.01:pi;h=freqs(num,den,omega);                          %

11、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 p

12、art ');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 ri

13、pple 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=

14、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(

15、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數(shù)字帶通濾波器% Program to Design iir digital bandpass Filterclear allWp = 0.4 0.6;&

16、#160;                                                 &

17、#160;     %Wp and Ws 都是向量組，帶通的性質(zhì)Ws = 0.3 0.7;Rp = 0.6;Rs = 35;N,Wn = buttord(Wp, Ws, Rp, Rs);b,a = butter(N,Wn);                            

18、       %巴特臥茲濾波器的分子和分母h,w = freqz(b,a,256);                                 %求出在頻域的畫圖點，256，即要花的點數(shù)gain = 20*log1

19、0(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 on 10% Program t

20、o 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 = ');&

21、#160;                                                 &

22、#160;                      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 Mini

23、mum 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 lo

24、wpass filter)gird on  %num=1,-0.2,0.5;den=1,3.2,1.5,-0.8,1.4;y= impz(num,den,30);                                   

25、0;      %inverse transformn=1:30;stem(n,y);grid on  %clearN=5;Rp=0.8;Rs=35;Wn=0.35;                               &

26、#160;                                %去掉了pi的輸入num,den = ellip(N,Rp,Rs,Wn); F = 0:0.01:0.3;GdH,w = grpdelay(num,den,512);%gd=grpdelay(num,den

27、,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);                  &#

28、160;                                                % Equalize th

29、e 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)                  &#

30、160;                                                 &#

31、160;   %下面是另外一種方法%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

32、,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

33、 index n'); ylabel('amplitude')%generate moving average filterb=ones(4,1)/4;                                    &

34、#160;                 %window length is 4y=filter(b,1,x);                            

35、                            %let signal go through the filterfigure(4)stem(n,y)title('signal through filter yn')xlabel('time index n'); ylabel('a

36、mplitude')figure(5)plot(n,d,'g:',n,s,'b',n,x,'c-',n,y,'r-.');legend('noise','sn','xn','yn')                       

37、60;            %在圖上加方框?qū)ζ錁俗rid on14%N=10;x=ones(1,10);h,w=freqz(x,1,16);                             

38、0;                        subplot(2,1,1)plot(w,abs(h)xlabel('omega/'); ylabel('amplitude');title('amplitude response')subplot(2,1,2)plot(w,angle(h)xlabel(&#

39、39;omega/'); ylabel('phase(rad)');title('phase response')%notice:  here omega hasnt been歸一化 15%已知系統(tǒng)的系統(tǒng)函數(shù)為：%用MATLAB的filter函數(shù)求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;           

40、60;                                                 

41、60;         %x=1,zeros(1,19)h=filter(num,den,x);y=impz(num,den,N);                                &#

42、160;             %用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低通濾波器的階數(shù)。該濾波器的性能指標為：通帶截止頻率為1500Hz，阻帶截至頻率為1800Hz，通帶紋波為=0.015，阻帶紋波為=0.021，抽樣頻率為5000Hz%fed

43、ge=1500 1800;                                                

44、;         % unit Hzmval=1 0;                                      

45、0;                                 %desired magnitude values in each banddev=0.015 0.021;         &#

46、160;                                               %通帶和阻帶波紋Ft=5000;N=r

47、emezord(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%    

48、; 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)

49、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，受干擾的序列為：* MERGEFORMAT，分別繪制長度為40的原始序列，噪聲序列和受干擾序列，以及滑動平均濾波器的輸出。

50、%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','

51、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)

52、,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('請按降冪輸入分子系數(shù)=');b=input('請按降冪輸入分母系數(shù)=');y=filter(a,b,n) 16%fedge = input('Type in the bandedges = &

53、#39;);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的中

54、值濾波器程序，原始未受干擾的序列為：sn=3n(0.8)n,加性噪聲信號d【n】為隨機序列，幅度為0.6，分別繪制長度為40的受干擾序列以及中值濾波器的輸出clearn=1:40;s=3*n.*power(0.8,n);            %power表示指數(shù)函數(shù)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)     

55、;           %中值濾波器設計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);          &#

56、160;                                                 &#

57、160;              %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頁

58、 例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

59、);    if abs(y) < 0.000001, break, end    x = 0;    S = S + abs(y);    y0 = y;zi = zf;end if n < 1000    disp('Stable Transfer Function');else    disp('Unstable Transfer Function');en

60、d%法二定義法，從極點位置看穩(wěn)定性clear allnum=1 -0.2 0.5;den=1 3.2 1.5 0.8 1.4;z,p,k=tf2zp(num,den);zplane(z,p,k)p_m=abs(p);if max(p_m)>=1    disp('The system is not stable!')else    disp('The system is stable!')end 20 % Program 10_2% Design of Equiripple Linear

61、-Phase FIR Filters等波紋線性相位濾波器設計%詳細信息見書P449 例10.15%format long                                        

62、                                         %作用何在呢？%fedge = input('Band edges in Hz = ');%mval

63、= 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;

64、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); &

65、#160;                     %展示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 = le

66、ngth(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)                   &#

67、160;       %調(diào)用函數(shù)fft實現(xiàn)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)

68、title('DFT')subplot(2,1,2)plot(k/1024,X_dtft)title('DTFT')xlabel('omega/pi');ylabel('magnitude')     23%已知FIR濾波器的系統(tǒng)函數(shù)為：%H(Z)=2.4+3.2z-1 +1.5z-2 +0.8z-3 +1.4z-4 +3.6z-5 +5.2z-6%用MATLAB將系統(tǒng)函數(shù)分解為二次多項式之積，并寫出各二次多項式的表達式。%clear P=2.4,3.2,1.5,0.8,

69、1.4,3.6,5.2; r=roots(P);%調(diào)用函數(shù)計算 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=t

70、f2sos(b,1)                                        24%已知FIR數(shù)字低通濾波器的性能指標為通帶截止頻率0.35 ，阻帶截止頻率0.45 ，通帶和阻帶

71、波紋 =0.01，設計滿足該濾波器的Kaisers窗函數(shù)，繪制出Kaisers窗函數(shù)的增益響應。% Lowpass Filter Design Using the Kaiser Window%fpts = 0.3 0.4;                           %輸入wp和ws,不帶Pi的值哦mag = 1 0;  

72、;                               %輸入各頻段期望的幅度值dev = 0.01 0.01;                

73、0;      %輸入通帶和阻帶波紋N,Wn,beta,ftype = kaiserord(fpts, mag, dev);w = kaiser(N+1, beta);b=w/sum(w);                              

74、           %歸一化嗎？%b = fir1(N, Wn, w)h,W = freqz(b,1,512);              %調(diào)用函數(shù)得出Kaiser的增益響應plot(W/pi,20*log10(abs(h);grid onxlabel('omega/pi');ylabel('Gain, dB'); 25%

75、program_25    cleark=500;                      % number of frequency samples is 500num=1 -.2 .5 2 -.6;           %Numerator coeffi

76、cientsden=1 3.2 1.5 -.8 1.4;         %Denominator coefficientsw=0:pi/(k-1):pi;                 h=freqz(num,den,w);          &

77、#160;  %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 pha

78、se functiontitle('Phase Spectrum')xlabel('omega/pi');ylabel('Phase,radians')sos,g=tf2sos(num,den)                               26%

79、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); 

80、 % 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&

81、#39;);ylabel('Amplitude'); 27%已知IIR濾波器的系統(tǒng)函數(shù)為:                          %用matlab將系統(tǒng)函數(shù)表示為級聯(lián)型結(jié)構(gòu)形式,并寫出各級聯(lián)子系統(tǒng)的表達式.%b=2 5 1 -3 4 6;     &#

82、160;             %Numeratora=1 3 -5 2 -4 3;                  %Denominatorsyms z;            

83、60;                                                 

84、60;                                    %定義符號變量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;n

85、umerator=factor(num);                                              

86、60;                                                 

87、60;              denominator=factor(den);         %方法二：clear allnum=2 5 1 -3 4 6;den=1 3 -5 2 -4 3;sos,G=tf2sos(num,den)   %級聯(lián)型結(jié)構(gòu)形式,并寫出各級聯(lián)子系統(tǒng)的表達式      

88、;  28%用kaisers窗函數(shù)設計FIR數(shù)字高通濾波器，去濾波器的性能指標為：通帶截止頻率0.55，阻帶截止頻率0.45，通帶和阻帶波紋=0.04。繪制出該濾波器的增益響應。%program_28fpts=.45 .55;                           %the bandedgesmag=0 1;

89、60;                           %the desired magnitude valuesdev=.04 .04;                  &

90、#160;       %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

91、25,26,28 25 %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);                %求出頻率響應函數(shù)subplot(2,1,1);plot(w,abs(H);subplot(2,1,2);plot(w/pi,angle(H);   

92、0;               %畫出相頻曲線 %x1=2.2 3 -1.5 4.2 -1.8;x2=0.8 -1 1.6 0.8;x2=x2 0;                          

93、0;                         % 把數(shù)據(jù)對齊Fx1=fft(x1,5);Fx2=fft(x2,5);H=Fx1.*Fx2;                  

94、0;                  %總的系統(tǒng)函數(shù)x3=ifft(H,5);                              

95、          %傅里葉反變換n=0:4;stem(n,x3)      %n,wn,beta,type=kaiserord(0.45,0.55,0 1,0.04 0.04);                      b=f

96、ir1(n,wn,'high',kaiser(n+1,beta),'noscale');                           %其出單位沖擊響應h,omega=freqz(b,1,256);        

97、60;                                                 

98、60;   %求出系統(tǒng)函數(shù)plot(omega/pi,-20*log(h)                                          &

99、#160;                     %畫出增益函數(shù) 29%繪制六點滑動平均濾波器的幅頻特性和相頻特性%w = 0:pi/255:pi;                     

100、;                       %何解呢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')

101、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,使用重疊相加法實現(xiàn)5點滑動平均濾波器對xn的處理。繪制未受干擾的序列sn和濾波器輸出的有噪聲序列clearR = 50;d =

102、1.2*(rand(1,R)-0.5);m = 0:1:R-1;s = 3*m.*(0.8.m);x =s+d;                                        

103、0;  %重疊相加法b = ones(5,1)/5;                           %五階的滑動平均濾波器y= filter(b,1,x);