Continuous Wavelet Transform



In [1]:

    
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline



In [2]:

    
tremor = np.loadtxt('benchmark_signals/tremor.txt')



In [3]:

    
synthetic = np.loadtxt('benchmark_signals/synthetic.txt')



In [4]:

    
speech = np.loadtxt('benchmark_signals/speech.txt')



In [244]:

    
def tf(signal, fs, w=256, wtime=False, poverlap=None, ylim=None, colorbar=False, vmin=None, vmax=None):
    dt = 1./fs
    n = signal.size
    t = np.arange(0.0, n*dt, dt)
        
    if wtime:
        # Then the window length is time so change to samples
        w *= fs
        w = int(w)
        
    if poverlap:
        # Then overlap is a percentage
        noverlap = int(w * poverlap/100.)
    else:
        noverlap = w - 1
        
    plt.figure(figsize=(12,8))
    plt.subplot(211)
    plt.plot(t, signal)
    ax1 = plt.subplot(212)
    Pxx, freqs, bins, im = plt.specgram(signal, NFFT=w, Fs=fs, noverlap=noverlap, cmap='Greys', vmin=vmin, vmax=vmax)
    if colorbar: plt.colorbar()
    if ylim:
        ax1.set_ylim((0,ylim))
    plt.show()

`scipy.signal.cwt`



In [6]:

    
from scipy import signal



In [7]:

    
sig = np.random.rand(20) - 0.5
wavelet = signal.ricker
widths = np.arange(1, 11)
cwtmatr = signal.cwt(sig, wavelet, widths)



In [8]:

    
plt.plot(cwtmatr)
plt.show()



In [9]:

    
widths = np.arange(1, 11)
cwtmatr = signal.cwt(synthetic, wavelet, widths)



In [10]:

    
plt.plot(cwtmatr)
plt.show()

I assume I at looking at responses at different scales... Not sure how to plot this. Let's try something else.

Roger Fearick, University of Cape Town

From his Pyscellania page.



In [11]:

    
import pycwt



In [237]:

    
#cwt = pycwt.Paul(synthetic, largestscale=40, notes=24, order=12, scaling="log") # 'order' specifies order :)
#cwt = pycwt.DOG4(synthetic, largestscale=200, notes=36, scaling="log")
#cwt = pycwt.Haar(synthetic, largestscale=40, notes=24, scaling="log") # Has type bug
#cwt = pycwt.MexicanHat(synthetic, largestscale=200, notes=36, scaling="log")
#cwt = pycwt.MorletReal(synthetic, largestscale=100, notes=24, scaling="log")

# This is the only one that looks decent... 
cwt = pycwt.Morlet(synthetic, largestscale=100, notes=12, scaling="log")



In [241]:

    
plt.figure(figsize=(15,5))
plt.imshow(cwt.getpower(), aspect="auto", cmap="Greys", vmax=20)
plt.colorbar()
plt.show()



In [250]:

    
tf(synthetic, 800, w=128, ylim=256, vmin=-30)

There are several other wavelets:

Morlet
MorletReal
MexicanHat
Paul2 : Paul order 2
Paul4 : Paul order 4
DOG1 : 1st Derivative Of Gaussian
DOG4 : 4th Derivative Of Gaussian
Haar : Unnormalised version of continuous Haar transform
HaarW : Normalised Haar

Let's wrap everything:



In [65]:

    
def cwt(data):
    w = pycwt.Morlet(data, largestscale=12, notes=24, order=2, scaling="log")
    
    scales = w.getscales()
    pwr = w.getpower()
    
    scalespec = np.sum(pwr, axis=1) / scales
    y = w.fourierwl * scales
    x = np.arange(0, data.size, 1)
    
    plt.figure(figsize=(15,5))
    plt.imshow(pwr, aspect="auto", cmap='Greys')
    plt.colorbar()
    plt.show()



In [66]:

    
cwt(speech)



In [25]:

    
tf(speech, 8000)



In [ ]: