数字信号处理基于计算机的方法第三版matlab程序.docx
《数字信号处理基于计算机的方法第三版matlab程序.docx》由会员分享，可在线阅读，更多相关《数字信号处理基于计算机的方法第三版matlab程序.docx（23页珍藏版）》请在冰豆网上搜索。
数字信号处理基于计算机的方法第三版matlab程序
%Program2_1
%Generationoftheensembleaverage
%
R=50;
m=0:
R-1;
s=2*m.*（0.9.^m）;%Generatetheuncorruptedsignal
d=rand（R,1）-0.5;%Generatetherandomnoise
x1=s+d';
stem（m,d）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;title（'Noise'）;
pause
forn=1:
50;
d=rand（R,1）-0.5;
x=s+d';
x1=x1+x;
end
x1=x1/50;
stem（m,x1）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;title（'Ensembleaverage'）;
%Program2_2
%Generationofcomplexexponentialsequence
%
a=input（'Typeinrealexponent='）;
b=input（'Typeinimaginaryexponent='）;
c=a+b*i;
K=input（'Typeinthegainconstant='）;
N=input（'Typeinlengthofsequence='）;
n=1:
N;
x=K*exp（c*n）;%Generatethesequence
stem（n,real（x））;%Plottherealpart
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
title（'Realpart'）;
disp（'PRESSRETURNforimaginarypart'）;
pause
stem（n,imag（x））;%Plottheimaginarypart
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
title（'Imaginarypart'）;
%Program2_3
%Generationofrealexponentialsequence
%
a=input（'Typeinargument='）;
K=input（'Typeinthegainconstant='）;
N=input（'Typeinlengthofsequence='）;
n=0:
N;
x=K*a.^n;
stem（n,x）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
title（['\alpha=',num2str（a）]）;
%Program2_4
%SignalSmoothingbyaMoving-AverageFilter
%
R=50;
d=rand（R,1）-0.5;
m=0:
1:
R-1;
s=2*m.*（0.9.^m）;
x=s+d';
plot（m,d,'r-',m,s,'b--',m,x,'g:
'）
xlabel（'Timeindexn'）;ylabel（'Amplitude'）
legend（'d[n]','s[n]','x[n]'）;
pause
M=input（'Numberofinputsamples='）;
b=ones（M,1）/M;
y=filter（b,1,x）;
plot（m,s,'r-',m,y,'b--'）
legend（'s[n]','y[n]'）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）
%Program2_5
%IllustrationofMedianFiltering
%
N=input（'MedianFilterLength='）;
R=50;a=rand（1,R）-0.4;
b=round（a）;%Generateimpulsenoise
m=0:
R-1;
s=2*m.*（0.9.^m）;%Generatesignal
x=s+b;%Impulsenoisecorruptedsignal
y=medfilt1（x,3）;%Medianfiltering
subplot（2,1,1）
stem（m,x）;axis（[050-18]）;
xlabel（'n'）;ylabel（'Amplitude'）;
title（'ImpulseNoiseCorruptedSignal'）;
subplot（2,1,2）
stem（m,y）;
xlabel（'n'）;ylabel（'Amplitude'）;
title（'OutputofMedianFilter'）;
%Program2_6
%IllustrationofConvolution
%
a=input（'Typeinthefirstsequence='）;
b=input（'Typeinthesecondsequence='）;
c=conv（a,b）;
M=length（c）-1;
n=0:
1:
M;
disp（'outputsequence='）;disp（c）
stem（n,c）
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
%Program2_7
%ComputationofCross-correlationSequence
%
x=input（'Typeinthereferencesequence='）;
y=input（'Typeinthesecondsequence='）;
%Computethecorrelationsequence
n1=length（y）-1;n2=length（x）-1;
r=conv（x,fliplr（y））;
k=（-n1）:
n2';
stem（k,r）;
xlabel（'Lagindex'）;ylabel（'Amplitude'）;
v=axis;
axis（[-n1n2v（3:
end）]）;
%Program2_8
%ComputationofAutocorrelationofa
%NoiseCorruptedSinusoidalSequence
%
N=96;
n=1:
N;
x=cos（pi*0.25*n）;%Generatethesinusoidalsequence
d=rand（1,N）-0.5;%Generatethenoisesequence
y=x+d;%Generatethenoise-corruptedsinusoidalsequence
r=conv（y,fliplr（y））;%Computethecorrelationsequence
k=-28:
28;
stem（k,r（68:
124））;
xlabel（'Lagindex'）;ylabel（'Amplitude'）;
%Program3_1
%Discrete-TimeFourierTransformComputation
%
%Readinthedesirednumberoffrequencysamples
k=input（'Numberoffrequencypoints='）;
%Readinthenumeratoranddenominatorcoefficients
num=input（'Numeratorcoefficients='）;
den=input（'Denominatorcoefficients='）;
%Computethefrequencyresponse
w=0:
pi/（k-1）:
pi;
h=freqz（num,den,w）;
%Plotthefrequencyresponse
subplot（2,2,1）
plot（w/pi,real（h））;grid
title（'Realpart'）
xlabel（'\omega/\pi'）;ylabel（'Amplitude'）
subplot（2,2,2）
plot（w/pi,imag（h））;grid
title（'Imaginarypart'）
xlabel（'\omega/\pi'）;ylabel（'Amplitude'）
subplot（2,2,3）
plot（w/pi,abs（h））;grid
title（'MagnitudeSpectrum'）
xlabel（'\omega/\pi'）;ylabel（'Magnitude'）
subplot（2,2,4）
plot（w/pi,angle（h））;grid
title（'PhaseSpectrum'）
xlabel（'\omega/\pi'）;ylabel（'Phase,radians'）
%Program3_2
%Generatethefiltercoefficients
h1=ones（1,5）/5;h2=ones（1,14）/14;
%Computethefrequencyresponses
[H1,w]=freqz（h1,1,256）;
[H2,w]=freqz（h2,1,256）;
%Computeandplotthemagnituderesponses
m1=abs（H1）;m2=abs（H2）;
plot（w/pi,m1,'r-',w/pi,m2,'b--'）;
ylabel（'Magnitude'）;xlabel（'\omega/\pi'）;
legend（'M=5','M=14'）;
pause
%Computeandplotthephaseresponses
ph1=angle（H1）*180/pi;ph2=angle（H2）*180/pi;
plot（w/pi,ph1,w/pi,ph2）;
ylabel（'Phase,degrees'）;xlabel（'\omega/\pi'）;
legend（'M=5','M=14'）;
%Program3_3
%Setupthefiltercoefficients
b=[-6.7619513.456335-6.76195];
%Setinitialconditionstozerovalues
zi=[00];
%Generatethetwosinusoidalsequences
n=0:
99;
x1=cos（0.1*n）;
x2=cos（0.4*n）;
%Generatethefilteroutputsequence
y=filter（b,1,x1+x2,zi）;
%Plottheinputandtheoutputsequences
plot（n,y,'r-',n,x2,'b--',n,x1,'g-.'）;grid
axis（[0100-1.24]）;
ylabel（'Amplitude'）;xlabel（'Timeindexn'）;
legend（'y[n]','x2[n]','x1[n]'）
%Program4_1
%4-thOrderAnalogButterworthLowpassFilterDesign
%
formatlong
%Determinezerosandpoles
[z,p,k]=buttap（4）;
disp（'Polesareat'）;disp（p）;
%Determinetransferfunctioncoefficients
[pz,pp]=zp2tf（z,p,k）;
%Printcoefficientsindescendingpowersofs
disp（'Numeratorpolynomial'）;disp（pz）
disp（'Denominatorpolynomial'）;disp（real（pp））
omega=[0:
0.01:
5];
%Computeandplotfrequencyresponse
h=freqs（pz,pp,omega）;
plot（omega,20*log10（abs（h）））;grid
xlabel（'Normalizedfrequency'）;ylabel（'Gain,dB'）;
%Program4_2
%ProgramtoDesignButterworthLowpassFilter
%
%Typeinthefilterorderandpassbandedgefrequency
N=input（'Typeinfilterorder='）;
Wn=input（'3-dBcutoffangularfrequency='）;
%Determinethetransferfunction
[num,den]=butter（N,Wn,'s'）;
%Computeandplotthefrequencyresponse
omega=[0:
200:
12000*pi];
h=freqs（num,den,omega）;
plot（omega/（2*pi）,20*log10（abs（h）））;
xlabel（'Frequency,Hz'）;ylabel（'Gain,dB'）;
%Program4_3
%ProgramtoDesignType1ChebyshevLowpassFilter
%
%Readinthefilterorder,passbandedgefrequency
%andpassbandrippleindB
N=input（'Order='）;
Fp=input（'PassbandedgefrequencyinHz='）;
Rp=input（'PassbandrippleindB='）;
%Determinethecoefficientsofthetransferfunction
[num,den]=cheby1（N,Rp,2*pi*Fp,'s'）;
%Computeandplotthefrequencyresponse
omega=[0:
200:
12000*pi];
h=freqs（num,den,omega）;
plot（omega/（2*pi）,20*log10（abs（h）））;
xlabel（'Frequency,Hz'）;ylabel（'Gain,dB'）;
%Program4_4
%ProgramtoDesignEllipticLowpassFilter
%
%Readinthefilterorder,passbandedgefrequency,
%passbandrippleindBandminimumstopband
%attenuationindB
N=input（'Order='）;
Fp=input（'PassbandedgefrequencyinHz='）;
Rp=input（'PassbandrippleindB='）;
Rs=input（'MinimumstopbandattenuationindB='）;
%Determinethecoefficientsofthetransferfunction
[num,den]=ellip（N,Rp,Rs,2*pi*Fp,'s'）;
%Computeandplotthefrequencyresponse
omega=[0:
200:
12000*pi];
h=freqs（num,den,omega）;
plot（omega/（2*pi）,20*log10（abs（h）））;
xlabel（'Frequency,Hz'）;ylabel（'Gain,dB'）;
functiony=circonv（x1,x2）
%Developsasequenceyobtainedbythecircular
%convolutionoftwoequal-lengthsequencesx1andx2
L1=length（x1）;
L2=length（x2）;
ifL1~=L2,error（'Sequencesofunequallength'）,end
y=zeros（1,L1）;
x2tr=[x2
（1）x2（L2:
-1:
2）];
fork=1:
L1,
sh=circshift（x2tr',k-1）';
h=x1.*sh;
disp（sh）;
y（k）=sum（h）;
end
function[Y]=haar_1D（X）
%ComputestheHaartransformoftheinputvectorX
%ThelengthofXmustbeapower-of-2
%RecursivelybuildstheHaarmatrixH
v=log2（length（X））-1;
H=[11;1-1];
form=1:
v,
A=[kron（H,[11]）;
2^（m/2）.*kron（eye（2^m）,[1-1]）];
H=A;
end
Y=H*double（X（:
））;
%Program5_1
%IllustrationofDFTComputation
%
%ReadinthelengthNofsequenceandthedesired
%lengthMoftheDFT
N=input（'Typeinthelengthofthesequence='）;
M=input（'TypeinthelengthoftheDFT='）;
%Generatethelength-Ntime-domainsequence
u=[ones（1,N）];
%ComputeitsM-pointDFT
U=fft（u,M）;
%Plotthetime-domainsequenceanditsDFT
t=0:
1:
N-1;
stem（t,u）
title（'Originaltime-domainsequence'）
xlabel（'Timeindexn'）;ylabel（'Amplitude'）
pause
subplot（2,1,1）
k=0:
1:
M-1;
stem（k,abs（U））
title（'MagnitudeoftheDFTsamples'）
xlabel（'Frequencyindexk'）;ylabel（'Magnitude'）
subplot（2,1,2）
stem（k,angle（U））
title（'PhaseoftheDFTsamples'）
xlabel（'Frequencyindexk'）;ylabel（'Phase'）
%Program5_2
%IllustrationofIDFTComputation
%
%ReadinthelengthKoftheDFTandthedesired
%lengthNoftheIDFT
K=input（'TypeinthelengthoftheDFT='）;
N=input（'TypeinthelengthoftheIDFT='）;
%Generatethelength-KDFTsequence
k=0:
K-1;
V=k/K;
%ComputeitsN-pointIDFT
v=ifft（V,N）;
%PlottheDFTanditsIDFT
stem（k,V）
xlabel（'Frequencyindexk'）;ylabel（'Amplitude'）
title（'OriginalDFTsamples'）
pause
subplot（2,1,1）
n=0:
N-1;
stem（n,real（v））
title（'Realpartofthetime-domainsamples'）
xlabel（'Timeindexn'）;ylabel（'Amplitude'）
subplot（2,1,2）
stem（n,imag（v））
title（'Imaginarypartofthetime-domainsamples'）
xlabel（'Timeindexn'）;ylabel（'Amplitude'）
%Program5_3
%NumericalComputationofFouriertransformUsingDFT
升级会员 