Why does using a log scale ruin (warp) my scaleograms?

Question

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)

Results with log scaling (incorrect)

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 = (2abs(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')

Answer 1

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,