1、Advanced Digital Signal Processing(Modern Digital Signal Processing)Chapter 4 Power Spectrum Estimation Amplitude Spectrum&Power SpectrumnAmplitude spectrum densitynPower spectrum density(PSD)Finite-energy deterministic signal 4.1 IntroductionFourier transformAmplitude spectrum(including phase)Real
2、stationary random signal Fourier transformPower spectrum(without phase)Autocorrelation function Signal amplitude at each frequency Signal power at each frequency()()1()()2jj mxxxxmjj mxxxxPerm ermPeed xxrmE x n x nm Wiener-Khintchine TheoremFourier transform Inverse Fourier transform()*()21()lim()()
3、211lim()()211lim()()211lim()21Njj mxxNmnNNj njn mNnNmNj njn mNnNmNj nNnNPex n x nm eNx n ex nm eNx n ex nm eNx n eN1()lim()()()()21NxxNnNrmx n x nmE x n x nmN PSD of Ergodic Stationary Random Signal Power Spectrum Estimation(PSE)Estimating the PSD of real ergodic stationary random signal with finite
4、 observations(sample)Classical&Modern PSE MethodsnClassical(linear)PSE PSE based on Fourier transform,non-parametric model method nModern(nonlinear)PSE PSE based on signal model(parametric model method)01NNnxx Nwnnxx Desirable properties of PSEnUnbiasednConsistent nEfficient nHigh frequency resoluti
5、on Narrow main lobe,low side lobenSmall sample length 10111NmxxNNnNj mBTxxmNrmxn xnmNPrm e Blackman-Tukey(BT)Method for PSE 4.2 Classical PSEThe Fourier transform of the biased estimation of autocorrelation function.21()lim()21Njj nxxNnNPex n eN()NnxjxxPe Periodogram Method for PSE nPeriodogram 1NTh
6、e average energy spectrum of finite length sample is an estimation of PSD212011()()Njj njxxNNnPexn eXeNN2FFT()jNXe2()jNXe1(1)limlim1()Njj mjxxxxxxNNmNmE Perm ePeNnThe periodogram PSE is an asymptotically unbiased but not a consistent estimation2224sinlim varlim1sinjjxxxxNNjxxxNPePeNPeIf()is Gaussian
7、 white noise,thenx nBartlet window Limitations of Classical PSEnLow frequency resolution Caused by the effects of data window:(1)The degradation in resolution by main lobe;(2)The power leakage by side lobe(inter-spectrum interference).nInconsistent estimation00.20.40.60.81-40-30-20-10010Normalized F
8、requency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density Estimate via Periodogram00.20.40.60.81-40-30-20-10010Normalized Frequency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density Estimate via Periodogram00.20.40.60.81-40-30-20-1001020Normalized Frequency (rad/sample
9、)Power/frequency(dB/rad/sample)Power Spectral Density Estimate via Periodogram00.20.40.60.81-40-30-20-1001020Normalized Frequency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density Estimate via PeriodogramThe PSE of Harmonic Process With Periodogram96N 32N 256N 128N The PSE of White No
10、ise With Periodogram00.20.40.60.81-30-20-10010Normalized Frequency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density Estimate via Periodogram00.20.40.60.81-30-20-10010Normalized Frequency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density Estimate via Periodogram00.20.40
11、.60.81-30-20-10010Normalized Frequency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density Estimate via Periodogram00.20.40.60.81-30-20-10010Normalized Frequency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density Estimate via Periodogram64N 32N 256N 128N Modifications of C
12、lassical PSE Achieving low variance at the expense of bias and frequency resolutionnAveraging Periodogram(Bartlet method)AveragingN:data length,N=LMMMMMPeriodogramPeriodogramPeriodogramPeriodogramAveraging Periodogram1(1)1(1)1()1()MBjjj mxxixxmMNjj mxxxxmNmE PeE Ierm eMmE Perm eN1varvarvarvarjjBjjxx
13、ixxiPeIePeIeL21011()1()Mjj niNnLBjjxxiiIexMin eMPeI eL Its bias is larger than the periodogram while its variance is less than the Periodogram:1(1)1()()Mjj mxxxxmMmE Perm w m eNnModified Periodogram()1()()2jjjxxNPeIeW ed 210where()is the PSD of window(),and 1 ()jNjj nNNnW ew nIexn eN The window will
14、 smooth the PSD acquired by periodogram.Its function is similar to a lowpass filter.nAveraging modified periodogram(Welch method)AveragingN:data length,N=LMMMMMModified PeriodogramAveraging PeriodogramModified PeriodogramModified PeriodogramModified Periodogram00.20.40.60.81-40-30-20-100Normalized F
15、requency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density Estimate via Welch00.20.40.60.81-80-60-40-20020Normalized Frequency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density Estimate via Welch00.20.40.60.81-60-40-20020Normalized Frequency (rad/sample)Power/frequency(
16、dB/rad/sample)Power Spectral Density Estimate via Welch00.20.40.60.81-80-60-40-20020Normalized Frequency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density Estimate via WelchThe PSE of Harmonic Process With WelchThe sequence is divided into eight sections with 50%overlap,each section i
17、s windowed with a Hamming window96N 32N 256N 128N 00.20.40.60.81-30-20-10010Normalized Frequency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density Estimate via Welch00.20.40.60.81-30-20-10010Normalized Frequency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density Estimate
18、 via Welch00.20.40.60.81-30-20-10010Normalized Frequency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density Estimate via Welch00.20.40.60.81-30-20-10010Normalized Frequency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density Estimate via WelchThe PSE of White Noise With We
19、lchThe sequence is divided into eight sections with 50%overlap,each section is windowed with a Hamming window64N 32N 256N 128N 4.3 Parameter Model Methods for PSE Basic Principles Autocorrelation functionClassical PSEPSDFourier transformObservations xN(n)Estimation of signal model H(z)Parameter mode
20、l methodsPSDObservations xN(n)22()()jjxxwP eH e The Time Series Model of Stationary Random Signal 20,var()wwmw nLinear system with transfer function H(z)White noise w(n)Stationary random sequence x(n)1111111()1()()11qiiqqipppiiibzB zbzb zH zA za za za z22()()jjxxwPeH e21()()()xxwSzH z H znMA(q)model
21、(all-zero model)Suitable for signals whose power spectra have vales but no peaksnAR(p)model(all-pole model,most widely used)Suitable for signals whose power spectra have peaks but no vales,but be widely used since the linear relation between its parameters and the signal autocorrelation function 1()
22、()1qiiiH zB zb z 22111jxxwpj iiiPeae2211qjj ixxwiiPebe111()()1piiiH zA za znARMA(p,q)model(zero-pole model)Suitable for signals whose power spectra have vales and peaks221111qj iijixxwpj iiibePeae111()()()1qiiipiiibzB zH zA za znModel parameters to be estimated21,pwaa111()()()1qiiipiiibzB zH zA za z
23、211,pqwaabb21,qwbbAR(p):MA(q):ARMA(p,q):21 xxwSz A zB z H z111()()()1qiiipiiibzB zH zA za z 211200(),1pqxxk xxwkkkqwkkrma rmkhmb h kmb h kmb 2121 xxwwB zSzH z H zH zA z The Relation between the Autocorrelation Function&the Model Parameters nARMA(p,q)model Inverse z transform0 h kmkmIf()is causal,i.e
24、.()0 when 0H zh nn 00,0,1,0,0,1,0,qqkk mkkq mkmkb h kmmqb h kmmqbh kmqmq 201,0,1,0,q mpwk mkxxk xxkbh kmqrma rmkmq 2101,0,1,pq mk xxwk mkkxxpk xxka rmkbh kmqrma rmkmq Generalized Yule-Walker equations:a nonlinear equation set,but the equations are linear when mq.nMA(q)model 20,0,1,0,q mwk mkxxkbbmqr
25、mmqnAR(p)model0 1b 211,0,0pk xxwkxxpk xxka rmkmrma rmkm 2011(0),0,0pk xxwkxxpk xxka rmkb hmrma rmkmYule-Walker equation:a linear equation set 110limlim11pzzkkhH za k zInitial value theorem nAR model power spectrum estimation(AR PSE)211 ,0,1,pk xxwkxxpk xxka rmkmrma rmkmp 0,1,xxxxxxrrrpEstimation of
26、autocorrelation functionObservations xN(n)22()()jjxxwPeH eEstimation of AR model parametersAR model estimation ()H z 211,0 ,1,pxxwkxxpxxka k rmkmrma k rmkmp The AR model parameter estimation is obtained by solving the Yule-Walker equation(m=0,1,p):0,1,xxxxxxrrrp Properties of AR PSEnThe implied auto
27、correlation function extensionWith the p+1 samples of autocorrelation function(ACF)estimationFor mp,the can be extrapolated from those ACF estimations of mp by 1,pxxxxkrma k rmkmp i.e.,extrapolating from xxrm 0,1,xxxxxxrrrp ,xxrmmp tolniiiHpp Continuous value random variable()ln()ln()Hp xp x dxEp x
28、nMaximum entropy spectral estimation(MESE)&AR PSEnEntropyDiscrete value random variableEntropy of random variable XUncertainty of random variable XMaximum entropyPDF with maximum uncertainty(least restriction),xxrmmNKnown ACF estimation 0,1,xxxxxxrrrNnMESEMaximum entropy extrapolationUnknown ACF est
29、imationACFs with maximum uncertaintynMESE of zero-mean Gaussian random sequence PDF of N-dimension Gaussian random sequence 01,(0)(1)()(1)(0)(1)()()(1)(0)TNxxxxxxxxxxxxxxxxxxxxXx xxrrrNrrrNRNrNrNr11122121(,)(2)det()exp()2NTNxxxxp x xxRNXRNX1122ln(2)det()NxxHRN(0)(1)(1)(1)(0)()(1)(1)()(0)xxxxxxxxxxxx
30、xxxxxxxxrrrNrrrNRNrNrNr(1)(0)(1)(2)(1)(2)0(1)()(1)xxxxxxxxxxxxxxxxxxrrrNrrrNrNrNrmaxmax det()xxHRN(1)det(1)maxdet(1)0(1)xxxxxxrNxxdRNRNd rN(1)(0)(1)(2)(1)(2)0(1)()(1)xxxxxxxxxxxxxxxxxxrrrNrrrNrNrNr 0,1,xxxxxxrrrN(1)xxrN and 0,1,1xxxxxxrrrN (2)maxdet(2)xxxxrNRN(2)xxrN and so on.nThe equivalence betwe
31、en the AR PSE and the MESE of Gaussian random sequence AR PSE for p=N:121212 1011 2102 111xxxxxxN xxxxxxxxN xxxxxxxxN xxra ra ra rNra ra ra rNrNa rNa rNa r 211 ,0,0Nk xxwkxxNk xxka rmkmrma rmkm1,2,1mN(1)(0)(1)(2)(1)(2)0(1)()(1)xxxxxxxxxxxxxxxxxxrrrNrrrNrNrNr 121212 10110 21020 1110 xxxxxxN xxxxxxxxN
32、 xxxxxxxxN xxra ra ra rNra ra ra rNrNa rNa rNa rextrapolating(1)from(0),(1),()xxxxxxxxrNrrrNMESE of Gaussian random sequenceFor Gaussian random sequenceMaximum entropy extrapolationAR PSE implied extrapolation=MESEAR PSE=There are no poles of its AR model outside the unit circle,else 22()xxn2 lim()x
33、nn nThe stability of AR modelStationary random sequence x(n)The AR model of stationary random sequence is stable(minimum phase model)nThe relationship between the AR PSE and the linear prediction One-step pure linear optimal prediction filter(1)(0)()(1)(1)(1)(1)x nhx nhx nh px npAR model:12()(1)(2)(
34、)()px na x na x na x npw n(1)(1)()x nx ne n(1)(0)()(1)(1)(1)(1)()x nhx nhx nh px npe n(1)kh ka or12(1)()(1)(1)(1)px na x na x na x npw n()(1)e nw nhence,A one-step pure linear optimal prediction filter is the solution of the Yule-Walker equation:min11(0)(1)()(1)(0)(1)0()(1)(0)0 xxxxxxxxxxxxpxxxxxxrr
35、rparrrparprpr0here,1;(1),1,2,1kaah kkp 20,0 0,1,2,pwk xxkma rmkmpmin0,00,1,2,pk xxkma rmkmpAR(p)model The AR(p)parameters could be obtained as the coefficients that minimized the prediction error power of a p-th order linear predictor.21,1,11,1,11011101001011100 xxxxxxxxpwpxxxxxxxxppxxxxxxxxppxxxxxx
36、xxrrrprparrrprparprprrarprprr 2,1,1011010100 xxxxxxp wpxxxxxxp pxxxxxxrrrparrrparprpr Methods of AR PSE(Solutions of Y-K Equation)nLevinson-Durbin recursive algorithmp order AR model equationp+1 order AR model equation 2,1,1011010100 xxxxxxp wpxxxxxxp pxxxxxxrrrparrrparprpr Let 2,1,0111101001011100
37、xxxxxxxxp wxxxxxxxxpxxxxxxxxp pxxxxxxxxprrrprprrrprparprprrarprprrD,001,1ppp i xxpiDa rpia ,12,0110010101011101pxxxxxxxxxxxxxxxxp pxxxxxxxxpp wxxxxxxxxDrrrprprrrprparprprrarprprr 2,1,0111101001011100 xxxxxxxxp wxxxxxxxxpxxxxxxxxp pxxxxxxxxprrrprprrrprparprprrarprprrD 2,11,1000p wppp ppaRaD Let 11011
38、101101110 xxxxxxxxxxxxxxxxTppxxxxxxxxxxxxxxxxrrrprprrrprpRRrprprrrprprr,1,12,0001pp pppp wDaRa Expanded equationPreparative equation22,1,12,0000000pp wpwpp wpDD 2,1,111,12,10000001pp wpp ppppp ppp wpDaaRaaD1,1,1,11,11,111001ppp ppppp ppppaaaaaaaIf 21,1,111,1,11000pwppppppaRaa1,1,1pip ipp piaaa Then
39、i.e.1:reflection coefficientp22,1,12,0000000pp wpwpp wpDD 12,22221,1,1 1ppp wpwp wppp wpDD221,121,0 pwp wppppp wDD2221,1221,1pwp wppwp wThe predictive error is reduced gradually as p increases.2,21,1 and is known,andwhat is to be determined is and p wppwpD1,01,1(1)1;(0)xxxxraar 1,1,1pip ipp piaaa,00
40、1,1ppp i xxpiDa rpia 121,1,0 wi xxia ri12,22221,1,1 1ppp wpwp wppp wpDDLevinson-Durbin recursive algorithm2,until 1 reach at the predetermined order or k wpfrom 1p nBurgs recursive algorithmComputing AR model parameters from observation data directly 2000,(0)wxxenenx nr1,1,Bpip ipp p iaaa 11111,1BBp
41、pppppppenenenenenen 1111221211Nppn pBpNppn pen enenen 22221,1,11Bpwp wpkp wpDfor 0,1,2,1 pp00.20.40.60.81-400-300-200-1000Normalized Frequency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density Estimate via Burg00.20.40.60.81-400-300-200-1000Normalized Frequency (rad/sample)Power/frequ
42、ency(dB/rad/sample)Power Spectral Density Estimate via Burg00.20.40.60.81-400-300-200-1000Normalized Frequency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density Estimate via Burg00.20.40.60.81-350-300-250-200-150-100-50Normalized Frequency (rad/sample)Power/frequency(dB/rad/sample)Pow
43、er Spectral Density Estimate via BurgThe PSE of Harmonic Process With Burgs Method96N 32N 256N 128N 00.20.40.60.81-30-20-10010Normalized Frequency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density Estimate via Burg00.20.40.60.81-30-20-10010Normalized Frequency (rad/sample)Power/freque
44、ncy(dB/rad/sample)Power Spectral Density Estimate via Burg00.20.40.60.81-30-20-10010Normalized Frequency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density Estimate via Burg00.20.40.60.81-30-20-10010Normalized Frequency (rad/sample)Power/frequency(dB/rad/sample)Power Spectral Density E
45、stimate via BurgThe PSE of White Noise With Burgs Method64N 32N 256N 128N Selection of AR Models OrdernBasic principles 2,p wp2,is a balance between the and the total parameter estimation errorp wpmodel parameterstotal parameter estimation errorp Higher order has lower prediction error or higher app
46、roximation accuracy nCriteria to select AR model ordernFinal prediction error(FPE)criterionnAkaike information criterion(AIC)nCriterion autoregressive transfer(CAT)nExperiential method2,min()minp wppNpFPE pNp21,111min()minpppiip wCAT pN2,min()minln2p wppAIC pNp2,ii wNNi32NNp MA&ARMA PSEnMA PSE Durbins method:based on AR PSEnARMA PSE 1,1,2,pxxk xxkrma rmkmqqqp 212,pwb bb12,pa aaDurbins method Textbook：n习题：2、3、5、6、7、8n上机作业：1、2Exercises
