[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

*Subject*: Re: cross-spectral analysis*From*: Craig Markwardt <craigmnet(at)cow.physics.wisc.edu>*Date*: 09 Jul 1999 19:24:53 -0500*cc*: t.osborn(at)uea.ac.uk*Newsgroups*: comp.lang.idl-pvwave*Organization*: U. Wisc. Madison Physics -- Compact Objects*References*: <7m202t$ld8$1@cpca14.uea.ac.uk>*Reply-To*: craigmnet(at)cow.physics.wisc.edu*Xref*: news.doit.wisc.edu comp.lang.idl-pvwave:15614

t.osborn@uea.ac.uk (T.Osborn) writes: > > > I need to do some cross-spectral analysis (preferably in IDL). > In particular, I need to compute the "squared coherence" > between two time series. Is there a built-in routine to do > this in IDL, or does someone have their own? I'm not very > au fait with spectral analysis, so I hope it's not too > difficult! Yes, this is definitely do-able, I have a specialized routine to do it for light curves in X-ray astronomy. There are unfortunately many definitions for the squared coherence in the engineering literature. Sometimes it can be difficult to understand which one to use. Generally speaking, if F1 and F2 are the complex fourier coefficients of two time series at a given frequency, then the squared coherence is gamma^2 = |<F1* F2>|^2 / ( <|F1|^2> <|F2|^2> ) where |x| is the absolute value of x, <x> is the ensemble average of x, and x* is the complex conjugate of X. It is very important to realize that the coherence cannot be computed from a single measurement. You need to measure the fourier coefficients many times, and average over that entire sample. In X-ray astronomy, this often means breaking the observation up into small segments and doing lots of FFTs. The classic reference for this kind of thing is Bendat and Piersol, "Random Data: Analyses and Measurement Procedures." How to do it in IDL? Well, it's probably not appropriate to regurgitate all of my code. Here is a snippet of the important bits. for i = 0, nmeasurements - 1 do begin ;; You must construct the time series here ts1 = [Time series 1, measurement i] ts2 = [Time series 2, measurement i] ;; Apply your window function as needed. ;; Compute the FFT of each time series f1 = fft(ts1, +1) f2 = fft(ts2, +1) ;; f1 and f2 are complex vectors ;; Extract only the useful elements of FFT output nfbins = n_elements(ts1) / 2 f1 = f1(1:nfbins) ;; Ignores DC value, but keeps Nyquist value f2 = f2(1:nfbins) ;; Compute complex cross amplitude c12 = CONJ(f1) * f2 ;; Cross amplitude 1 vs. 2 p1 = abs(f1)^2 ;; Power spectrum series 1 p2 = abs(f2)^2 ;; Power spectrum series 2 ;; Accumulate many samples if i EQ 0 then begin sum_c12 = c12 sum_p1 = p1 sum_p2 = p2 nsamples= 1L endif sum_c12 = sum_c12 + c12 sum_p1 = sum_p1 + p1 sum_p2 = sum_p2 + p2 nsamples= nsamples+ 1 endelse endfor ;; Compute squared coherence value gam2 = abs(sum_c12)^2 / (sum_p1 * sum_p2) ;; Uncertainty in the mean of squared coherence value egam2 = sqrt( (2D * gam2*(1D - gam2)^2) / nsamples ) The result, gam2, is the squared coherence value at each frequency, and egam2 is the standard error of the mean of gam2. Of course, there are tweaks. Several corrections can apply if your signal is noise-dominated. Computing the lag between two signals is also straightforward: ;; Lag at each freq., computed in radians lag = atan(imaginary(sum_c12), double(sum_c12)) ;; Uncertainty in the lag elag = sqrt( (1D - gam2)/ (2D * gam2) / nsamples ) The resulting value, lag, is positive when the first series leads the second. Good luck! Craig -- -------------------------------------------------------------------------- Craig B. Markwardt, Ph.D. EMAIL: craigmnet@astrog.physics.wisc.edu Astrophysics, IDL, Finance, Derivatives | Remove "net" for better response --------------------------------------------------------------------------

- Prev by Date:
**Array index bug....** - Next by Date:
**Re: Array index bug....** - Prev by thread:
**Re: help with relational operators** - Next by thread:
**Too many users...** - Index(es):