# Chapter 12. Morlet Wavelets and Wavelet Convolution


```python
import numpy as np
import scipy.io
from matplotlib import pyplot as plt

from numpy import pi, sin, cos, exp, sqrt, log, random  #import basic functions from numpy that we'll need
```


```python
data=scipy.io.loadmat('../data/sampleEEGdata')
```


```python
EEGdata=data["EEG"][0,0]["data"]
srate = float(data["EEG"][0,0]["srate"][0,0])
```

### Figure 12.1
The real component of a Morlet wavelet is a cosine multiplied by a Gaussian window:
\begin{equation}
Re(\text{Morlet wavelet}) = \cos(2\pi f t) e^{-t^2/(2s^2)}
\end{equation}


```python
time = np.arange(-1,1+1/srate,1/srate)

f=4 #sinewave frequency, Hz

#create a sinewave (cosine wave)
sine_wave = cos(2*pi*f*time)

#create a Gaussian

s=4/(2*pi*f) #standard deviation

gaussian_win = exp(-time**2/(2*s**2))

#plot our first wavelet!
_=plt.plot(time,sine_wave*gaussian_win)
```

### Figure 12.2


```python
fig, ax = plt.subplots(5,1)
ax[0].plot(np.squeeze(EEGdata[46,:,0]))

sine_wave = cos(2*pi*12*time) # 12Hz cosine wave
ax[1].plot(time, sine_wave)

#boxcar envelope
boxcar = np.zeros(len(sine_wave))
midpoint = (len(time))//2
boxcar[midpoint-int(srate//12//5):midpoint+int(srate//12//2)] = 1
ax[2].plot(time, sine_wave*boxcar)

#boxcar of different length
boxcar = np.zeros(len(sine_wave))
midpoint = int(len(time))//2
boxcar[midpoint-50:midpoint+50] = 1
ax[3].plot(time, sine_wave*boxcar)

s = 1.5/(2*pi*f)
gaussian_win = exp(-time**2/(2*s**2))
ax[4].plot(time, sine_wave*gaussian_win)

fig.tight_layout()
```

### Figure 12.3


```python
srate = 500. #sample rate in Hz
f = 10
time = np.arange(-1,1,1/srate)

#complex sinusoid
sine_wave = exp(2*pi*1j*f*time)

#Gaussian window
s = 6/(2*pi*f)
gaussian_win = exp(-time**2/(2*s**2))

#together they make a complex morlet wavelet!
wavelet = sine_wave*gaussian_win

#create plots for each component
fig, ax = plt.subplots(3,1,sharex='all', sharey='all')

ax[0].plot(time,np.real(sine_wave))
ax[0].set(title="sine wave")

ax[1].plot(time,gaussian_win)
ax[1].set(title="gaussian window")


ax[2].plot(time, np.real(wavelet))
ax[2].set(title="my first wavelet", xlabel='time')

fig.tight_layout()
```

### Figure 12.4


```python
num_wavelets = 80 # number of frequency bands
lowest_frequency = 2 #in Hz
highest_frequency = 100 # in Hz

#(linear) equally spaced frequencies for our wavelet family
frequencies = np.linspace(lowest_frequency, highest_frequency, num_wavelets)

plt.figure()
plt.plot(frequencies)
plt.xlabel("Frequency order")
_=plt.ylabel("Frequency in Hz")
```


```python
#initialize our wavelet family
wavelet_family = np.zeros([num_wavelets,len(time)])*1j #1j is to create a complex array of zeros

#iterate through freqs and make a wavelet family

for fi in range(num_wavelets):
    
    #create a sine wave
    sinewave = exp(2*1j*pi*frequencies[fi]*time)
    
    #create gaussian window
    gaus_win = exp(-time**2/(2*(6/(2*pi*frequencies[fi]))**2))
    
    #create wavelet by multiplying our sine wave by the gaussian window
    wavelet_family[fi,:] = sinewave*gaus_win
    
    #this could be done in one line
#     wavelet_family[fi,:] = exp(2*1j*pi*frequencies[fi]*time) * exp(-time**2/(2*(6/(2*pi*frequencies[fi]))**2))

#plot some of our wavelet family
fig=plt.figure()
plt.subplots(2,1)
plt.subplot(211)
plt.plot(time,np.real(wavelet_family[::int(random.rand()*30),:].T))

plt.subplot(212)
plt.plot(time,np.real(wavelet_family[30,:]))
plt.plot(time,np.imag(wavelet_family[30,:]),'r:')
plt.title("real and imaginary parts of one wavelet")
plt.legend(["real","imaginary"])
plt.tight_layout()
```


```python
fig=plt.figure(figsize=(6,6))
plt.imshow(np.real(wavelet_family), 
           extent=[time[0], time[-1], frequencies[0], frequencies[-1]],
           aspect="auto",
           cmap=plt.get_cmap("hot"),
           origin="lower")

plt.xlabel("time (s)")
_=plt.ylabel("frequency (Hz)")
```

### Figure 12.5


```python
from numpy.fft import fft, ifft #import fft functions for ease of use
from scipy import signal as sig
#EEG data from one trial (electrode FCz)
eegdata = np.squeeze(EEGdata[46,:,9])
EEGpnts = data["EEG"][0,0]["pnts"][0,0] #number of points in EEG data
EEGtimes = data["EEG"][0,0]["times"][0]
EEGsrate = float(data["EEG"][0,0]["srate"][0])

#create wavelet
time = np.arange(-1,1 + 1/EEGsrate,1/EEGsrate)
f = 6 #frequency in Hz
sine_wave = exp(2*1j*pi*f*time)

#compute gaussian
s=4.5/(2*pi*f)
gaussian_win = exp(-time**2/(2*s**2))
#window the sinewave by a gaussian to create complex morlet wavelet
wavelet = sine_wave * gaussian_win

#half of wavelet size, useful for chopping off edges after convolution
halfwaveletsize = len(wavelet)//2

#convolve with data
n_conv = len(wavelet) + EEGpnts - 1 #number of points in our convolution

fft_w = fft(wavelet,n_conv)
fft_e = fft(eegdata,n_conv)
#convolution theorem -- convolution = pointwise multiplication in frequency-space
ift = ifft(fft_e*fft_w,n_conv)*sqrt(s)/10 #sqrt(s)/20 is empirical scaling factor (sqrt(s)/10 in the book)
wavelet_conv_data = np.real(ift[halfwaveletsize: - halfwaveletsize]) #take middle portion of convolution
wavelet_conv_data = np.real(ift[halfwaveletsize:-halfwaveletsize]) #take middle portion of convolution

#create a filter to apply to data
nyquist = EEGsrate/2
transition_width = 0.2 #percent
filter_low = 4 #Hz
filter_high = 8 #Hz

ffrequencies = np.array([0 ,filter_low*(1-transition_width), 
                         filter_low, filter_high, filter_high*(1+transition_width), nyquist])/nyquist
ideal_response = np.array([0, 0, 1, 1, 0, 0])

#there doesn't seem to be a python equivalent to MATLAB's firls function,
#so I am going to use butterworth filter as a close approximation. 

b, a = sig.butter(5, np.array([filter_low*(1-transition_width),filter_high*(1+transition_width)])/nyquist,btype="bandpass")
eeg_4to8 = sig.filtfilt(b, a, eegdata, padlen=150)

plt.plot(EEGtimes,eegdata)
plt.plot(EEGtimes,wavelet_conv_data,'r')
plt.plot(EEGtimes,eeg_4to8,'g')
plt.axis([-200,1200,-40,40])
plt.xlabel("time (ms)")
plt.ylabel("voltage (mV)")
_=plt.legend(["raw","wavelet conv","band-passed"])
```

### Figure 12.6


```python
time = np.arange(-1,1+1/EEGsrate,1/EEGsrate)
n_conv = EEGpnts + len(time) -1
n2p1 = n_conv//2 + 1

f = 6 #hz
s = 6/(2*pi*f)
wavelet = exp(2*pi*1j*f*time) * exp(-time**2/(2*s**2))
halfwaveletsize = len(wavelet)//2

eegdata = np.squeeze(EEGdata[46,:,9])

plt.figure()
plt.subplot(311)
plt.plot(EEGtimes,eegdata)
plt.xlim([-500,1200])
plt.title("raw")

plt.subplot(323)
fft_w = fft(wavelet,n_conv)
hz    = np.linspace(0,EEGsrate/2.,n2p1)
plt.plot(hz,np.absolute(fft_w[:n2p1])/np.max(np.absolute(fft_w[:n2p1])),'b')

fft_e = fft(eegdata,n_conv)
plt.plot(hz,np.absolute(fft_e[:n2p1])/np.max(np.absolute(fft_e[:n2p1])),'g')
plt.axis([0,40,0,1.05])
plt.title("individual power spectra")

plt.subplot(324)
plt.plot(hz,np.absolute(fft_e[:n2p1]*np.absolute(fft_w[:n2p1])))
plt.xlim([0, 40])
plt.title("convolved power spectrum")

plt.subplot(313)
plt.plot(EEGtimes,eegdata)
ift = ifft(fft_e*fft_w,n_conv)*sqrt(s)/10 #sqrt(s)/20 is empirical scaling factor (sqrt(s)/10 in the book)
plt.plot(EEGtimes,np.real(ift[halfwaveletsize:-halfwaveletsize]),'r')
plt.title("wavelet filtered")

plt.tight_layout()
```

### Figure 12.7


```python
#create 10Hz wavelet kernel

time = np.arange(-(EEGpnts/EEGsrate/2),EEGpnts/EEGsrate/2 + 1/EEGsrate,1/EEGsrate)
f = 10. #hz
s = 4/(2*pi*f) #sd of gaussian
wavelet = cos(2*pi*f*time) * exp(-time**2/(2*s**2))

#signal is one sine cycle
timeS = np.arange(0,1/f + 1/EEGsrate,1/EEGsrate)
signal = sin(2*pi*f*timeS)

#zeropad the signal
zz = np.zeros(EEGpnts//2 - len(timeS)//2)
signal = np.concatenate([zz,signal,zz])

plt.figure(figsize=(6,6))

#plot waves
plt.subplot(321)
plt.plot(wavelet,'r')
plt.xlim(200, len(time) - 200)
plt.title("wavelet")

plt.subplot(323)
plt.plot(signal)
plt.xlim([200, len(time)-200])
plt.title("1 cycle of signal")

plt.subplot(325)
plt.plot(np.convolve(wavelet,signal,mode="same"),'purple')
plt.axis([200,len(time)-200,-12,12])
plt.title("convolved wavelet and signal")

#plot the dot products at selected phase lags
plt.subplot(322)
plt.plot(wavelet[int(100/f)-2-1:],'r')
plt.plot(signal)
plt.xlim([200,len(time)-200])
plt.title("dot product: " + str( np.fix(np.sum(wavelet[int(100/f)-2-1:]*signal[:-int(100/f)+3]))))
plt.legend(["wavelet","signal"])

plt.subplot(324)
plt.plot(wavelet[int(2.3*100/f)-2-1:],'r')
plt.plot(signal)
plt.xlim([200,len(time)-200])
plt.title("dot product: " + str( np.fix(np.sum(
                wavelet[int(2.3*100/f)-2-1:]*signal[:-int(2.3*100/f)+3]))))

plt.subplot(326)
plt.plot(wavelet,'r')
plt.plot(signal)
plt.xlim([200,len(time)-200])
plt.title("dot product: " + str( np.fix(np.sum(
                wavelet*signal))))

plt.tight_layout()
```

6d189dfb61d1d9cd243722a2ee0656b80ddab083
ipynb
Jupyter Notebook
notebooks/chapter12.ipynb
danielemarinazzo/Analyzing_Neural_Time_Series
838d521492665b5a26b41ce44c58637a8fc4967a
MIT
notebooks/chapter12.ipynb
danielemarinazzo/Analyzing_Neural_Time_Series
838d521492665b5a26b41ce44c58637a8fc4967a
MIT
notebooks/chapter12.ipynb
danielemarinazzo/Analyzing_Neural_Time_Series
838d521492665b5a26b41ce44c58637a8fc4967a
MIT
2020-10-25T05:02:32.000Z
2020-10-25T05:02:32.000Z
Qwen/Qwen-72B
1. YES
2. YES