import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy
import scipy.io as sio
import os
import sys
root_path = os.path.join(os.path.dirname(__file__), ".." )
sys.path.append(root_path)
signal_proc_path = os.path.join(root_path, "Signal_processing_toolbox/src" )
sys.path.append(signal_proc_path)
[docs]
def welch_example():
"""
Spectral analysis of an EEG signal using welch method
**Welch method:**
Cut the signal in n window and compute fft on each of them.
Imagine a time-domain signal represented as a sinusoidal wave.
At the core of the Welch method, you divide this signal into partially overlapping segments.
For each segment, you apply a window to mitigate edge effects. T
hen, instead of taking the FFT of the entire signal, you take the FFT of each modified segment and average them.
The idea is that this overlapping approach reduces the variance of the estimation,
which can be particularly useful when the signal changes over time.
**Advantage over fft:**
More adapt if signal characteristics change over time.
E.g Noise increase/ decrease/ change of shape over time.
**Example:**
.. image:: _static/images/EEG_signal.png
.. image:: _static/images/welch_versus_fft.png
"""
file_data_path = os.path.join(root_path, r"data\EEGrestingState.mat")
if not os.path.isfile(file_data_path):
print("\n\tFile data could not be found. Please check that you have access to it\n")
return
matdat = sio.loadmat(file_data_path)
eegdata = matdat['eegdata'][0]
srate = matdat['srate'][0]
# time vector
N = len(eegdata)
timevec = np.arange(0,N)/srate
# plot the data
plt.figure()
plt.title("EEG signal")
plt.plot(timevec,eegdata,'k')
plt.xlabel('Time (seconds)')
plt.ylabel('Voltage (\muV)')
plt.show()
## one big FFT (not Welch's method)
# "static" FFT over entire period, for comparison with Welch
eegpow = np.abs( scipy.fftpack.fft(eegdata)/N )**2
hz = np.linspace(0,srate/2,int(np.floor(N/2)+1))
# create Hann window
winsize = int( 2*srate ) # 2-second window
hannw = .5 - np.cos(2*np.pi*np.linspace(0,1,winsize))/2
# number of FFT points (frequency resolution)
nfft = srate*100
f, welchpow = scipy.signal.welch(eegdata,fs=srate,window=hannw,
nperseg=winsize,noverlap=winsize/4,
nfft=nfft)
plt.figure()
plt.title("Welch versus fft")
plt.plot(f,welchpow, label="welch")
plt.plot(hz,np.log(eegpow[0:len(hz)]), label="fft")
plt.xlim([0,40])
plt.xlabel('frequency [Hz]')
plt.ylabel('Power')
plt.legend()
plt.show()
[docs]
def spectrogramm_with_morlet_wavelet():
r"""
Spectogramm by doing a convolution between the raw signal and complex morlet wavelet:
**Method:**
Create multiple complex morlet wavelet from different frequency.
The frequency range of the morlet wavelets correspond to the frequency axis of the spectogramm.
**Example:**
**EEG raw signal**
.. image:: _static/images/eeg_signal_for_spetogramm.png
1. **Create complex morlet wavelet**
- Example of one complex molet wavelet (one for the range of the frequency)
.. image:: _static/images/complex_wavelet_morvelet.png
- All the morlet wavelet shape/amplitude represented on a colormap
.. image:: _static/images/time_frequency_morlet_wavelet.png
2. **Convolutions**
The most efficient way is to loop on the different complex morlet wavelet
and multiply fft(raw_signal) with fft(complex_morlet_wavelet[i])
3. **Plot Spectogramm**
.. image:: _static/images/spectrogramm_morlet_wavelet.png
**Comparison Spectrogramm with Morlet Wavelet vs FFT:**
- **scipy.signal.spectrogram(data, fs) (FFT-based Spectrogram)**:
**Method:** Utilizes the Fast Fourier Transform (FFT) to compute the frequency content of the signal.
**Features:** The spectrogram generated by this method has a fixed frequency resolution and a time resolution based on the size of the temporal window used.
**Advantages:** It is a fast and efficient method to obtain a time-frequency representation of the signal, suitable for relatively stationary signals.
**Disadvantages:** Less suitable for non-stationary signals or short-duration events.
- **Time-Frequency Analysis with Morlet Wavelet Transform**:
**Method:** Uses a family of Morlet wavelets to compute the frequency contribution at different temporal moments.
**Features:** Offers variable time and frequency resolution, which is particularly useful for non-stationary signals. It can better adapt to frequency variations in the signal over time.
**Advantages:** Suitable for non-stationary signals, provides representation with better temporal localization for short-duration events.
**Disadvantages:** May be computationally more intensive than the FFT method.
In summary, the choice between these two approaches depends on the type of signal you are analyzing and your specific goals. If your signal is relatively stationary and you need a quick analysis, the FFT method may be appropriate. On the other hand, if your signal is non-stationary with significant frequency variations, the Morlet wavelet transform may be more suitable despite potentially higher computational complexity.
"""
##### load in data #####
file_data_path = os.path.join(root_path, r"data\data4TF.mat")
if not os.path.isfile(file_data_path):
print("\n\tFile data could not be found. Please check that you have access to it\n")
return
braindat = sio.loadmat(file_data_path)
timevec = braindat['timevec'][0]
srate = braindat['srate'][0]
data = braindat['data'][0]
# plot the signal
plt.figure()
plt.plot(timevec,data)
plt.xlabel('Time (s)'), plt.ylabel('Voltage (\muV)')
plt.title('EEG Time-domain signal')
plt.show()
##### create complex Morlet wavelets #####
# wavelet parameters
nfrex = 50 # 50 frequencies
frex = np.linspace(8,70,nfrex)
fwhm = .2 # full-width at half-maximum in seconds
# time vector for wavelets
wavetime = np.arange(-2,2,1/srate)
# initialize matrices for wavelets
wavelets = np.zeros( (nfrex,len(wavetime)) ,dtype=complex)
# create complex Morlet wavelet family
for wi in range(0,nfrex):
# Gaussian
gaussian = np.exp( -(4*np.log(2)*wavetime**2) / fwhm**2 )
# complex Morlet wavelet
wavelets[wi,:] = np.exp(1j*2*np.pi*frex[wi]*wavetime) * gaussian
###### show the wavelets #####
plt.figure()
plt.title("Complex Morlet Wavelet")
plt.plot(wavetime,np.real(wavelets[10,:]),label='Real part')
plt.plot(wavetime,np.imag(wavelets[10,:]),label='Imag part')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.xlim([-.5, .5])
plt.legend()
plt.show()
plt.figure()
plt.pcolormesh(wavetime,frex,np.real(wavelets))
plt.xlabel('Time (s)'), plt.ylabel('Frequency (Hz)')
plt.title('Real part of wavelets')
plt.xlim([-.5,.5])
plt.show()
####### run convolution using spectral multiplication #####
# convolution parameters
nconv = len(timevec) + len(wavetime) - 1 # M+N-1
halfk = int( np.floor(len(wavetime)/2) )
# Fourier spectrum of the signal
dataX = scipy.fftpack.fft(data,nconv)
# initialize time-frequency matrix
tf = np.zeros( (nfrex,len(timevec)) )
# convolution per frequency
for fi in range(0,nfrex):
# FFT of the wavelet
waveX = scipy.fftpack.fft(wavelets[fi,:],nconv)
# amplitude-normalize the wavelet
waveX = waveX / np.abs(np.max(waveX))
# convolution
convres = scipy.fftpack.ifft( waveX*dataX )
# trim the "wings"
convres = convres[halfk-1:-halfk]
# extract power from complex signal
tf[fi,:] = np.abs(convres)**2
##### plot the results #####
plt.figure()
plt.subplot(211)
plt.plot(timevec,data)
plt.xlabel('Time (s)'), plt.ylabel('Voltage (\muV)')
plt.title('EEG Time-domain signal')
plt.subplot(212)
plt.pcolormesh(timevec,frex,tf,vmin=0,vmax=1e3)
plt.xlabel('Time (s)'), plt.ylabel('Frequency (Hz)')
plt.title('Spectrogramm (Morlet Wavelet)')
plt.xlim(-0.5, 1.5)
plt.show()
# plt.figure()
# plt.pcolormesh(timevec,frex,tf,vmin=0,vmax=1e3)
# plt.xlabel('Time (s)'), plt.ylabel('Frequency (Hz)')
# plt.title('Spectrogramm (Morlet Wavelet)')
# plt.show()
##### compare with normal fft spectrogramm ####
frexfft,timefft,pwrfft = scipy.signal.spectrogram(data,srate, nperseg=10, noverlap=4)
plt.figure()
plt.subplot(311)
plt.plot(timevec,data)
plt.xlabel('Time (s)'), plt.ylabel('Voltage (\muV)')
plt.title('EEG Time-domain signal')
plt.subplot(312)
plt.pcolormesh(timevec,frex,tf,vmin=0,vmax=1e3)
plt.xlabel('Time (s)'), plt.ylabel('Frequency (Hz)')
plt.title('Spectrogramm (Morlet Wavelet)')
plt.subplot(313)
plt.title("Spectrogramm (fft)")
plt.pcolormesh(timefft-0.5,frexfft,pwrfft,vmin=0,vmax=9)
plt.xlabel('Time (s)'), plt.ylabel('Frequency (Hz)')
plt.ylim(10, 70)
plt.xlim(-0.5, 1.5)
plt.show()
if __name__ == "__main__":
welch_example()
spectrogramm_with_morlet_wavelet()