calculating signal coherence with fft's

Share icon

Hello,

I am trying to calculate the spectral coherence of two signals using Igor's built in FFT functions. I am roughly following this example I found in a python package, trying to reproduce the final chart on the page: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.coher…

I was having trouble until I stumbled across Igor's DSPPeriodogram function and I can now calculate the signal coherence with the following code:

function coherence()
    variable n = 4096
    make/O/N=(n) signal_x
    setscale/P x, 0, 0.01, signal_x
    signal_x = gnoise(1)
    duplicate/O signal_x signal_y
   
    FilterFIR/LO={0,0.5,4096} signal_y
    signal_x += sin(2*pi*10*x)
    signal_y += gnoise(0.1)
   
    FFT/OUT=1/DEST=fft_x signal_x
    FFT/OUT=1/DEST=fft_y signal_y
    MatrixOp/O psd_xx = mag(fft_x*conj(fft_x))
    MatrixOp/O psd_yy = mag(fft_y*conj(fft_y))
    MatrixOp/O psd_xy = mag((fft_x)*conj(fft_y))
    MatrixOp/O coh_xy = mag(psd_xy/psd_xx/psd_yy)

    dspPeriodogram/COHR/SEGN={2048,0} signal_x, signal_y
    wave W_periodogram
    MatrixOp/O coh_xy_dsp = magsqr(W_periodogram)
end


Plotting W_periodogram gives a qualitatively similar plot to the one in the python documentation.

My question is: Why do I need to use the /SEGN flag to get a reasonable plot of the coherence. If I omit this flag, I find that the coherence is ~1.000 at almost every frequency. This is the same result that I find when I try to use the raw signal FFT's to estimate the power spectra and then the coherence:

function fft_coherence()
    wave signal_x, signal_y
   
    FFT/OUT=1/DEST=fft_x signal_x
    FFT/OUT=1/DEST=fft_y signal_y
    MatrixOp/O psd_xx = mag(fft_x*conj(fft_x))
    MatrixOp/O psd_yy = mag(fft_y*conj(fft_y))
    MatrixOp/O psd_xy = magsqr((fft_x)*conj(fft_y))
    MatrixOp/O coh_xy = psd_xy/psd_xx/psd_yy  // Following the convention in the DSPPeriodogram documentation??

end


Is there some fundamental reason you cannot use the fft magnitude squared to estimate power spectra between two signals for coherence? I see that in the python example they are using the welch psd estimate, which is averaging the power spectrum in a way similar to the /SEGN flag. I want to be able to implement the coherence using FFT's because I want to calculate coherence for a 2D wave.

I guess this is more of a signal processing question than an Igor question, so thank you for your help.
Take a look at the equations for coherence and you will see that they involve some averages. Without the /SEGN you are basically evaluating a ratio where the numerator and denominator are effectively the same. With /SEGN you are introducing some averaging, though it looks as insufficient for your example.

April 2, 2017 at 11:04 am - Permalink

Thanks for your reply,

From what I gather, the magnitude squared of the FFT is a poor estimate of the power spectral density of a discretely sampled signal. Using this method, the SNR of the power is ~1, very poor. To improve the estimate of the PSD, you can divide up the signal into shorter segments (e.g. by using the /SEGN flag) and averaging the magnitude squared FFTs together. Thus, you loose some resolution in frequency, but gain accuracy in the estimation of the power.

April 12, 2017 at 09:48 pm - Permalink