Scribd
Upload
14 views

temp_imf

The document contains a Python script that performs Empirical Mode Decomposition (EMD) on a simulated signal. It computes instantaneous frequency and amplitude of intrinsic mode functions (IMFs) derived from the signal using the Hilbert transform. The script also visualizes the original signal, IMFs, and their instantaneous frequencies through various plots.

Uploaded by

gunjanc080
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as TXT, PDF, TXT or read online on Scribd
14 views

temp_imf

The document contains a Python script that performs Empirical Mode Decomposition (EMD) on a simulated signal. It computes instantaneous frequency and amplitude of intrinsic mode functions (IMFs) derived from the signal using the Hilbert transform. The script also visualizes the original signal, IMFs, and their instantaneous frequencies through various plots.

Uploaded by

gunjanc080
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as TXT, PDF, TXT or read online on Scribd
You are on page 1/ 3

import numpy as np

import matplotlib.pyplot as plt


from scipy.signal import find_peaks
from scipy import interpolate
from scipy.signal import hilbert
from scipy.interpolate import CubicSpline
import pandas as pd
import math
import random
import emd

def instantaneous_frequency(signal, sampling_rate):


# Compute the analytic signal using the Hilbert transform
analytic_signal = hilbert(signal)

# Compute phase of the analytic signal


phase = np.angle(analytic_signal)

# Unwrap the phase to handle phase wrapping


phase = np.unwrap(phase)

# Compute instantaneous frequency


delta_phase = np.diff(phase)
instantaneous_frequency = (delta_phase / (2 * np.pi)) * sampling_rate

return instantaneous_frequency

def instantaneous_amplitude(imfs):
amplitude = []
for imf in imfs:
# Apply Hilbert transform to obtain analytic signal
analytic_signal = hilbert(imf)
# Calculate magnitude of the analytic signal (instantaneous amplitude)
amp = np.abs(analytic_signal)
amplitude.append(amp)
return amplitude

pi=np.pi
max_iter =6
# imf=np.zero(100)
g=2560
t = np.linspace(0, 1, g)
# print(t)
signal=np.zeros(len(t))
# tol = 1

# df1=pd.read_csv(f'D:\\pythonProject1\\venv\\pipe.csv')
# d = np.array(df1)
# # print(d)
# t=d[:,0]
# print(t[-1])
# leak=d[:,1]
# no_leak=d[:,2]
# print(leak)
# print(no_leak)
# plt.plot(t,leak)
# plt.plot(t,no_leak)
# plt.legend()
# plt.show()
gr=1
# for k in range(len(t)):
# signal[k]=5*np.sin(2*pi**50*t[k])+6*np.sin(2*pi*80*t[k])
# while abs(gr) >0.01:
# for l in range(len(t)):
# dt=random.randint(1,2)
# if dt//2==0:
# signal[l] = np.cos(10*pi*t[l])**1+(-1) ** dt * random.random()**1
# # signal[l]=0
# else:
# signal[l] =np.sin(1*pi*t[l])**1*(-1) ** dt * random.random()**1
#
# gr=np.mean(signal)
# print(gr)

peak_freq = 15
sample_rate = 256
seconds = 10
noise_std = .4
x = emd.simulate.ar_oscillator(peak_freq, sample_rate, seconds,noise_std=noise_std,
random_seed=42, r=0.96)[:, 0]
signal = x*1e-8

fig1, axs = plt.subplots(max_iter+2)


fig2, axs1 = plt.subplots(max_iter)
fig3, axs2 = plt.subplots(max_iter)
# fig1.set_figwidth(5)
# fig1.set_figheight(3)
fig1.suptitle('Epirical Mode Decomposition')
fig2.suptitle('Instantaneous Frequency')
fig3.suptitle('Instantaneous Frequency Histograms')
# plt.subplot(1,1,1)

# plt.legend

df1=pd.read_csv('hist_temp.csv')
d1 = np.array(df1)

f = d1[:,1]

g=len(f)
t = np.linspace(0, 1, g)
print(f)
x = f
axs[0].set_title('Signal',x=1.05,y=0.5)
axs[0].plot(t,x,color='black',label='signal')
axs[max_iter+1].plot(t,x,color='black',label='signal')

d=np.zeros(len(t))
sum=0
for i in range(max_iter):
av = 1

while abs(av) > 1:


peaks, _ = find_peaks(x)
troughs, _ = find_peaks(-x)

P=x[peaks]
T=x[troughs]
x_p = t[peaks]
x_t= t[troughs]

fp = CubicSpline(x_p, P)
xp_new = np.linspace(0, 1, g)
yp_new = fp(xp_new)

ft = CubicSpline(x_t, T)
xt_new = np.linspace(0, 1, g)
yt_new = ft(xt_new)

m=0.5*(yp_new+yt_new)

d = x-m
av=np.mean(d)
print(av)
if abs(av) > 1:
x=d
else:
imf = d
# plt.subplot(i+1,i+1,1)
col = (np.random.random(), np.random.random(), np.random.random())
# Compute Instantaneous Frequency (IF)
if_signal = instantaneous_frequency(imf, sample_rate)
# plt.figure(figsize=(10, 6))
# plt.plot(t[:-1], if_signal) # Plot instantaneous frequency over time
# plt.xlabel('Time')
axs1[i + 0].set_title(f'IMF{i + 1}', x=1.05, y=0.5)
axs1[i + 0].plot(t[:-1], if_signal, c=col, label=f'imf {i + 1}')
axs2[i + 0].set_title(f'IMF{i + 1}', x=1.05, y=0.5)
axs2[i + 0].hist(imf,bins=100, color=col, label=f'imf {i + 1}')
# plt.hist(imf, bins=30, color='skyblue', edgecolor='black')
# plt.ylabel('Instantaneous Frequency (Hz)')
# plt.title('Instantaneous Frequency (IF) in Hilbert Spectral
Analysis')
axs[i+1].set_title(f'IMF {i + 1}',x=1.05,y=0.5)
axs[i+1].plot(xt_new, d,c=col, label=f'imf {i + 1}')
# plt.legend()

x=m
sum=sum+d
# plt.tight_layout()
# plt.plot(t,signal,label='signal')
# plt.plot(xt_new, sum, label='summed up imfs')
# plt.legend()
xt_new = np.linspace(0, 1, g)
col = (np.random.random(), np.random.random(), np.random.random())
axs[max_iter+1].set_title('summed up',x=1.05,y=0.5)
axs[max_iter+1].plot(xt_new, sum,c=col, label='summed up imfs')

plt.show()

You might also like