Page MenuHomec4science
Files F111637497

Utils_EMD.py
No OneTemporary
Actions

File Metadata

Created
Sat, May 3, 23:56

Utils_EMD.py
View Options

# -*- coding: utf-8 -*-
"""
Created on Wed Feb 1 23:21:53 2023
@author: srpv
"""
from scipy import signal
import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.patches as mpatches
import pandas as pd
import numpy as np
import torch, math
import pandas as pd
import numpy as np
import os
import emd
def EMDplot(x, sample_rate,path, num_imf):
t0=0
N=len(x)
dt=1/sample_rate
t = np.arange(0, N) * dt + t0
imfs = emd.sift.sift(x,max_imfs=num_imf)
# imfs=imfs.detach().cpu().numpy().squeeze()
imfs=imfs.T
num_imf=num_imf+1
plt.rcParams.update(plt.rcParamsDefault)
plt.rcParams['xtick.labelsize'] = 15
plt.rcParams['ytick.labelsize'] = 15
plt.figure(figsize=(5,12))
for i in range(num_imf):
plt.subplot(num_imf,1,i+1)
# min=np.min(x.ravel())
# max=np.max(x.ravel())
plt.plot(t,x.ravel(),'black',color='0.6')
plt.plot(t,imfs[i],'r')
# plt.ylim([min*1.1,max*1.1])
plt.ylabel('IMF '+str(i+1),fontsize=15)
if i != num_imf-1:
plt.tick_params(
axis='x', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
bottom=False, # ticks along the bottom edge are off
top=False, # ticks along the top edge are off
labelbottom=False)
if i == num_imf-1:
plt.ticklabel_format(axis='x', style='sci',scilimits=(0,0))
plt.xlabel('Time (s)',fontsize=15)
if i == 0:
plt.title('Raw signal vs IMFs',fontsize=15)
# plt.legend(['Rawsignal','IMF(s)'],loc='lower right',frameon=False, bbox_to_anchor=(1.0, 1.1))
plt.tight_layout()
graphname='IMF.png'
plt.savefig(os.path.join(path, graphname), bbox_inches='tight',dpi=100)
plt.show()
plt.close()
plt.rcParams.update(plt.rcParamsDefault)
plt.rcParams['xtick.labelsize'] = 15
plt.rcParams['ytick.labelsize'] = 15
fig, axs = plt.subplots(nrows=len(imfs), ncols=1, sharex=True,
figsize=(6, 12),dpi=100)
for i in range(num_imf):
ax = axs.flat[i]
Frequencyplot_(imfs,sample_rate,i,ax)
if i == 0:
skyblue = mpatches.Patch(color='#657e39', label='0-150 kHz')
# ax.legend(handles=[skyblue])
red = mpatches.Patch(color='#aa332d', label='150-300 kHz')
# ax.legend(handles=[red])
yellow = mpatches.Patch(color='#f0a334', label='300-450 kHz')
# ax.legend(handles=[yellow])
green = mpatches.Patch(color='#0080ff', label='450-600 kHz')
# ax.legend(handles=[green])
cyan = mpatches.Patch(color='#b05aac', label='600-750 kHz')
ax.set_title('IMF vs Frequency',fontsize=15)
if i != num_imf-1:
# ax.set_xticks([])
plt.setp( ax.get_xticklabels(), visible=False)
# plt.tick_params(top='off', bottom='off', left='off', right='off', labelleft='off', labelbottom='off')
if i == num_imf-1:
# plt.setp( ax.get_xticklabels(), visible=True)
plt.xlabel('Frequency(Hz)',fontsize=15)
# ax.legend(handles=[skyblue,red,yellow,green,cyan], bbox_to_anchor=(1, 0.35))
plt.tight_layout()
graphname='IMF_FFT.png'
plt.savefig(os.path.join(path, graphname), bbox_inches='tight',dpi=100)
plt.show()
def Frequencyplot_(rawspace,sample_rate,choicewindow,ax):
data_new=rawspace[choicewindow]
# Define window length (4 seconds)
win = 0.1 * sample_rate
freqs, psd = signal.welch(data_new, sample_rate, nperseg=win)
# Plot the power spectrum
# sns.set(font_scale=1.5, style='white')
ax.plot(freqs, psd, color='k', lw=0.001)
# ax.plot(freqs, psd)
sec1, sec2 = 0, 150000
sec3, sec4 = 150000, 300000
sec5, sec6 = 300000, 450000
sec7, sec8 = 450000, 600000
sec9, sec10 = 600000, 750000
# Find intersecting values in frequency vector
idx_delta1 = np.logical_and(freqs >= sec1, freqs <= sec2)
idx_delta2 = np.logical_and(freqs >= sec3, freqs <= sec4)
idx_delta3 = np.logical_and(freqs >= sec5, freqs <= sec6)
idx_delta4 = np.logical_and(freqs >= sec7, freqs <= sec8)
idx_delta5 = np.logical_and(freqs >= sec9, freqs <= sec10)
ax.set_ylabel('PSD \n(dB/Hz)',fontsize=15)
ax.yaxis.tick_right()
ax.yaxis.set_label_position("right")
ax.set_xlim([0, sample_rate//2])
ax.fill_between(freqs, psd, where=idx_delta1, color='#657e39')
ax.fill_between(freqs, psd, where=idx_delta2, color='#aa332d')
ax.fill_between(freqs, psd, where=idx_delta3, color='#f0a334')
ax.fill_between(freqs, psd, where=idx_delta4, color='#0080ff')
ax.fill_between(freqs, psd, where=idx_delta5, color='#b05aac')
ax.ticklabel_format(axis='x', style='sci',scilimits=(0,0))
ax.ticklabel_format(axis='y', style='sci',scilimits=(0,0))
ax.get_yaxis().get_offset_text().set_position((1.15,0))

Event Timeline

c4science · Help