How to get Autocorrelation from power spectrum in MATLAB? 2

[arhabib]

I'd like to compute the autocorrelation function of a time series from its power spectrum.
I am using MATLAB's pwelch function with a blackman window and 50% overlap.
How should I inverse the output from pwelch to get the autocorrelation?
sample code:

[px fx]=pwelch(x,blackman(512),256,512,100);
acfx=real(ifft(px*100));

The output from the ifft function has values other than 1.0 at t_lag=0!
I don't have an electrical background and little info about the math behind
this!
What kind of scaling does pwelch use?
Thanks in advance.

 

[LiteYear]

This is a pretty intricate calculation to be doing - in essence the autocorrelation and the PSD are fourier pairs, but to mathematically ensure that ac = ifft(psd) would require some careful conditioning. I think you'd be better off writing your own psd function so you can precisely control the scaling and treatment. Accounting for the intimate details of the built-in pwelch might be a frustrating experience. If you do manage it though, I'd be interested in seeing the details!

 

[GregLocock]

I may have picked the wrong end of the snake up, but surely a psd has destroyed all time relevant information (the phase) and therefore the autocorrelation will be meaningless.


Cheers

Greg Locock


New here? Try reading these, they might help FAQ731-376

http://eng-tips.com/market.cfm?

 

[LiteYear]

Greg, I think that's almost true. If you define the PSD as |F(w)|^2 then certainly the phase information gets lost in the modulus operation and the autocorrelation is unrecoverable.

More precisely however, the PSD is defined on infinite stationary processes as lim(T->inf) (1/T) * E(|F(w)|^2) where E is the expectation operator. This expands to:

lim(T->inf) (1/T) * E(F(w)*F'(w))

where F' is the complex conjugate of F.

Now provided the original process was real, F(w) = integral(f(t)*e^(-2*pi*i*w*t), w.r.t. t) and F'(w) = integral(f(s)*e^(2*pi*i*w*s), w.r.t. s)

So F(w)*F'(w) is just the double integral of f(t)*f(s)*e^(-2*pi*i*w*(t-s))

If you're really careful with the limits and the double integral, you can change variables to tau and another arbitrary variable, evaluate the arbitrary variable integral and arrive at

F(w)*F'(w) = integral( (T-tau) * f(t)*f(t-tau)*e^(-2*pi*i*w*(tau), w.r.t. tau)

If you apply the expectation operator, it only applies to the f(t)*f(t-tau) term and if you multiply by (1/T) you get:

(1/T) * E(F(w)*F'(w)) = integral( (1-tau/T) * E(f(t)*f(t-tau))*e^(-2*pi*i*w*(tau), w.r.t. tau)

Of course if we now take the limit as T->inf the tau/T term disappears the right hand side is just the Fourier transform of E(f(t)*f(t-tau)).

If we do a bit more handwaving, E(f(t)*f(t-tau)) is the autocorrelation of f(t).

So if all these conditions are true and you squint hard enough, the PSD of f(t) is equal to the Fourier transform of the autocorrelation of f(t). In fact, this is sometimes how it is defined. Theoretically then, the inverse fourier transform of the PSD gives you back the autocorrelation.

But I think by the time you convert to discrete, limit to a finite signal, and incorporate the fancy estimation stuff in MATLAB's pwelch function, this tenuous link to the autocorrelation might be a bit hard to recover!

PS. I'm sorry the math notation is going to be very hard to follow. I tried to make it as clean as possible. It is possible to render LaTeX or something instead?

 

[GregLocock]

Thanks- I think I understand the general idea, and as you say it is going to get a bit hairy once you start trying to solve it arithmetically rather than symbolically. star for trying to do real maths in a text box!

Cheers

Greg Locock


New here? Try reading these, they might help FAQ731-376

http://eng-tips.com/market.cfm?

 

[hacksaw]

Liteyear,

that's how it was defined for us in the day, where the FT of the autocorrelation was the recommended way of reducing the discrete time events to PSD, we later used the FT directly, getting the same result as long as we took care to insure that the avg value of the time series was subtracted.