DOC PREVIEW
RIT SIMG 713 - Power Spectrum Estimation

This preview shows page 1-2-24-25 out of 25 pages.

Save
View full document
View full document
Premium Document
Do you want full access? Go Premium and unlock all 25 pages.
Access to all documents
Download any document
Ad free experience
View full document
Premium Document
Do you want full access? Go Premium and unlock all 25 pages.
Access to all documents
Download any document
Ad free experience
View full document
Premium Document
Do you want full access? Go Premium and unlock all 25 pages.
Access to all documents
Download any document
Ad free experience
View full document
Premium Document
Do you want full access? Go Premium and unlock all 25 pages.
Access to all documents
Download any document
Ad free experience
Premium Document
Do you want full access? Go Premium and unlock all 25 pages.
Access to all documents
Download any document
Ad free experience

Unformatted text preview:

Power Spectrum EstimationLecture 15Spring 2002Wiener-Khintchine TheoremLet x(n) be a WSS random process with autocorrelation sequencerxx(m)=E[x(n + m)x∗(n)]The power spectral density is defined as the Discrete Time FourierTransform of the autocorrelation sequencePxx(f)=T∞n=−∞rxx(m)e−i2πfmTwhere T is the sampling interval.The signal is assumed to be bandlimited in frequency to ±1/2T andis periodic in frequency with period 1/T .Lecture 15 1Wiener-Khintchine TheoremThe inverse DTFT isrxx(m)=1/2T−1/2TPxx(f)ei2πfmTdfAnd, rxx(0) is the average powerrxx(0) =1/2T−1/2TPxx(f)dfDue to the property rxx(−m)=r∗xx(m), the PSD must be a strictlyreal, nonnegative function.Lecture 15 2White NoiseA white noise process is zero for all lags except m = 0. Thus,rww(m)=σ2wδ(m)The PSD is thereforePww(f)=σ2wTIf y(n) is the sequence that is produced by exciting a discrete linearsystem with white noise, thenPyy(f)=σ2wH(z)H(1/z) z=ei2πfTLecture 15 3Ergodic Random ProcessesIf a process is ergodic then ensemble averages can be replaced with“time” averages. We will assume ergodic sequences from here on.Power spectrum definition using time averages:Pxx(f) = limm→∞E1(2M +1)T TMn=−Mx(n)e−i2πfnT 2This is the basis of the periodogram method for estimating the PSD.The Fourier transform of x(n) is computed, which can be computedvia the FFT. To reduce the variance in the estimate, many spectracan be averaged.Lecture 15 4Correlogram EstimateBegin with the definition:Pxx(f) = limm→∞E1(2M +1)T TMn=−Mx(n)e−i2πfnT 2After some algebra this can be converted to a DTFT of the auto-correlation sequencePxx(f)=T∞n=−∞rxx(m)e−i2πfmTThe Correlogram method of estimating the PSD first estimates theautocorrelation sequence and then transforms it to estimate thePSD.Lecture 15 5Autocorrelation Sequence EstimationAssume that n data samples indexed n =0toN − 1 are available.A possible estimator of the ACS isˆrxx(m)=1N − mN− m−1n=0x(n + m)x∗(n)The values of m for which the ACS is computed are know as the“lag” indexes.This estimator is unbiased, sinceE[ˆrxx(m)] =1N − mN− m−1n=0E[x(n + m)x∗(n)] = rxx(m)Lecture 15 6ACS EstimationThe variance is approximatelyvar[ˆrxx(m)] ≈N(N − m)2∞−∞(r2xx(k)+rxx(k + m)rxx(k − m))for N  m. Note that the variance increases for lags close to N.Lecture 15 7Biased ACS EstimatorAn alternative estimate isˇrxx(m)=1NN− m−1n=0x(n + m)x∗(n)This estimate is related to ˆrxx(m)byˇrxx(m)=N − mNˆrxx(m)It is a biased estimator withE[ˇrxx(m)=N − mNrxx(m)and variancevar[ˇrxx(m)] =N − mNvar[ˆrxx(m)]Lecture 15 8Numerical IssuesThe biased estimator can never produce an autocorrelation matrixthat is not positive semi-definite, whereas the biased estimator some-times does.Hence, there are numerical advantages to the biased estimator insome cases.Lecture 15 9Correlogram Estimate of PSDAssume that the ACS has been estimated for lags −L ≤ m ≤ L.Then the PSD can be estimated byˆPxx(f)=TLm=−Lˆrxx(m)e−i2πfmTThe value of the maximum lag index L should be much less thanthe number of data samples, N. Typically, L ≈ N/10.The biased estimator can be substituted for the unbiased estimator.Lecture 15 10Correlogram ProgramFUNCTION CORRELOGRAMPSD,X,Y,LAG,T,NF,R;+; PSD=CORRELOGRAMPSD(X,Y,LAG,T,NF); computes the power spectrum density by; the Blackman-Tukey correlogram method.;; INPUTS; X - A Sequence of length NX; Y - A Sequence of length NY; LAG - The maximum lag to be used.; T - The interval between samples; NF - The number of frequency points;; OUTPUT; PSD = Power Spectral Density estimate at NF points;; Notes:; The length of the correlation sequence is; N=NX<NY. The minimum is used.;; H. Rhody; May, 1998;; REFERENCE; S. L. Marple,; Digital Spectral Analysis with Applications,; Prentice Hall, 1987; Appendix 5.B;-Correlogram program, continuedlags=FINDGEN(LAG)+1lags=[lags-LAG-1,0,lags]w=0.538+0.462*cos(!pi*lags/LAG)NX=N_ELEMENTS(X)NY=N_ELEMENTS(Y)N=NX<NYR=C_CORRELATE(X[0:N-1],Y[0:N-1],lags,/COVARIANCE)R=R*N/(N-ABS(lags))R=R*wS=[R[LAG:2*LAG],FLTARR(NF-2*LAG-1),R[0:LAG-1]]PSD=ABS(FFT(S))RETURN,SHIFT(PSD,NF/2)ENDLecture 15 12Periodogram PSD EstimatorThe periodogram is based on the DTFT of the sample sequence.From the definition,Pxx(f) = limm→∞E1(2M +1)T TMn=−Mx(n)e−i2πfnT 2we can form an estimator like˜P(f)=TN N− 1n=0x(n)e−i2πfnT 2Lecture 15 13Improved Periodogram PSD EstimatorThis estimator can be improved by• Smoothing the sample sequence before computing the transform• Segmenting the sample sequence and computing more estima-tors that can be averaged• Using overlapping sample sequencesThese techniques can be combined to construct the Welch Peri-odogram estimator.Lecture 15 14Periodogram Estimator ProgramFUNCTION PERIODOGRAM,X,Y,D,S,WINDOW=window,WINFUN=w;+; PXY=PERIODOGRAM(X,Y,D,S); computes the cross power spectral density between; sequences X and Y. If X and Y are the same sequence; then a direct PSD is computed.;;; INPUT PARAMETERS; X An array of data samples; Y An array of data samples; D The size of the analysis window; S The shift of the analysis window;; The window function is determined by the keyword WINDOW:; WINDOW=0 (Default) All values equally weighted (no window); WINDOW=1 Hamming Window; WINDOW=2 Bartlett Window; WINDOW=3 Hann Window; WINDOW=4 Truncated Gaussian;Lecture 15 15Periodogram Estimator Program (continued); The window function can be accessed via WINFUN keyword;; If X and Y are not the same length then the shorter; one sets the length for the analysis.;; Example: Find the PSD of a sequence X at 500 points; using a window of size 80 and shift of size 40.;; PSD=PERIODOGRAM(X,X,500,80,Window=1); PLOT,PSD; Note that the PSD plot contains 500 points.;; H. Rhody May, 1998; Revised April 2000 to include window functions and; remove extraneous parameter NF. Added access to; window function via WINFUN.;; Reference:; S. L. Marple,; Digital Spectral Analysis with Applications,; Prentice Hall, 1987; Appendix 5.C;-Lecture 15 16Periodogram Estimator Program (continued)PSD=FLTARR(D)XS=FLTARR(D)YS=FLTARR(D)IF NOT KEYWORD_SET(window) THEN window=0k=(findgen(D)-(D-1)/2.)/(D-1)CASE window OF0: w=FLTARR(D)+1 ; No window. Default1: BEGIN ; Hamming Windoww=0.538+0.462*COS(2*!pi*k)END2: BEGIN ; Bartlett


View Full Document

RIT SIMG 713 - Power Spectrum Estimation

Download Power Spectrum Estimation
Our administrator received your request to download this document. We will send you the file to your email shortly.
Loading Unlocking...
Login

Join to view Power Spectrum Estimation and access 3M+ class-specific study document.

or
We will never post anything without your permission.
Don't have an account?
Sign Up

Join to view Power Spectrum Estimation 2 2 and access 3M+ class-specific study document.

or

By creating an account you agree to our Privacy Policy and Terms Of Use

Already a member?