# -*- coding: utf-8 -*-
"""
Personal Processing Functions 1 (ppfunctions_1)
Created on Mon Apr 9 11:48:37 2018
@author: Kevin MAchado
Module file containing Python definitions and statements
"""
# Libraries
import numpy as np
from scipy.signal import kaiserord, lfilter, firwin
from scipy.fftpack import fft
#import peakutils # Librery to help in peak detection
# Functions
# Normal energy = np.sum(pcgFFT1**2)
# Normalized average Shannon Energy = sum((l**2)*np.log(l**2))/l.shape[0]
# Third Order Shannon Energy = sum((l**3)*np.log(l**3))/l.shape[0]
def vec_nor(x):
"""
Normalize the amplitude of a vector from -1 to 1
"""
lenght=np.size(x) # Get the length of the vector
xMax=max(x); # Get the maximun value of the vector
nVec=np.zeros(lenght); # Initializate derivate vector
nVec = np.divide(x, xMax)
nVec=nVec-np.mean(nVec);
nVec=np.divide(nVec,np.max(nVec));
return nVec
def derivate (x):
"""
Derivate of an input signal as y[n]= x[n+1]- x[n-1]
"""
lenght=x.shape[0] # Get the length of the vector
y=np.zeros(lenght); # Initializate derivate vector
for i in range(lenght-1):
y[i]=x[i-1]-x[i];
return y
# -----------------------------------------------------------------------------
# Energy
# -----------------------------------------------------------------------------
def Energy_value (x):
"""
Energy of an input signal
"""
y = np.sum(x**2)
return y
def E_VS (pcgFFT1, vTfft1, on):
"""
Energy of PCG Vibratory Spectrum
(frequency components, frequency value vector, on = on percentage or not)
According with [1] The total vibratory spectrum can be divided into 7 bands:
1. 0-5Hz, 2. 5-25Hz; 3. 25-120Hz; 4. 120-240Hz; 5. 240-500Hz; 6. 500-1000Hz; 7. 1000-2000Hz
The PCG signal producess vibrations in the spectrum between 0-2k Hz.
[1] Abbas, Abbas K. (Abbas Khudair), Bassam, Rasha and Morgan & Claypool Publishers Phonocardiography signal processing. Morgan & Claypool Publishers, San Rafael, Calif, 2009.
"""
c1 = (np.abs(vTfft1-5)).argmin()
c2 = (np.abs(vTfft1-25)).argmin()
c3 = (np.abs(vTfft1-120)).argmin()
c4 = (np.abs(vTfft1-240)).argmin()
c5 = (np.abs(vTfft1-500)).argmin()
c6 = (np.abs(vTfft1-1000)).argmin()
c7 = (np.abs(vTfft1-2000)).argmin()
# All vector energy
xAll = Energy_value(pcgFFT1)
# Procesando de 0.01-5 Hz
pcgFFT_F1 = pcgFFT1[0:c1]
x1 = Energy_value(pcgFFT_F1)
# Procesando de 5-25 Hz
pcgFFT_F2 = pcgFFT1[c1:c2]
x2 = Energy_value(pcgFFT_F2)
# Procesando de 25-120 Hz
pcgFFT_F3 = pcgFFT1[c2:c3]
x3 = Energy_value(pcgFFT_F3)
# Procesando de 120-240 Hz
pcgFFT_F4 = pcgFFT1[c3:c4]
x4 = Energy_value(pcgFFT_F4)
# Procesando de 240-500 Hz
pcgFFT_F5 = pcgFFT1[c4:c5]
x5 = Energy_value(pcgFFT_F5)
# Procesando de 500-1000 Hz
pcgFFT_F6 = pcgFFT1[c5:c6]
x6 = Energy_value(pcgFFT_F6)
# Procesando de 1000-2000 Hz
pcgFFT_F7 = pcgFFT1[c6:c7]
x7 = Energy_value(pcgFFT_F7)
x = np.array([xAll, x1, x2, x3, x4, x5, x6, x7])
if (on == 'percentage'):
x = 100*(x/x[0])
return x
def shannonE_value (x):
"""
Shannon energy of an input signal
"""
y = sum((x**2)*np.log(x**2))/x.shape[0]
return y
def shannonE_vector (x):
"""
Shannon energy of an input signal
"""
mu = -(x**2)*np.log(x**2)/x.shape[0]
y = -(((x**2)*np.log(x**2)) - mu)/np.std(x)
return y
# -----------------------------------------------------------------------------
# Filter Processes
# -----------------------------------------------------------------------------
def Fpass(X,lp):
"""
Fpass is the function to pass the coefficients of a filter trough a signal'
"""
llp=np.size(lp) # Get the length of the lowpass vector
x=np.convolve(X,lp); # Disrete convolution
x=x[int(llp/2):-int(llp/2)];
x=x-(np.mean(x));
x=x/np.max(x);
y=vec_nor(x); # Vector Normalizing
return y
def FpassBand(X,hp,lp):
"""
FpassBand is the function that develop a pass band filter of the signal 'x' through the
discrete convolution of this 'x' first with the coeficients of a High Pass Filter 'hp' and then
with the discrete convolution of this result with a Low Pass Filter 'lp'
"""
llp=np.shape(lp) # Get the length of the lowpass vector
llp=llp[0]; # Get the value of the length
lhp=np.shape(hp) # Get the length of the highpass vector
lhp=lhp[0]; # Get the value of the length
x=np.convolve(X,lp); # Disrete convolution
x=x[int(llp/2):-int(llp/2)];
x=x-(np.mean(x));
x=x/np.max(x);
y=np.convolve(x,hp); # Disrete onvolution
y=y[int(lhp/2):-int(lhp/2)];
y=y-np.mean(y);
y=y/np.max(y);
x=np.convolve(y,lp); # Disrete convolution
x=x[int(llp/2):-int(llp/2)];
x=x-(np.mean(x));
x=x/np.max(x);
y=np.convolve(x,hp); # Disrete onvolution
y=y[int(lhp/2):-int(lhp/2)];
y=y-np.mean(y);
y=y/np.max(y);
y=vec_nor(y); # Vector Normalizing
return y
def FpassBand_1(X,Fs,H_cutoff_hz, L_cutoff_hz):
"""
Ref: http://scipy-cookbook.readthedocs.io/items/FIRFilter.html
http://lagrange.univ-lyon1.fr/docs/scipy/0.17.1/generated/scipy.signal.firwin.html
FpassBand_1 is a function to develop a passband filter using 'firwin'
and 'lfilter' functions from the "Scipy" library
"""
# The Nyquist rate of the signal.
nyq_rate = Fs / 2.0
# The desired width of the transition from pass to stop,
# relative to the Nyquist rate. We'll design the filter
# with a 5 Hz transition width.
width = 5.0/nyq_rate
# The desired attenuation in the stop band, in dB.
ripple_db = 60.0
# Compute the order and Kaiser parameter for the FIR filter.
N, beta = kaiserord(ripple_db, width)
# Use firwin with a Kaiser window to create a lowpass FIR filter.
taps = firwin(N, L_cutoff_hz/nyq_rate, window=('kaiser', beta))
taps_2 = firwin(N, H_cutoff_hz/nyq_rate, pass_zero=False)
# Use lfilter to filter x with the FIR filter.
X_l= lfilter(taps, 1.0, X)
X_pb= lfilter(taps_2, 1.0, X_l)
return X_pb[N-1:]
def FhighPass(X, Fs, H_cutoff_hz):
"""
Ref: http://scipy-cookbook.readthedocs.io/items/FIRFilter.html
http://lagrange.univ-lyon1.fr/docs/scipy/0.17.1/generated/scipy.signal.firwin.html
FhighPass is a function to develop a highpass filter using 'firwin'
and 'lfilter' functions from the "Scipy" library
"""
# The Nyquist rate of the signal.
nyq_rate = Fs / 2.0
# The desired width of the transition from pass to stop,
# relative to the Nyquist rate. We'll design the filter
# with a 5 Hz transition width.
width = 5.0/nyq_rate
# The desired attenuation in the stop band, in dB.
ripple_db = 60.0
# Compute the order and Kaiser parameter for the FIR filter.
N, beta = kaiserord(ripple_db, width)
# Use firwin with a Kaiser window to create a lowpass FIR filter.
taps_2 = firwin(N, H_cutoff_hz/nyq_rate, pass_zero=False)
# Use lfilter to filter x with the FIR filter.
X_h= lfilter(taps_2, 1.0, X)
return X_h[N-1:]
def FlowPass(X, Fs, L_cutoff_hz):
"""
Ref: http://scipy-cookbook.readthedocs.io/items/FIRFilter.html
http://lagrange.univ-lyon1.fr/docs/scipy/0.17.1/generated/scipy.signal.firwin.html
FlowPass is a function to develop a lowpass filter using 'firwin'
and 'lfilter' functions from the "Scipy" library
"""
# The Nyquist rate of the signal.
nyq_rate = Fs / 2.0
# The desired width of the transition from pass to stop,
# relative to the Nyquist rate. We'll design the filter
# with a 5 Hz transition width.
width = 5.0/nyq_rate
# The desired attenuation in the stop band, in dB.
ripple_db = 60.0
# Compute the order and Kaiser parameter for the FIR filter.
N, beta = kaiserord(ripple_db, width)
# Use firwin with a Kaiser window to create a lowpass FIR filter.
taps = firwin(N, L_cutoff_hz/nyq_rate, window=('kaiser', beta))
# Use lfilter to filter x with the FIR filter.
X_l= lfilter(taps, 1.0, X)
return X_l[N-1:]
# -----------------------------------------------------------------------------
# Peak Detection
# -----------------------------------------------------------------------------
def PDP(Xf, samplerate):
"""
Peak Detection Process
"""
timeCut = samplerate*0.25 # time to count another pulse
vCorte = 0.6
Xf = vec_nor(Xf)
dX=derivate(Xf); # Derivate of the signal
dX=vec_nor(dX); # Vector Normalizing
size=np.shape(Xf) # Rank or dimension of the array
fil=size[0]; # Number of rows
positive=np.zeros((1,fil+1)); # Initializating Vector
positive=positive[0]; # Getting the Vector
points=np.zeros((1,fil)); # Initializating the Peak Points Vector
points=points[0]; # Getting the point vector
points1=np.zeros((1,fil)); # Initializating the Peak Points Vector
points1=points1[0]; # Getting the point vector
vCorte = 0.6
'''
FIRST! having the positives values of the slope as 1
And the negative values of the slope as 0
'''
for i in range(0,fil):
if dX[i]>0:
positive[i]=1;
else:
positive[i]=0;
'''
SECOND! a peak will be found when the ith value is equal to 1 &&
the ith+1 is equal to 0
'''
for i in range(0,fil):
if (positive[i]==1 and positive[i+1]==0):
points[i]=Xf[i];
else:
points[i]=0;
'''
THIRD! Define a minimun Peak Height
'''
p=0;
for i in range(0,fil):
if (points[i] > vCorte and p==0):
p=i
points1[i]=Xf[i]
else:
points1[i]=0;
if (p+timeCut < i):
p=0
return points1
# -----------------------------------------------------------------------------
# Fast Fourier Transform
# -----------------------------------------------------------------------------
def fft_k(data, samplerate, showFrequency):
'''
Fast Fourier Transform
Ref: https://docs.scipy.org/doc/scipy/reference/tutorial/fftpack.html
'''
pcgFFT = fft(data) # FFT Full Vector
short_pcgFFT = 2.0/np.size(data) * np.abs(pcgFFT[0:np.size(data)//2]) # FFT positives values
vTfft = np.linspace(0.0, 1.0/(2.0*(1/samplerate)), np.size(data)//2) # Vector of frequencies (X-axes)
idx = (np.abs(vTfft-showFrequency)).argmin() # find the value closest to a value
return short_pcgFFT[0:idx], vTfft[0:idx]
def fft_k_N(data, samplerate, showFrequency):
'''
Normalized Fast Fourier Transform
Ref: https://docs.scipy.org/doc/scipy/reference/tutorial/fftpack.html
'''
pcgFFT = fft(data) # FFT Full Vector
short_pcgFFT = 2.0/np.size(data) * np.abs(pcgFFT[0:np.size(data)//2]) # FFT positives values
vTfft = np.linspace(0.0, 1.0/(2.0*(1/samplerate)), np.size(data)//2) # Vector of frequencies (X-axes)
idx = (np.abs(vTfft-showFrequency)).argmin() # find the value closest to a value
return vec_nor(short_pcgFFT[0:idx]), vTfft[0:idx]
Datasets
Models
