%load_ext autoreload %autoreload 2 %reload_ext autoreload

%matplotlib inline import IPython import _pickle as pickle import numpy as np import scipy as sp import scipy.signal as spsig import pandas as pd import matplotlib.pyplot as plt

#Resting EMG # Sampling rate is 3.3kHz F_0 = 3.3 * 1e3 #load the data csv_path = 'Fouzia resting.txt' emg_rest = pd.read_csv(csv_path, usecols=[0], header=None)[0] # create a time axis t_ax = np.linspace(0, len(emg_rest) / F_0, len(emg_rest)) # plot the signal plt.plot(t_ax,emg_rest) #plt.xlim([t_ax[0], t_ax[-1]]) plt.title("Resting EMG Data") plt.xlabel("Time (s)") plt.ylabel("Amplitude") plt.show()

#Resting FFT resting = np.fft.fft(emg_rest) #f_ax = np.linspace(0, 1, len(emg_side)) f_ax_rest = np.linspace(-F_0/2, F_0/2, len(emg_rest)) fft_csv_roll_rest = np.fft.fftshift(resting) plt.plot(f_ax_rest, np.abs(fft_csv_roll_rest), 'b') plt.xlabel(r"Frequency") plt.ylabel("Amplitude") plt.title("Fourier transform of resting EMG") #plt.ylim(0, 50000) #plt.xlim(0,20)

#Max Clenching EMG # Sampling rate is 3.3kHz F_0 = 3.3 * 1e3 #load the data csv_path = 'Fouzia max clenching.txt' emg = pd.read_csv(csv_path, usecols=[0], header=None)[0] # create a time axis t_ax = np.linspace(0, len(emg) / F_0, len(emg)) # plot the signal plt.plot(t_ax,emg) plt.xlim([t_ax[], t_ax[-1]]) plt.title("Max clenching") plt.xlabel("Time (s)") plt.ylabel("Amplitude") plt.show()

#Head rotation # Sampling rate is 3.3kHz F_0 = 3.3 * 1e3 #load the data csv_path = 'Fouzia head rotation.txt' emg = pd.read_csv(csv_path, usecols=[0], header=None)[0] # create a time axis t_ax = np.linspace(0, len(emg) / F_0, len(emg)) # plot the signal plt.plot(t_ax,emg) #plt.xlim([t_ax[0], t_ax[-1]]) plt.title("Head rotation") plt.xlabel("Time (s)") plt.ylabel("Amplitude") plt.show()

#Rest + Grinding bouts # Sampling rate is 3.3kHz F_0 = 3.3 * 1e3 #load the data csv_path = 'Fouzia sets of rest and grinding.txt' emg_side = pd.read_csv(csv_path, usecols=[0], header=None)[0] # create a time axis t_ax = np.linspace(0, len(emg_side) / F_0, len(emg_side)) # plot the signal plt.plot(t_ax,emg_side) plt.xlim([t_ax[15000], t_ax[30000]]) plt.title("Rest and Grinding") plt.xlabel("Time (s)") plt.ylabel("Amplitude") plt.show()

#Resting and grinding FFT grinding_one_side = np.fft.fft(emg_side) #f_ax = np.linspace(0, 1, len(emg_side)) f_ax = np.linspace(-F_0/2, F_0/2, len(emg_side)) fft_csv_roll_grind = np.fft.fftshift(grinding_one_side) plt.plot(f_ax, np.abs(fft_csv_roll_grind), 'b') plt.xlabel(r"Frequency") plt.ylabel("Amplitude") plt.title("Fourier transform of resting+grinding") plt.ylim(0, 200000) plt.xlim(0,20)

#Sets of resting and grinding both sides # Sampling rate is 3.3kHz F_0 = 3.3 * 1e3 #load the data csv_path = 'Fouzia sets of resting and grinding both sides.txt' emg_both = pd.read_csv(csv_path, usecols=[0], header=None)[0] # create a time axis t_ax = np.linspace(0, len(emg_both) / F_0, len(emg_both)) # plot the signal plt.plot(t_ax,emg_both) #plt.xlim([t_ax[0], t_ax[-1]]) plt.title("Resting + Grinding Both Sides signal") plt.xlabel("Time (s)") plt.ylabel("Amplitude") plt.show()

#Resting and grinding FFT grinding_both_side = np.fft.fft(emg_both) #f_ax = np.linspace(0, 1, len(emg_both)) f_ax_both = np.linspace(-F_0/2, F_0/2, len(emg_both)) fft_csv_roll_grind_both = np.fft.fftshift(grinding_both_side) plt.plot(f_ax_both, np.abs(fft_csv_roll_grind_both), 'b') plt.xlabel(r"Frequency") plt.ylabel("Amplitude") plt.title("Fourier transform of EMG grinding on both sides") plt.ylim(0, 300000) plt.xlim(0,100)

#Swallowing # Sampling rate is 3.3kHz F_0 = 3.3 * 1e3 #load the data csv_path = 'Fouzia swallowing.txt' emg = pd.read_csv(csv_path, usecols=[0], header=None)[0] # create a time axis t_ax = np.linspace(0, len(emg) / F_0, len(emg)) # plot the signal plt.plot(t_ax,emg) #plt.xlim([t_ax[0], t_ax[-1]]) plt.title("Swallowing") plt.xlabel("Time (s)") plt.ylabel("Amplitude") plt.show()

#Talking # Sampling rate is 3.3kHz F_0 = 3.3 * 1e3 #load the data csv_path = 'Fouzia talking.txt' emg = pd.read_csv(csv_path, usecols=[0], header=None)[0] # create a time axis t_ax = np.linspace(0, len(emg) / F_0, len(emg)) # plot the signal plt.plot(t_ax,emg) #plt.xlim([t_ax[0], t_ax[-1]]) plt.title("Talking") plt.xlabel("Time (s)") plt.ylabel("Amplitude") plt.show()

#Yawning # Sampling rate is 3.3kHz F_0 = 3.3 * 1e3 #load the data csv_path = 'Fouzia yawning.txt' emg = pd.read_csv(csv_path, usecols=[0], header=None)[0] # create a time axis t_ax = np.linspace(0, len(emg) / F_0, len(emg)) # plot the signal plt.plot(t_ax,emg) #plt.xlim([t_ax[0], t_ax[-1]]) plt.title("Yawning") plt.xlabel("Time (s)") plt.ylabel("Amplitude") plt.show()

#Rest (clean) # Sampling rate is 3.3kHz F_0 = 3.3 * 1e3 #load the data csv_path = 'fouziarestnoblink.txt' emg = pd.read_csv(csv_path, usecols=[0], header=None)[0] # create a time axis t_ax = np.linspace(0, len(emg) / F_0, len(emg)) # plot the signal plt.plot(t_ax[1000:-1],emg[1000:-1]) #plt.xlim([t_ax[1000], t_ax[-1]]) plt.title("Rest (no blinking") plt.xlabel("Time (s)") plt.ylabel("Amplitude") plt.show()

#Clenching (5 off 5 on x 6), no blink # Sampling rate is 3.3kHz F_0 = 3.3 * 1e3 #load the data csv_path = 'fouziaclenching5on5off.txt' emg = pd.read_csv(csv_path, usecols=[0], header=None)[0] # create a time axis t_ax = np.linspace(0, len(emg) / F_0, len(emg)) # plot the signal plt.plot(t_ax[1000:-15000],emg[1000:-15000]) #plt.xlim([t_ax[1000], t_ax[-1]]) plt.title("On and off clenching") plt.xlabel("Time (s)") plt.ylabel("Amplitude") plt.show()

#Grinding (15 off 30 sec on), no blink # Sampling rate is 3.3kHz F_0 = 3.3 * 1e3 #load the data csv_path = 'fouzia15srest30secondsgrind.txt' emg = pd.read_csv(csv_path, usecols=[0], header=None)[0] # create a time axis t_ax = np.linspace(0, len(emg) / F_0, len(emg)) # plot the signal plt.plot(t_ax[1000:-15000],emg[1000:-15000]) #plt.xlim([t_ax[1000], t_ax[-1]]) plt.title("Continuous teeth grinding") plt.xlabel("Time (s)") plt.ylabel("Amplitude") plt.show()

#Grinding (10 off 10 sec on), no blink # Sampling rate is 3.3kHz F_0 = 3.3 * 1e3 #load the data csv_path = 'fouzia10xgrind10s.txt' emg = pd.read_csv(csv_path, usecols=[0], header=None)[0] # create a time axis t_ax = np.linspace(0, len(emg) / F_0, len(emg)) # plot the signal plt.plot(t_ax[1000:-1000],emg[1000:-1000]) #plt.xlim([t_ax[1000], t_ax[-1]]) plt.title("Continuous teeth grinding") plt.xlabel("Time (s)") plt.ylabel("Amplitude") plt.show()

# enter your code here F_0 = 3.3e3

csv_path = 'Head side to side.txt' csv = pd.read_csv(csv_path, usecols=[0], header=None)[0] csv_path2 = 'Clenching (5 off 5 on).txt' csv2 = pd.read_csv(csv_path2, usecols=[0], header=None)[0] # create a time axis t_ax = np.linspace(0, len(csv) / F_0, len(csv)) t_ax2 = np.linspace(0, len(csv2) / F_0, len(csv2)) # plot the signal plt.plot(t_ax, csv) plt.title("Head side to side") plt.xlabel("Time (s)") plt.ylabel("Amplitude") # plot the signal plt.figure() plt.plot(t_ax2, csv2) #plt.ylim([50, 90]) plt.title("Clenching") plt.xlabel("Time (s)") plt.ylabel("Amplitude")

# Resting data #Resting EMG # Sampling rate is 3.3kHz F_0 = 3.3 * 1e3 #load the data csv_path = 'Rest 2 (THIS ONE).txt' emg_rest = pd.read_csv(csv_path, usecols=[0], header=None)[0] # create a time axis t_ax = np.linspace(0, len(emg_rest) / F_0, len(emg_rest)) # plot the signal plt.plot(t_ax,emg_rest) plt.xlim([5, t_ax[-1]]) plt.ylim([40, 100]) plt.title("Resting EMG") plt.xlabel("Time (s)") plt.ylabel("Amplitude") plt.ylim([50,100]) plt.show()

# Maximal #Resting EMG # Sampling rate is 3.3kHz F_0 = 3.3 * 1e3 #load the data csv_path = 'Max clench.txt' emg_rest = pd.read_csv(csv_path, usecols=[0], header=None)[0] # create a time axis t_ax = np.linspace(0, len(emg_rest) / F_0, len(emg_rest)) # plot the signal plt.plot(t_ax,emg_rest) plt.xlim([5, t_ax[-1]]) plt.ylim([40, 100]) plt.title("Max clench EMG") plt.xlabel("Time (s)") plt.ylabel("Amplitude") plt.ylim([50,100]) plt.show()

#Ethan Data # Sampling rate is 3.3kHz F_0 = 3.3 * 1e3 #load the data csv_path = 'Clenching (5 off 5 on).txt' emg_ethan_grind = pd.read_csv(csv_path, usecols=[0], header=None)[0] # create a time axis t_ax = np.linspace(0, len(emg_ethan_grind) / F_0, len(emg_ethan_grind)) # plot the signal plt.plot(t_ax,emg_ethan_grind) #plt.xlim([0, 20]) plt.ylim([50, 100]) plt.title("Clenching trials EMG") plt.xlabel("Time (s)") plt.ylabel("Amplitude") plt.show()

def binary_filter(y, x0, sig0, N, sigma2): ''' A binary filter that takes y - data x0 - initial condition on mean sig0 - initial condition on varianace N - number of samples in data sigma2 - parameter sigma^2 and outputs the filtered estimate p_n|n ''' # initialize p, and set first probability (zero index) according to function p = np.zeros(DATLENGTH) p[0] = np.exp(x0)/(1 + np.exp(x0)) # implement the recursion xcurrent = x0 sigcurrent = sig0 for n in range(1,DATLENGTH - 1): xlast = xcurrent #xlast is xn-1_n-1 is the same as xn_xn-1 by choice siglast = sigcurrent #siglast is sign-1_n-1 is the same as sign_n-1 sigmid = siglast + sigma2 #sigmid is sign_n-1 pmid = np.exp(xlast)/(1 + np.exp(xlast)) #pmid is pn_n-1 sigcurrent = 1/(1/sigmid + N * pmid*(1 - pmid)) xcurrent = xlast + sigcurrent * (y[n] - N * pmid) p[n] = np.exp(xcurrent)/(1 + np.exp(xcurrent)) return p

DATLENGTH = np.size(csv2) #constant

grindbf1 = binary_filter(csv2, -3.5, 0, 50, 10**-1) plt.figure(figsize=(12, 8), dpi=80) plt.plot(t_ax2[0:50000], grindbf1[0:50000])

BSS_data = pickle.load(open('BSSdata.p','rb'),encoding='latin1') t = BSS_data['t'] X = BSS_data['X'] Fs = BSS_data['Fs']

import scipy.io.wavfile as wav wav.write('p1_mix1.wav',Fs,X[0,:]) wav.write('p1_mix2.wav',Fs,X[1,:])

IPython.display.Audio("p1_mix1.wav")

IPython.display.Audio("p1_mix2.wav")