The goal

I wish to compute the coherence estimate using the continuous wavelet transform (CWT) of a real-valued signal with the complex Morlet (a.k.a. Gabor) wavelet. I compute the cwt of the signal in a very similar way as in the pycwt library, but I use the pyfftw library for faster FFT computation. I'm confident that the cwt is computed correctly, because I get nearly identical results to pycwt.cwt and the scalograms compare well to STFT spectrograms. Although I use Python, I think it is not language specific.

Coherence calculation resources

However, I'm not so confident as to how to compute the coherence estimate. I based my coherence calculation on the classic Torrence & Webster 1999 paper (as does pycwt) where prior to computing the coherence fraction they smooth the cross-power density (cwt(sig)1 * cwt(sig2).conj()) and power densities (abs(cwt(sig))**2) in time with a Gaussian kernel with a standard deviation equal to the scale for each scale. This makes sense to me, because the spectral characteristic changes at a slower rate for larger scales, so a larger smoothing window is necessary than for smaller scales. They also smooth them in the scale domain with a simple moving average, which makes less sense to me.

Issues with the definition and implementation

Here are my issues with the implementation and/or definition of this smoothing operator for which I seek explanation and/or help

  1. Smoothing in the time domain

    • The Gaussian smoothing in time doesn't seem to have a wide-enough smoothing window at small scales (I suspect because the std. deviation is close to the sampling step so it doesn't smooth much) and the resulting coherence at small scales (high frequencies close to Nyquist) is almost always unity. I compute this convolution via multiplication with FFT (as does pycwt). Is this Gaussian window inappropriate for some reason in my case or is there something wrong with my implementation?
    • Smoothing with a moving average produces much more reasonable results at high frequencies, but understandably worse at low frequencies (large scales). It seems the MATLAB wavelet toolbox also uses a moving window for time domain smoothing in wcoher. Is this the recommended approach then? What about lower frequencies?
  2. Smoothing in the scale domain

    • The moving average smoothing in the scale domain seems to be empirical, but I can live with that. I's just that there seems to be now real explanation why it is necessary (it isn't done in Welch's method of coherence estimation as far as I know). Is there a good reason to perform smoothing in the scale domain?
    • What really confuses me is that pycwt computes it on logarithmically-spaced scales, which is a non-linear filter and I'm not sure whether that makes any sense. Does it make sense even as a non-linear filter?

Results so far

So far I've used just a simple moving average filter only in the time domain and that gives reasonable results above low frequencies. However, I don't feel confident about the results and I will need coherence at lower frequencies in the future, so I'm asking here. I also tried using the Gaussian window with a much larger standard deviation to get results similar to the moving average window smoothing results.

The thick lines in every first subplot are lowpass-filtered versions of the original, semi-transparent signals. I added these lowpassed signals to make it simpler to distinguish the two signals.

With logarithmic frequency axes (shows low frequencies better): example of xwt with log freq axes With linear frequency axes (shows general coherency): enter image description here This is the part of the code that I use for smoothing, fftmod is pyfftw.interfaces.numpy_fft (which provides the same functions as the numpy.fft module)

# smoothing in time domain
freqs = fftmod.fftfreq(self.time_axis.data.shape[0],
                       self.axis_step(0))[:,np.newaxis]
freqs0 = self.axis1.data
if smoothing_window:  # moving average
    if isinstance(smoothing_window, complex):
        # in time units, convert to number of points
        smoothing_window = abs(smoothing_window) / self.axis_step(0)
    smoothing_window = np.empty(smoothing_window)[:,np.newaxis]
    smoothing_window.fill(1.0 / smoothing_window.size)
    smoothing_ft = fftmod.fft(smoothing_window, n=freqs.size,
                              axis=0, **fftw_kwargs)
else: # gaussian smoothing, gaussian_std_mul multiplies the std dev
    # FT of smoothing window which is a gaussian with std dev 1/f0
    # is a gaussian with std dev f0 without f0 in normalization const
    smoothing_ft = np.exp(-(freqs/freqs0 * gaussian_std_mul)**2 / 2.0) / _SQRT_2PI
# smooth by convolution with window as multiplication in Fourier domain
csd = fftmod.ifft(smoothing_ft * fftmod.fft(
    csd, axis=0, **fftw_kwargs), axis=0, **fftw_kwargs) # TODO maybe fft of vstack
psd_self = fftmod.ifft(smoothing_ft * fftmod.fft(
    psd_self, axis=0, **fftw_kwargs), axis=0, **fftw_kwargs).real
psd_signal = fftmod.ifft(smoothing_ft * fftmod.fft(
    psd_signal, axis=0, **fftw_kwargs), axis=0, **fftw_kwargs).real

More recent resources

I did some further digging and found an article from 2007 by Bernard Cazelles. The code available on his website (includes the article in the zip file) seems to suggest that he also uses smoothing both in the scale and frequency domain with various different windows (Barttlet, Hamming, Boxcar,...), but does use different window lengths for different scales in the time domain and just one length in the scale domain. The code shows that the scale domain is again logarithmic and yet is smoothed with a linear window. However, at least he cites Chatfield 1989 who cites a suggestion from P. J. Daniell in 1946 that smoothing in the frequency domain is also needed.

I'm not sure if I can help you with a perfect solution, but hopefully at least with some hints.

As far as I can tell the kernel for the smoothing filter of the wavelet (power/amplitude) spectrum has to scale with the size of the reproducing kernel of the wavelet, i.e. the window width in time scales with wavelet scale, and window width in scale depends on your 'sampling' in scale direction – for instance you use a fixed window width for log scaling, but you would need to adapt the smoothing filter window width if you would decide for a linear scaling in scale, which is also quite common in cwt analysis.

Regarding problems close to nyquist I can say that any analysis of frequencies close to nyquist have to be taken with care as your approximation of the continuous (!) wavelet transform get's pretty bad close to nyquist. Further you smoothing filter will cause boundary effects especially close to nyquist as you cannot go above.

Finally I want to recommend you this PhD Thesis by Douglas Maraun which focuses exactly on this topic and extends the work done by Torrence and formalizes it in a mathematical more rigorous way. There is some software package available as well, not sure, however, if this is still maintained.

  • Thank you for your answer, it helped a lot. That thesis looks very well. However, I'm still not entirely sure why on log-scaled scales a "fixed window" is used. Specifically, I'm not sure what is meant by "fixed". Could you please expand on that? – Ondřej Grover Jun 19 '15 at 4:45
  • With fixed I mean a constant window width along the scale axis. As mentioned, the smoothing window width should scale in accordance with the width of the wavelet's reproducing kernel in scale and time, respectively. The reproducing kernel is the wavelet transform of the wavelet itself. As you move one wavelet along the scale axis and look at its width in scale you will notice that this width stays constant if a log scale is chosen. This reflects simply the 'natural' scaling property of the wavelet transform. If you sample differently along the scale axis the shape gets 'distorted'. – André Bergner Jun 21 '15 at 18:42

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.