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

```python 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

```python 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

```python 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.