ssqueezy notes

[OverLordGoldDragon]

Thread for various implementation or theoretical notes / ideas.

[OverLordGoldDragon]

Shifting center frequency shifts wavelet (Morlet), but doesn't change its width (of underlying continuous-time function that we sample). Higher center frequency thus suffers worse max scale behavior, as second bin, the smallest possible for symmetric psih, is at peak, but bins 1 and 3 are near zero, thus rendering psi a near-pure sinusoid, with awful localization. (Localization is awful either way, but in pure-sine case std_t=inf; asymmetric but greater-valued side-bins are preferred over perfect symmetry but low side-values, as lower scales show)

N, mu, scale = 1024, 5, 407.44
wavelet = Wavelet(('morlet', {'mu': mu}))
w1 = np.linspace(0, 20, 4096)
w2 = np.linspace(0, np.pi, N//2) * scale

vline = (mu, {'color': 'tab:red', 'linestyle': '--'})
plot(w1, wavelet(w1), vlines=vline)
scat(w2, wavelet(w2), color='r', show=1)

plot(wavelet(w2)[:10])
scat(wavelet(w2)[:10], show=1)

psi = wavelet.psifn(scale=scale, N=N)
plot(psi, complex=1)
plot(np.abs(psi), color='k', linestyle='--', show=1)

N, mu, scale = 1024, 5*2, 407.44*2

---

or... two-point symmetry about a "virtual bin" is actually more minimal. At which bin such symmetry becomes possible increases with increased center frequency, but is still more permissive than three-point symmetry.

N, mu, scale = 1024, 5, 543.25

N, mu, scale = 1024, 5*2, 325.95*2

[OverLordGoldDragon]

Greater center frequency leaves fewer "well-behaved" time-domain scales (or rather shifts the range of well-behaved-ness) and makes edge scale behavior messier; larger N helps with low scales but much less so with large scales. Note suboptimal behavior in time-domain is directly parred with suboptimal in freq-domain.

(animation idea: sweep scales and show above three plus dot on below to show current scale)

[OverLordGoldDragon]

1. Logscale CWT toward max possible scale is very wasteful:

can remedy by downsampling high scales, and maybe sampling in linspace:

but introduces the burden of treating parts of CWT separately, e.g. in synchrosqueezing and reconstruction.

2. The redundancy seems worse with more frequency-localized (lesser freq-domain width of CT-function) wavelets, which operate on effectively a single point in frequency domain, thus rescaling the exact same time-domain wavelet, as opposed to a sum of two or more as with more time-localized:

3. cwt_scalebounds automatic min and max scale auto-setting can be improved; currently maximal manages to evade k=1 bin, which is undesired.
4. kymatio's scales sampling appears to be optimally efficient, but as-is does it in a NSGT manner, i.e. both STFT & CWT, also does poorly near Nyquist or k=1, but all may be easy to remedy once understood.

---

Code

import numpy as np
from ssqueezepy.toolkit import cos_f
from ssqueezepy.visuals import plot, plotscat, imshow
from ssqueezepy import cwt, Wavelet

#%%###########################################################################
def fn(x, wavelet, scales, color='tab:blue'):
    Wx, scales, *_ = cwt(x, wavelet, scales=scales)
    imshow(Wx, abs=1, cmap='jet')
    Psih = wavelet(scale=scales)
    plot(Psih[-90//3:, :13].T, title="last 30 wavelets", color=color, show=1)
    return scales

#%%###########################################################################
x  = cos_f([1], 513, endpoint=1)
x += cos_f([8], 513, endpoint=1)
wavelet = Wavelet(('gmw', {'beta': 60}))

#%%###########################################################################
scales = fn(x, wavelet, 'log:maximal')

#%%###########################################################################
a, b = 1, 100
d = b - a
scalesl = np.linspace(scales[-b], scales[-a], d//3, endpoint=1)
scales2 = np.hstack([scales[:-b][:,0], scalesl[:,0]])
scales3 = np.hstack([scales[:-b][:,0], scales[-b::3][:,0]])

plot(scales, xlims=(-5, len(scales)*1.08))
plotscat(scales2, color='tab:orange')
plotscat(scales3, color='tab:green', show=1)

#%%
fn(x, wavelet, scales=scales2, color='tab:orange')
fn(x, wavelet, scales=scales3, color='tab:green')

[OverLordGoldDragon]

Re: prev comment. Problem: waste at high scales. Solution: downsample high scales. Problem: inversion & synchrosqueezing expect constant sampling rate of scales. Solution: do parts separately. Problem: where to separate? Solution: at red dot:

Code

import numpy as np
from ssqueezepy.visuals import plot, scat

scales = np.exp(np.linspace(0, 6, 256))
b = 120  # number of scales to downsample
d = 3    # downsample factor
scalesd = np.hstack([scales[:-b], scales[-b::d]])

# index of where downsampling begins;
# second finite difference nonzero only at one point
idx = np.argmax(np.diff(np.diff(np.log(scalesd)))) + 1

plot(scalesd)
scat(idx, scalesd[idx], color='tab:red', s=40)

This has been implemented as scales='log-piecewise'. The algorithm's logic is also lot more reliable than that currently behind cwt_scalebounds; it should probably used for min and max scale determination.

[OverLordGoldDragon]

Low frequency wavelet design

Seems best accomplished by sampling freq-domain wavelets in an upsampled length, then unpadding the convolution more (or trimming the wavelet for storage).

Below is a Morlet sampled at length 16384, and its time-domain waveform (centered):

then sampled at length 65536

and trimmed to same length

This has same frequency content but with better time decay. The spectrum of this trimmed segment is:

It has bits of negative and complex bin values, along few negative frequencies, and even DC; there's no way we'd sample this directly from a continuous-time function.

If one doesn't mind boundary effects, this presents a superior design alternative: sample longer then trim the wavelet. While the trimmed wavelet is entirely legitimate in the original (untrimmed) frame, it contains a slight dc after trimming which should be adjusted for somehow.

(off-topic upload: jtfs_fbank.zip)

D.zip