现代数字信号处理AdvancedDigitalSignalProcessingch4PSE

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 phas

2、e)real 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 ()*()

3、21()lim( ) ()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 r

4、andom signal with finite 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 psenunbiasednconsist

5、ent nefficient nhigh frequency resolution 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

6、 periodogram method for pse nperiodogram 1nthe 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 estimation2224sinli

7、m varlim1sinjjxxxxnnjxxxnpepenpeif ( ) is gaussian 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).ninconsis

8、tent estimation00.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-10010normalized frequency ( rad/sample)power/frequency (db/rad/sample)power spectral density estimate via periodogram00.

9、20.40.60.81-40-30-20-1001020normalized 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)power/frequency (db/rad/sample)power spectral density estimate via periodogramthe pse of harmon

10、ic process with periodogram96n 32n 256n 128n the pse of white noise 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 (d

11、b/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.60.81-30-20-10010normalized frequency ( rad/sample)power/frequency (db/rad/sample)power spect

12、ral density estimate via periodogram64n 32n 256n 128n modifications of classical pse achieving low variance at the expense of bias and frequency resolutionnaveraging periodogram (bartlet method) averagingn: data length, n=lmmmmmperiodogramperiodogramperiodogramperiodogramaveraging periodogram1(1)1(1

13、)1( )1( )mbjjj mxxixxmmnjj mxxxxmnme pee ierm emme perm en1varvarvarvarjjbjjxxixxipeiepeiel21011()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 e

14、d 210where () is the psd of window ( ), and 1 ( )jnjj nnnnw ew niexn en the window will 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 periodogrammodifie

15、d periodogrammodified periodogrammodified periodogram00.20.40.60.81-40-30-20-100normalized frequency ( 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 densi

16、ty estimate via welch00.20.40.60.81-60-40-20020normalized frequency ( 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 welchthe pse of h

17、armonic process with welchthe sequence is divided into eight sections with 50% overlap, each section is 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-3

18、0-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 welch00.20.40.60.81-30-20-10010normalized frequency ( rad

19、/sample)power/frequency (db/rad/sample)power spectral density estimate via welchthe pse of white noise with welchthe 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 autocorr

20、elation functionclassical psepsdfourier transformobservations xn(n)estimation of signal model h(z)parameter model 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 rando

21、m sequence x(n)1111111( )1 ( )( )11qiiqqipppiiibzb zb zb zh za za za za z22()()jjxxwpeh e21( )( )()xxwszh z h znma(q) model (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 pe

22、aks but no vales, but be widely used since the linear relation between its parameters and the signal autocorrelation function 1 ( )( )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

23、 peaks221111qj iijixxwpj iiibepeae111( ) ( )( )1qiiipiiibzb zh za za znmodel parameters to be estimated21,pwaa111( ) ( )( )1qiiipiiibzb zh za za z211,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 xx

24、wwb 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. ( )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 m

25、k 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 mkxxkbbmqrmmqnar (p) model0 1b 211,0, 0pk xxwkxxpk xxka rmkmrma rmkm 2011(0),0, 0pk xxwkxxpk xxka rmkb hmrma rmkmyule-walker equation

26、: 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 autocorrelation functionobservations xn(n)22()()jjxxwpeh eestimation of ar model parametersar model estimation (

27、)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 autocorrelation function extensionwith the p+1 samples of autocorrelation function (acf) estimationfor mp, the

28、 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 nmaximum entropy spectral estimation (mese) & ar psenentropydiscrete value random variable

29、entropy of random variable xuncertainty of random variable xmaximum entropypdf with maximum uncertainty (least restriction) , xxrmmnknown acf estimation 0 ,1 ,xxxxxxrrrnnmesemaximum entropy extrapolationunknown acf estimationacfs with maximum uncertaintynmese of zero-mean gaussian random sequence pd

30、f 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)xxxxxxxxxxxxxxxxxxxxrrrnrrrnrnrnrnr(1)(0)(1)(2)(1)(2)0(1)()(1)xxxxxxxxxxxxxxxxxxrrrnrrrnrn

31、rnrmaxmax 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 between the ar pse and the mese of gaussian random sequence ar pse for p=n: 1

32、21212 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 xxxxxxxxn xxra ra ra rnra ra ra rnrna rna rna rextrapolating (1) fro

33、m (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( )xnn nthe stability of ar modelstationary random sequence x(n)the

34、 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)()( )px na x na x na x npw n(1)(1)( )x nx ne n(1)(0) ( )(1

35、) (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 xxxxxxxxxxxxpxxxxxxrrrparrrparprpr0here, 1; (1), 1,2,1kaah kkp 20,

36、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 xxxxxxxxpwpxxxxxxxxppxxxxxxxxppxxxxxxxxrrrprparrrprparprprrarprprr 2,1,1011

37、010100 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 xxxxxxxxp wxxxxxxxxpxxxxxxxxp pxxxxxx

38、xxprrrprprrrprparprprrarprprrd ,001, 1ppp i xxpida rpia ,12,0110010101011101pxxxxxxxxxxxxxxxxp pxxxxxxxxpp wxxxxxxxxdrrrprprrrprparprprrarprprr 2,1,0111101001011100 xxxxxxxxp wxxxxxxxxpxxxxxxxxp pxxxxxxxxprrrprprrrprparprprrarprprrd 2,11,1000p wppp pparad let 11011101101110 xxxxxxxxxxxxxxxxtppxxxxxx

39、xxxxxxxxxxrrrprprrrprprrrprprrrprprr,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 i.e. 1:reflection coefficientp22,1,

40、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 ,001, 1ppp i xxpida rpia 121,1,0

41、 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 , 1bbppppppppenenenenenen 11112

42、21211nppn 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/frequency (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-50normaliz