frequency2scale?

[Kaare Mikkelsen]

Hello

I am seeing a factor-of-two error when I try to create a spectrogram using pywt.cwt. 

 

MWE:

 

timeVec=np.linspace(0,10,1000)

srate=1/(timeVec[1]-timeVec[0])

freqVec=np.linspace(1,10,1000)

signal=np.sin(2*np.pi*freqVec*timeVec)

 

freqs=np.linspace(0.1,20,30)

widths = np.sort(pywt.frequency2scale('cmor1.5-1.0',freqs/srate))

cwtmatr,cwtFreqs = pywt.cwt(signal, widths, 'cmor1.5-1.0',sampling_period=1/srate)

 

plt.imshow(np.abs(cwtmatr), extent=[0, 10, 1, 30], cmap='PRGn', aspect='auto',

           vmax=abs(cwtmatr).max(), vmin=-abs(cwtmatr).max())  

plt.yticks(ticks=np.arange(len(widths))+1,labels=np.round(freqs,1));

plt.colorbar()

plt.ylabel('Frequency (Hz)')

 

This shows the frequency rising linearly to 20 Hz, even though it is clearly only rising to 10 hz in the signal definition. Naturally, if I scale the frequencies by the Nyquist frequency instead of the sample rate, everything works out, but that seems to be a mistake?

 If have tried it with multiple different wavelets now.

Regards

Kaare Mikkelsen

 

 

[Deepu]

Hello Kaare Mikkelsen,

The factor-of-two error you are encountering in your code is related to the interpretation of frequencies in the wavelet transform. The wavelet transform is a time-scale analysis, where the scale parameter determines the frequency resolution.

In your code, when you calculate the widths using `widths = np.sort(pywt.frequency2scale('cmor1.5-1.0', freqs/srate))`, you are converting the frequencies to scales using the sample rate (`srate`). However, the scale parameter in the wavelet transform is not directly equivalent to the frequency. It represents the number of samples per cycle, which is inversely related to frequency.

To correct this, you need to use the Nyquist frequency (half the sample rate) instead of the sample rate when calculating the widths. Here's the modified code snippet:

```python
timeVec = np.linspace(0, 10, 1000)
srate = 1 / (timeVec[1] - timeVec[0])
freqVec = np.linspace(1, 10, 1000)
signal = np.sin(2 * np.pi * freqVec * timeVec)

freqs = np.linspace(0.1, 20, 30)
nyquist_freq = srate / 2
widths = np.sort(pywt.frequency2scale('cmor1.5-1.0', freqs / nyquist_freq))

cwtmatr, cwtFreqs = pywt.cwt(signal, widths, 'cmor1.5-1.0', sampling_period=1 / srate)

plt.imshow(np.abs(cwtmatr), extent=[0, 10, 1, 30], cmap='PRGn', aspect='auto',
           vmax=abs(cwtmatr).max(), vmin=-abs(cwtmatr).max())
plt.yticks(ticks=np.arange(len(widths)) + 1, labels=np.round(freqs, 1))

plt.colorbar()
plt.ylabel('Frequency (Hz)')

plt.show()
```

By using the Nyquist frequency (`nyquist_freq`) in the calculation of widths, the frequency range in the resulting spectrogram will be correctly aligned with the original signal frequency range.

I hope this helps! Let me know if you have any further questions.