Freq-to-scale, scale-to-freq

[scottfleming]

My understanding is that the ssq_cwt function optionally accepts a set of wavelet scales rather than corresponding frequencies. While I understand that there are a few different interpretations of "corresponding frequencies" for wavelet scales, I've noted that there is a simple, analytical relationship (see Equation A3 of "An Introduction to Wavelet Analysis in Oceanography and Meteorology", also referenced in "A Practical Guide to Wavelet Analysis") between wavelet scale and frequency for certain wavelet functions. Would it be difficult to add functionality for the user to pass in frequences rather than wavelet scales and have those converted to wavelet scales on the backend?

[OverLordGoldDragon]

Good point. Indeed the conversion will be simple via a two-way center_frequency method. I'll add this to the todo list but I don't intend to work on this repository anytime soon. For now I recommend something like:

Freq-to-scale

code

import numpy as np
from ssqueezepy import Wavelet, center_frequency
from ssqueezepy.utils import cwt_scalebounds
from ssqueezepy.visuals import plot


def logb(x, base=2):
    return np.log(x) / np.log(base)

def freqs_to_scales(freqs, kind='peak', base=2):
    def log(x):
        return logb(x, base)
        
    assert np.all(freqs >= 0), "frequencies must be positive"
    assert freqs.max() <= np.pi,     "max frequency must be <=pi"
    assert freqs.max() == freqs[-1], "max frequency must be last sample"
    assert freqs.min() == freqs[0],  "min frequency must be first sample"

    smin, smax = cwt_scalebounds(wavelet, N, preset='maximal')
    search_scales = np.logspace(log(smin), log(smax), M, base=base)
    
    f_from_scales = []
    for scale in search_scales:
        f = center_frequency(wavelet, scale, N, kind=kind)
        f_from_scales.append(min(max(f, 0), np.pi))
        
    # pick closest match
    fmin, fmax = freqs.min(), freqs.max()
    smax = search_scales[np.argmin(np.abs(f_from_scales - fmin))]
    smin = search_scales[np.argmin(np.abs(f_from_scales - fmax))]
    # make scales between found min and max
    scales = np.logspace(log(smax), log(smin), M, base=base)
    
    return scales

#%%###########################################################################
N = 1024
wavelet = Wavelet(('morlet', {'mu': 14}))

M, base = 100, 2
fmin, fmax = .01, np.pi/2
freqs = np.logspace(logb(fmin, base), logb(fmax, base), M, base=base)
scales = freqs_to_scales(freqs, kind='peak', base=base)

#%%
plot(freqs,  show=1, title="Frequencies")
plot(scales, show=1, title="Scales")

This assumes freqs is logspaced, which can be adjusted for linear spacing, and it is approximate.

[OverLordGoldDragon]

Scale-to-freq 1.0

Guaranteed to work*, but only kind='peak' supported (* though not tested exhaustively)

code

import numpy as np
import warnings
from ssqueezepy import Wavelet

def scale_to_freq(scales, wavelet, N, fs=None):
    # process args
    if isinstance(scales, float):
        scales = np.array([scales])    
    wavelet = Wavelet._init_if_not_isinstance(wavelet)
    
    # evaluate wavelet at `scales`
    Npad = int(2**np.ceil(np.log2(N)))
    psis = wavelet(scale=scales, N=Npad)
    # find peak indices
    idxs = np.argmax(psis, axis=-1)
    
    # check
    # https://github.com/OverLordGoldDragon/ssqueezepy/issues/41
    if np.any(idxs > Npad//2) or 0 in idxs:
        warnings.warn("found potentially ill-behaved wavelets (peak indices at "
                      "negative freqs or at dc); will round idxs to 1 or N/2")
        n_psis = len(psis)
        for i, ix in enumerate(idxs):
            if ix > Npad//2 or ix == 0:
                if i > n_psis // 2:  # low freq
                    idxs[i] = 1
                else:  # high freq
                    idxs[i] = Npad//2
    # convert
    freqs = idxs / Npad  # [0, ..., .5]
    assert freqs.min() >= 0,   freqs.min()
    assert freqs.max() <= 0.5, freqs.max()
    
    if fs is not None:
        freqs *= fs   # [0, ..., fs/2]
    return freqs

scales = np.array([1, 100, 1000])
wavelet, fs = 'gmw', 400
freqs = scale_to_freq(scales, wavelet, N=2048, fs=fs)

print("wavelet, fs: %s, %s" % (wavelet, fs))
print("scales: %s" % scales)
print("freqs:  %s" % freqs)

wavelet, fs: gmw, 400
scales: [   1  100 1000]
freqs:  [172.8515625   1.7578125   0.1953125]

Scale-to-freq 2.0

Might work for any kind

code

import numpy as np
from ssqueezepy import Wavelet, center_frequency
from ssqueezepy.utils import cwt_scalebounds


def scale_to_freq(scales, wavelet, N, fs=None, kind='peak'):
    # process args
    if isinstance(scales, float):
        scales = np.array([scales])    
    wavelet = Wavelet._init_if_not_isinstance(wavelet)
    # find bounds
    smin, smax = cwt_scalebounds(wavelet, N, preset='maximal')
    
    # ensure user input is within bounds
    # https://github.com/OverLordGoldDragon/ssqueezepy/issues/41
    if scales.min() < smin:
        raise ValueError("min `scales` ({:.1e}) exceeds bound ({:.1e}); "
                         "wavelet here is ill-defined.".format(scales.min(), smin))
    if scales.max() > smax:
        raise ValueError("max `scales` ({:.1e}) exceeds bound ({:.1e}); "
                         "wavelet here is ill-defined.".format(scales.max(), smax))
    
    # convert
    freqs = np.array([center_frequency(wavelet, s, N, kind=kind)
                      for s in scales])
    # radian to linear
    freqs /= (2*np.pi)
    assert freqs.min() >= 0,   freqs.min()
    assert freqs.max() <= 0.5, freqs.max()
    # physical units (else result is 0 to 0.5)
    if fs is not None:
        freqs *= fs
    return freqs

scales = np.array([1, 100, 1000])
wavelet, fs = 'gmw', 400
freqs = scale_to_freq(scales, wavelet, N=2048, fs=fs)

print("wavelet, fs: %s, %s" % (wavelet, fs))
print("scales: %s" % scales)
print("freqs:  %s" % freqs)

wavelet, fs: gmw, 400
scales: [   1  100 1000]
freqs:  [172.8515625   1.7578125   0.1953125]

[ghost]

Just question, I'm plotting y-axis with my function below,
is that conversion correct?

ssqueezepy_scale2freq

[OverLordGoldDragon]

@fairlight-kr Please open a separate issue. And code should be standalone (I don't have your data.txt, use randn). But no, freq != 1 / scale - see article.