数字信号处理机考dsp上机参考

CodingandAnswerbyhplonline,notingbyhplonline,typesettingbyvaa4第14页共14页NorightsreservedandWearenotresponsibleforanythingTOC\o"1-1"\h\z\u1. DFT定义，矩阵形式运算，循环运算 12. fir1,kaiserord,Kaiser,freqz 23. cheb1ord,cheb1 24. 系统函数的二阶分解，tf2sos,tf2zp,sort,unique 35. 零极点形式到直接型的转换 36. 阶数估计（表10.2）,hann,fir1，根据定义加窗 47. freqz,real,imag,abs,angle,unwrap,phase 48. fir1,hamming窗 59. 直接型转零极点型，tf2zp 610. butter,buttord 611. 部分分式展开（逆），从传输函数计算冲击响应，residuez,impz 712. freqz,real,imag,abs,phase,angle,unwrap 713. ellipord,ellip,freqz 814. fir1,hann,分贝转换 915. 信号能量计算 916. 各IIR阶数估计，buttord,cheb1ord,cheb2ord,ellipord 917. square,freqz 1018. 线性卷积，此处用fft实现 1119. 圆周卷积，要求用基于DFT的方法实现 1120. 部分分式展开，逆z变换，residuez 1221. freqz 1322. 信号通过系统的输出，filter 1323. 信号通过系统的输出，filter 1324. 重要提示 14DFT定义，矩阵形式运算，循环运算%WritingaMATLABprogramtocompute128-pointDFTof%thefollowingsequence,youmustfirstlyuse%DFTdefinition(directlycomputingDFT)tocompute%anduseMATLABfunctiontotesttheresult.%Plotthetworesultsinonefigure.%g[n]=sin(pi*n/4),0<=n<=31N=32;M=128;n=0:(N-1);m=0:(M-1);g=sin(pi*n/4);%DFTdefinition

%matrixexpressionofDFTG1=g*exp(-2*j*pi/M).^(n'*m);

%followingisthedirectconversion%fromdefinitiontocode%forii=0:M-1%G1(ii+1)=0;%forjj=0:N-1%G1(ii+1)=G1(ii+1)+g(jj+1)*exp(-2*j*pi/M*ii*jj);%end%end

subplot(2,1,1);plot(m,abs(G1));

%fftrealizedbyMATLABG2=fft(g,M);subplot(2,1,2);plot(m,abs(G2));fir1,kaiserord,Kaiser,freqz%Usingthefunctionfir1andwindowofKaiser,%designalinear-phaseFIRlowpassfilter%meetingthefollowingspecifications:%passbandedgefrequency=2kHz,%stopbandedgefrequency=2.5kHz,%passbandrippleδp=0.005,%stopbandrippleδs=0.005,%andsamplingrateof10kHz.%Plotitsgainandphaseresponses%andcheckifitmeetsthespecifications?fsamp=10000;fcuts=[20002500];mags=[10];devs=[0.0050.005];[n,Wn,beta,ftype]=kaiserord(fcuts,mags,devs,fsamp);hh=fir1(n,Wn,kaiser(n+1,beta),'noscale');freqz(hh)cheb1ord,cheb1%DesignaType1ChebyshevIIRlowpassfiltermeeting%thespecificationsasbelow:%samplingrateof12kHz,%passbandedgefrequencyof2.1kHz,%stopbandedgefrequencyof2.7kHz,%passbandrippleof0.6dB,%andaminimumstopbandattenuationof45dB.%Writedowntheexactexpressionfor%thetransferfunctiongenerated.%Doesyourdesignmeetthespecifications?Fs=12000;Fn=Fs/2;%NyquistcriterionWp=2100/Fn;Ws=2700/Fn;Rp=0.6;Rs=45;[n,Wn]=cheb1ord(Wp,Ws,Rp,Rs)[b,a]=cheby1(n,Rp,Wn)freqz(b,a,512,Fs);

%b=%0.00010.00090.00320.0064...%0.00800.00640.00320.00090.0001%a=%1.0000-4.681811.1095-16.8220...%17.5480-12.83006.4101-2.00620.3038系统函数的二阶分解，tf2sos,tf2zp,sort,unique%UsingMATLABdeterminethesecond-orderfactoredform%ofthefollowingz-transforms.Andshowtheirpole-zeroplots.%ThendetermineallpossibleROCsofabovez-transform.%4*z^4+15.6*z^3+6*z^2+2.4*z-6.4%G1(z)=%3*z^4+2.4*z^3+6.3*z^2-11.4*z+6num=[4,15.6,6,2.4,-6.4];den=[3,2.4,6.3,-11.4,6];[sosg]=tf2sos(num,den)[zpk]=tf2zp(num,den);pp=sort(unique([0,abs(p)',inf]));forii=1:length(pp)-1fprintf('possibleROC%d:(%f,%f)\n',ii,pp(ii),pp(ii+1));end零极点形式到直接型的转换%WriteaMATLABprogramtocomputeanddisplay%therationalz-transformfromitszero,%polesandgainconstant.IftheDTFTof%thecausalsystemexists,trytoplot%it'smagnitudeandphasespectra.%Usingthisprogramtodeterminetherationalform%ofaz-transformwhosezeroareat%ξ1=1.2,ξ2=2.3-j0.5,ξ3=-0.4+j0.2,ξ4=-0.4-j0.2,ξ5=2.3+j0.5;%theploesareat%λ1=0.5,λ2=-0.75+j0.2,λ3=0.6+j0.7,λ4=0.6-j0.7,λ5=-0.75-j0.2;%andthegainconstantkis2.1.

functionmain

z=[1.2,2.3-j*0.5,-0.4+j*0.2,-0.4-j*0.2,2.3+j*0.5];p=[0.5,-0.75,+j*0.2,0.6+j*0.7,0.6-j*0.7,-0.75-j*0.2];k=2.1;rzt(z,p,k);return;

functionrzt(z,p,k)%computeanddisplaynum=k*poly(z)den=poly(p)freqz(num,den);return;阶数估计（表10.2）,hann,fir1，根据定义加窗%DesignanFIRhanningdigitallowpassfilter%withthefollowingspecifications:%passbandedgefrequency0.3*pi,%stopbandedgefrequency0.4*pi,%stopbandattenuation40dB.%Pleaseplotfrequency-magnitituderesponse%andcheckifyourdesignfitthespecification。wp=0.3*pi;ws=0.4*pi;wn=(wp+ws)/2;rs=40;%choosehanningwindowaccordingtotable10.2

%basicparameterdw=ws-wp;M=ceil(3.11*pi/dw);N=M*2;%lengthis(N+1)

%designusingfir1h=fir1(N,wn/pi,hann(N+1));freqz(h,1);

figure;%designaccordingtoCHAP10wc=(wp+ws)/2;n=-M:M;hd=wc/pi*sinc(n*wc/pi);w=0.5+0.5*cos(2*pi*n/(N+1));ht=hd.*w;freqz(ht,1);freqz,real,imag,abs,angle,unwrap,phase%Givenasignal%x(t)=sin(0.1*pi*t)+2*cos(0.3*pi*t)+3*sin(0.5*pi*t)%whenusingasamplingfrequencyfT=20KHz,%plotthemagnitudeandphasespectrumofthesampledsequence%(givenlength-64).N=64;fT=20*1000;n=0:N-1;t=n/fT;x=sin(0.1*pi*t)+2*cos(0.3*pi*t)+3*sin(0.5*pi*t);[hw]=freqz(x,1);subplot(2,1,1);plot(w,abs(h));title('Magnitude');xlabel('\omega/\pi');ylabel('Magnitude');subplot(2,1,2);plot(w,angle(h));%angle()wraptheresultsto[-pi,pi]title('Phase');xlabel('\omega/\pi');ylabel('Phase');fir1,hamming窗%Usingthefunctionfir1andwindowofHamming,%designalinear-phaseFIRlowpassfilter%meetingthefollowingspecifications:%passbandedgefrequency=2kHz,%stopbandedgefrequency=2.3kHz,%passbandrippleδp=0.002,%stopbandrippleδs=0.002%andsamplingrateof20kHz.%Plotitsgainandphaseresponses%andcheckifitmeetsthespecifications?rp=0.002;%thetwoparameterisonlyusedto%determinewhichwindowtochoose%accordingtotable10.2rs=0.002;%=26.9897dB,hammingwindowcando

fp=2*1000;fs=2.3*1000;ft=20*1000;fc=(fp+fs)/2;wc=2*pi*fc/ft;dw=2*pi*(fs-fp)/ft;M=ceil(3.32*pi/dw);N=2*M;h=fir1(N,wc/pi);%defaulthammingwindowfreqz(h,1);直接型转零极点型，tf2zp%UsingMATLABprogramtoDetermineallpossible%ROCsofthefollowingz-Transforms,%anddescribethetypeoftheirinversez-Transforms%(left-sided,right-sided,two-sidedsequences)%associatedwitheachoftheROCs.%num=[8.16.93-23.8210.5];%den=[11.520.18-0.1768];

functionmainnum=[8.16.93-23.8210.5];den=[11.520.18-0.1768];getroc(num,den);return;

functiongetroc(num,den)

[z,p,k]=tf2zp(num,den);p=[0,p',inf];p=sort(unique(abs(p)));iflength(p)==2disp('ROC:entirez-plane');elsel=length(p);fprintf(1,'ROC:|z|<%f,leftsidedsignal\n',p(2));forii=2:l-2fprintf(1,'ROC:%f<|z|<%f,twosidedsignal\n',p(ii),p(ii+1));endfprintf(1,'ROC:|z|>%f,rightsidedsignal\n',p(l-1));end

return;butter,buttord%designanIIRbutterworthdigitallowpassfilter%withthefollowingspecifications:%samplingrateof7kHz,%passbandedgefrequencyof1.05kHz,%stopbandedgefrequencyof1.4kHz,%passbandrippleof0.4dB,%andaminimumstopbandattenuationof30dB.%plotfrequency-magnitudeand%checkifyourdesignfitsthespecification.ft=7*1000;fp=1.05*1000;fs=1.4*1000;rp=0.4;%dBrs=30;%dBwp=fp*2*pi/ft;ws=fs*2*pi/ft;[n,wn]=buttord(wp/pi,ws/pi,rp,rs);[b,a]=butter(n,wn);freqz(b,a);部分分式展开（逆），从传输函数计算冲击响应，residuez,impz%WriteaMATLABprogramtocomputethefirst%Lsamplesoftheinverseofrationalz-transforms%wherethevalueofLisprovidedbytheuser%throughthecommandinput.Usingthisprogram%tocomputeandplotthefirst50samplesof%theinverseoffollowingG(z).Usethecommand%stemforplottingthesequencegenerated%bytheinversetransform.%G(z)=-2+10/(4+z^(-1))-8/(2+z^(-1)),|z|>0.5L=input('L=');

r=[10/4,-8/2];p=[-1/4,-1/2];k=[-2];[b,a]=residuez(r,p,k);[h,t]=impz(b,a,L)stem(t,h);freqz,real,imag,abs,phase,angle,unwrap%TrytogiveaprogramtoevaluatethefollowingDTFT%intherange:(0,pi)%H(ejw):%num=[1.35,4.95,8.55,4.95,1.8];%den=[0.9,-1.8,1.65,-0.75,0.15];%ThenruntheProgramandcomputetherealandimaginary%partsoftheDTFT,andthemagnitudeandphasespectra.%IstheDTFTaperiodicfunctionofw?%Ifitis,whatistheperiod?.%Canyouexplainthejumpinthephasespectrum?%ThejumpcanberemovedusingtheMATLABcommand"unwrap".%EvaluatethephasespectrumwiththejumpremovedN=512;

num=[1.35,4.95,8.55,4.95,1.8];den=[0.9,-1.8,1.65,-0.75,0.15];w=(0:N-1)*pi/N;[H,w]=freqz(num,den);subplot(2,2,1);plot(w/pi,real(H));title('realpartofH');xlabel('\omega/\pi');ylabel('magnitude');subplot(2,2,2);plot(w/pi,imag(H));title('imaginepartofH');xlabel('\omega/\pi');ylabel('magnitude');subplot(2,2,3);plot(w/pi,abs(H));title('magnitudeofH');xlabel('\omega/\pi');ylabel('magnitude');subplot(2,2,4);plot(w/pi,phase(H))%unwrap(angle(H)));%equivalenttitle('phaseofH');xlabel('\omega/\pi');ylabel('phase');ellipord,ellip,freqz%designanellipticIIRlowpassfilter%withthefollowingspecifications:%samplingrateof7kHz,%passbandedgefrequencyof1.05kHz,%stopbandedgefrequencyof1.4kHz,%passbandrippleof0.4dB,%andaminimumstopbandattenuationof30dB.%plotfrequency-magnitudeand%checkifyourdesignfitsthespecificationft=7*1000;fp=1.05*1000;fs=1.4*1000;rp=0.4;rs=30;ws=fs*2*pi/ft;wp=fp*2*pi/ft;[n,wn]=ellipord(wp/pi,ws/pi,rp,rs);[b,a]=ellip(n,rp,rs,wn);[h,w]=freqz(b,a);plot(w/pi,log10(abs(h))*20);xlabel('\omega/\pi');ylabel('magnitude,dB');title('ellipticfilter');fir1,hann,分贝转换%DesignusingthewindowedFourierseriesapproach%alinear-phaseFIRlowpassfilterwiththefollowingspecifications:%passbandedgeat3.1rad/sec,%stopbandedgeat5rad/sec,%maximumpassbandattenuationof0.2dB,%minimumstopbandattenuationof45dB,%andasamplingfrequencyof40rad/sec.%UsethewindowofHanningtoperformit,%thenshowtheimpulseresponsecoefficients%andplotthegainresponseofthedesignedfilters.ft=40/(2*pi);wp=3.1/ft;ws=5/ft;rp=0.2;%notusebellowrs=45;%dBwn=(wp+ws)/2;wd=ws-wp;M=ceil(3.11*pi/wd)+1;N=2*M;h=fir1(N,wn/pi,hann(N+1))[Hw]=freqz(h,1);plot(w/pi,log10(abs(H))*20);title('gain');xlabel('\omega/\pi');ylabel('gain(dB)');信号能量计算%Generateandplotasequencex[n]=sin(5*pi/16*n),%with0<=n<=50.Computetheenergyofthesequence.n=0:50;x=sin(5*pi/16*n);stem(n,x);E=sum(x.*x)各IIR阶数估计，buttord,cheb1ord,cheb2ord,ellipord%UsingMATLABtodeterminethelowestorder%ofadigitalIIRlowpassfilterofallfourtypes%(Butterworth,Type1Cheby,Type2Cheby,elliptic).%Thespecificationsareasfollowing:%samplingrateof12kHz,%passbandedgefrequencyof2.1kHz,%stopbandedgefrequencyof2.7kHz,%passbandrippleof0.6dB,%andaminimumstopbandattenuationof45dB.%Commentonyourresults.

functionmain

ft=12*1000;fp=2.1*1000;fs=2.7*1000;rp=0.6;rs=45;wp=fp*2*pi/ft;ws=fs*2*pi/ft;iirord('buttord',wp,ws,rp,rs);iirord('cheb1ord',wp,ws,rp,rs);iirord('cheb2ord',wp,ws,rp,rs);iirord('ellipord',wp,ws,rp,rs);return;

functioniirord(fun,wp,ws,rp,rs)n=feval(fun,wp/pi,ws/pi,rp,rs);fprintf(1,'%s:%d\n',fun,n);return;square,freqz%Considerthefollowingperiodicsignalasbelow,%>>squarewave,period0.4s,dutycycle50%<<%whichissampledtoget512-pointsdiscrete-timesignal%ThesamplingperiodisT=0.02second.%Determinethesamplingtypeoftheprocess.%Sketchthetimedomainandfrequencydomainplot%ofcontinuousandsampledsignal.%Analysisthesamplingprocess.N=512;n=0:N-1;T=0.02;t=T*n;x=(square(2*pi/0.4*t)+1)/2;subplot(2,1,1);stem(t,x);title('timedomain');%ft=50,f0=2.5,sooversamplingsubplot(2,1,2);[hw]=freqz(x,1);plot(w,abs(h));title('frequencydomain');%itremainsaproblemhowtoplot%acontinuoussignal'sfrequencyspectrum线性卷积，此处用fft实现%WritingaMATLABprogramtocomputethe%linearconvolutionoftwolength-Nsequences.%Usingthisprogramtodeterminethe%followingpairofsequences:%g[n]=cos(0.2*pi*n);%h[n]=3.^n;%0<=n<=10;%Andplottheresultsequence.

functionmainn=0:10;g=cos(0.2*pi*n);h=3.^n;y=lconv(g,h);nn=0:length(y)-1;stem(nn,y);title('result');return;

functionre=lconv(x,y)N=length(x);M=length(y);L=N+M-1;X=fft(x,L);Y=fft(y,L);RE=X.*Y;re=ifft(RE);return;圆周卷积，要求用基于DFT的方法实现%WritingaMATLABprogramtocompute%thecircularconvolutionof%twolength-NsequencesviatheDFT-basedapproach.%Usingthisprogramtodeterminethefollowingpairofsequences:%g[n]={7,4,-9,0,2,-5},h[n]={1,-1,2,0,10,5}%org[n]=exp(j*pi*n/4),h[n]=2.^n,0<=n<=15%Andplottheresultsequence.

functionmain

g1=[7,4,-9,0,2,-5];h1=[1,-1,2,0,10,5];re1=cconv(g1,h1);stem(re1);

figure;n=0:15;g2=exp(j*pi*n/4);h2=2.^n;re2=cconv(g2,h2);stem(re2);

return;

functionre=cconv(x,y)N=length(x);X=fft(x,N);Y=fft(y,N);RE=X.*Y;re=ifft(RE);return;部分分式展开，逆z变换，residuez%WriteaMATLABprogramtodeterminethe%partial-fractionexpansionofthez-transform%listedthefollowingG(z)%num=[2,5,9,5,3];%den=[5,45,2,1,1];%|z|>1.2;%andthendeterminetheirinversez-transforms.num=[2,5,9,5,3];den=[5,45,2,1,1];[r,p,k]=residuez(num,den);

flag=0;

ftttorii=1:length(k)flag=1;fprintf(1,'(%s)*delta[n',num2str(k(ii)));ifii>1fprintf(1,'-%d]',ii-1);elsefprintf(1,']');endend

forii=1:length(r)ifflagfprintf(1,'+');elseflag=1;endfprintf(1,'(%s)*(%s).^n*u[n]',num2str(r(ii)),num2str(p(ii)));end

fprintf(1,'\n');freqz%WriteaMATLABprogramtocomputeand%plotthefrequencyresponseofacausalLTI%discrete-timesystemwithatransferfunctiongivenby%H(z)%num=[0.15,0,-0.15];%den=[1,-0.5,0.7];%for0≤ω≤π.Whattypeoffilterdoe