-1
$\begingroup$

I'm trying to make a scaleogram of neural data using the continuous wavelet transform (with Morlet wavelets). As a starting point, I use just a simple sine wave built from two frequencies, 50 and 80 Hz.

Everything looks great when using a linear scale, but since I expect most of the interesting results will be at low frequencies, I'd like to use a log scale. For some reason, once I do this the results look totally wrong -- either stretched or shifted on the graph, and definitely not corresponding to the expected y-axis values.

Does anyone know why this might be?

(Notice the graph in the bottom right)

Results with linear scaling (correct)

enter image description here

Results with log scaling (incorrect)

enter image description here

Code for dummy signal + parameters

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

from scipy import signal
from scipy import fftpack
from scipy.fft import fftshift

fs = 1000.
f1 = 80.
f2 = 50.
t = np.arange(5000)
signal1 = np.sin(2. * np.pi * f1 * t / fs) + np.sin(2. * np.pi * f2 * t / fs)

# Wavelet parameters 
min_freq = 1
max_freq = int(fs / 4)
num_freqs = int(max_freq / 2)

Scaleogram code (uncomment the 1-3 marked lines to change to linear or log scaling) (adapted from here)

# parameters
srate = 1000
time = t
pnts  = len(time)

### static power spectrum
hz = np.linspace(0,srate/2,int(pnts/2)+1)
powr = (2*abs(np.fft.fft(signal1)/pnts))**2

### time-frequency analysis 
nfrex = num_freqs 
##### FOR LOG SCALING UNCOMMENT
# frex = np.logspace(np.log10(min_freq), np.log10(max_freq), num_freqs) 
##### FOR LINEAR SCALING UNCOMMENT
# frex = np.linspace(min_freq, max_freq, num_freqs)
fwhm = np.linspace(.8,.7,nfrex)

# setup wavelet and convolution parameters
wavet = np.arange(-2,2,1/srate)
halfw = int(len(wavet)/2)
nConv = pnts + len(wavet) - 1

# initialize time-frequency matrix
tf = np.zeros((len(frex),pnts))

# spectrum of data
dataX = np.fft.fft(signal1,nConv)

# loop over frequencies
for fi,f in enumerate(frex):
    
    # create wavelet
    waveX = np.fft.fft( np.exp(2*1j*np.pi*f*wavet) * np.exp(-4*np.log(2)*wavet**2/fwhm[fi]**2),nConv )
    waveX = waveX / np.abs(max(waveX)) # normalize
    
    # convolve
    ast = np.fft.ifft( waveX*dataX )
    # trim and reshape
    ast = ast[halfw-1:-halfw]
    
    # power
    tf[fi,:] = np.abs(ast)**2
    
### plotting
fig = plt.figure(figsize=(8,7),tight_layout=True)
gs = gridspec.GridSpec(2,2)

ax = fig.add_subplot(gs[0,:])
ax.plot(time,signal1,'b',linewidth=.5)
ax.set_xlabel('Time (s)')
ax.set_ylabel('Amplitude')
ax.set_title('Time domain')

ax = fig.add_subplot(gs[1,0])
ax.plot(hz,powr[:len(hz)],'b',linewidth=1)
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('Power')
ax.set_xlim([0,100])
ax.set_title('Frequency domain')

##### FOR LOG SCALING UNCOMMENT
# rounded = np.logspace(np.log10(min_freq), np.log10(max_freq), 10).round(1)
ax = fig.add_subplot(gs[1,1])
ax.imshow(np.flip(tf),origin='upper',aspect='auto',extent=[time[0],time[-1],frex[0],frex[-1]],cmap='Spectral')
ax.set_xlabel('Time (sec.)')
ax.set_ylabel('Frequency (Hz)')
##### FOR LOG SCALING UNCOMMENT
# ax.set_yscale('log')
# ax.set_yticks(ticks=rounded, labels=rounded)
ax.set_title('Time-frequency domain')
$\endgroup$

1 Answer 1

0
$\begingroup$
filters = []
for fi,f in enumerate(frex):
    ...
    filters.append(waveX)
filters = np.abs(np.array(filters))
plt.plot(filters.T)

Your filterbank is very poorly designed, including the linear one (despite looking okay for this example). You're also not using wavelets since filters aren't zero-mean, and not doing wavelet transform since filter bandwidth is independent of f. Refer to this post.

The reason you're observing one peak in the log case is a direct consequence; there are, in fact, two peaks, and a now-deleted answer (whose undeletion I encourage) has shown that undoing the log scaling as a post-processing step will help see it. But it's not necessary; just plot a time slice:

plt.plot(tf[:, 0])

It's attenuated because the filterbank incompletely tiles the frequency domain. It's by sheer chance anything was captured at all; try e.g. f1, f2 = 150, 100. Regarding the y-axis values, simply remove np.flip(tf), or reverse rounded; also for better formatting,

ax.imshow(np.flip(tf),origin='upper',aspect='auto',cmap='Spectral')
idxs = np.linspace(0, len(rounded) - 1, 8).astype('int32')
yt = ["%.2f" % h for h in np.asarray(rounded)[idxs]][::-1]
ax.set_yticks(idxs)
ax.set_yticklabels(yt)

and using f1, f2 = 75, 50 so we're not unlucky,

$\endgroup$
3
  • $\begingroup$ Unless I misunderstood something, this really doesn't address my plotting issue. It also seems pretty unproductive to just keep replying to peoples issues saying that SciPy and PyWavelets are too bad to do this task well, while in the comments you link to you just reference yourself. $\endgroup$ Mar 22, 2022 at 20:47
  • $\begingroup$ @neverreally If you're not computing a scalogram, it won't look like a scalogram, is what I'm saying. My link is meant to be an aid to debugging your filterbank design; SP and PYWT just happen to suffer from the problems described in that post. And yes, they are bad and should be called out such, though one is free to keep using them if they please. Regardless I've added some detail. $\endgroup$ Mar 22, 2022 at 21:02
  • $\begingroup$ @neverreally I'll admit I missed the part on y-axis values in light of other issues; now added. $\endgroup$ Mar 22, 2022 at 21:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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