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
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?
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):
With linear frequency axes (shows general coherency):
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.
