# -*- 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, butter
from scipy.fftpack import fft
from scipy import signal as sg
#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]
# -----------------------------------------------------------------------------
# Statistic Variables
# -----------------------------------------------------------------------------
def _variance(x):
miu = 0.0
vari = 0.0
for i in range(x.shape[0]):
miu =+ x[i]
miu = miu/len(x)
for i in range(x.shape[0]):
vari = vari + (x[i]-miu)**2
vari = vari/len(x)
return vari
def _StandarD(x):
miu = 0.0
vari = 0.0
for i in range(x.shape[0]):
miu += x[i]
miu = miu/len(x)
for i in range(x.shape[0]):
vari = vari + (x[i]-miu)**2
vari = vari/len(x)
_Std = vari **(.5)
return _Std
def _CV(x):
# Coefficient of Variation
var = _variance(x)
std = _StandarD(x)
return 100*(var/std)
# -----------------------------------------------------------------------------
# PDS
# -----------------------------------------------------------------------------
def vec_nor(x):
"""
Normalize the amplitude of a vector from -1 to 1
"""
nVec=np.zeros(len(x)); # Initializate derivate vector
nVec = np.divide(x, max(x))
nVec=nVec-np.mean(nVec);
nVec=np.divide(nVec,np.max(nVec));
return nVec
def running_sum(x):
"""
Running Sum Algorithm of an input signal is y[n]= x[n] + y[n-1]
"""
y = np.zeros(len(x))
for i in range(len(x)):
y[i] = x[i] + y[i-1]
return vec_nor(y)
def derivate_1 (x):
"""
Derivate of an input signal as y[n]= x[n] - x[n-1]
"""
y=np.zeros(len(x)); # Initializate derivate vector
for i in range(len(x)):
y[i]=x[i]-x[i-1];
return vec_nor(y)
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
def derivate_positive (x):
"""
Derivate of an input signal as y[n]= x[n+1]- x[n-1]
for all values where the signal is positive
"""
lenght=x.shape[0] # Get the length of the vector
y=np.zeros(lenght); # Initializate derivate vector
for i in range(lenght-1):
if x[i]>0:
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 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
def shannonE_vector_1 (x):
"""
Shannon energy of an input signal
"""
N = x.shape[0]
# Se = -(1/N) *
mu = -(x**2)*np.log(x**2)/x.shape[0]
y = -(((x**2)*np.log(x**2)) - mu)/np.std(x)
return y
def E_VS_LF (pcgFFT1, vTfft1, on):
"""
Energy of PCG Vibratory Spectrum in low Frequencies (E_VS_LF)
(frequency components, frequency value vector, on = on percentage or not)
According with [1] The total vibratory spectrum can be divided into 7 bands.
This is a modification of this 7 bands
1. 0-5Hz, 2. 5-25Hz; 3. 25-60Hz; 4. 60-120Hz; 5. 120-400Hz
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()
# 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)
x = np.array([xAll, x1, x2, x3, x4, x5])
if (on == 'percentage'):
x = 100*(x/x[0])
return x
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 features_1(data_P1, fs):
# Defining the Vibratory Frequency Bands
bVec = [0.01, 120, 240, 350, 425, 500, 999]
# Initializing Vectors
band_matrix = np.zeros((len(bVec),len(data_P1))) # Band Matrix
power_matrix = np.zeros((len(bVec),int(1+(len(data_P1)/2)))) # Band Matrix
freqs_matrix = np.zeros((len(bVec),int(1+(len(data_P1)/2)))) # Band Matrix
SePCG = np.zeros(len(bVec)-1) # Shannon Energy in each Band
SePWR = np.zeros(len(bVec)-1) # Shannon Energy in each Band
for i in range(len(bVec)-1):
band_matrix[i,:] = butter_bp_fil(data_P1, bVec[i], bVec[i+1], fs)
freqs_matrix[i,:], power_matrix[i,:] = sg.periodogram(band_matrix[i,:], fs, scaling = 'density')
SePCG[i] = Energy_value(band_matrix[i,:])
SePWR[i] = Energy_value(power_matrix[i,:])
SePWR = 1*np.log10(SePWR)
SePCG = 1*np.log10(SePCG)
return SePWR, SePCG
# -----------------------------------------------------------------------------
# Filter Processes
# -----------------------------------------------------------------------------
def recursive_moving_average_F(X, Fs, M):
"""
The recursive Moving Average Filter its an algorithm to implement the typical
moving average filter more faster. The algorithm is written as:
y[i] = y[i-1] + x[i+p] - x[i-q]
p = (M - 1)/2, q = p + 1
Ref: Udemy Course "Digital Signal Processing (DSP) From Ground Up™ in Python",
available in: https://www.udemy.com/python-dsp/
---------------------------------------------------------------------------
X: input signal
Fs: Sampling Frequency
M: number of points in M range of time the moving average
"""
p = int(((M*Fs)-1)/2)
q = p + 1
Y = np.zeros(len(X))
for i in range(len(X)):
Y[i] = Y[i-1] + X[i-p] - X[i-q]
return vec_nor(Y)
def butter_bp_coe(lowcut, highcut, fs, order=1):
"""
Butterworth passband filter coefficients b and a
Ref:
[1] https://timsainb.github.io/spectrograms-mfccs-and-inversion-in-python.html
[2] https://gist.github.com/kastnerkyle/179d6e9a88202ab0a2fe
"""
nyq = 0.5 * fs
low = lowcut / nyq
high = highcut / nyq
b, a = butter(order, [low, high], btype='band')
return b, a
def butter_bp_fil(data, lowcut, highcut, fs, order=1):
"""
Butterworth passband filter
Ref:
[1] https://timsainb.github.io/spectrograms-mfccs-and-inversion-in-python.html
[2] https://gist.github.com/kastnerkyle/179d6e9a88202ab0a2fe
"""
b, a = butter_bp_coe(lowcut, highcut, fs, order=order)
y = lfilter(b, a, data)
return vec_nor(y)
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 bandpass_filter(data, lowcut, highcut, signal_freq, filter_order):
"""
Method responsible for creating and applying Butterworth filter.
:param deque data: raw data
:param float lowcut: filter lowcut frequency value
:param float highcut: filter highcut frequency value
:param int signal_freq: signal frequency in samples per second (Hz)
:param int filter_order: filter order
:return array: filtered data
"""
nyquist_freq = 0.5 * signal_freq
low = lowcut / nyquist_freq
high = highcut / nyquist_freq
b, a = butter(filter_order, [low, high], btype="band")
y = lfilter(b, a, data)
return y
def FpassBand_1(X,Fs, f1, f2):
"""
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)
#
taps = firwin(100, [f1, f2], pass_zero=False)
# taps = firwin(N, L_cutoff_hz/nyq_rate, window=('kaiser', beta))
# taps_2 = firwin(N, H_cutoff_hz/nyq_rate, pass_zero=True)
# Use lfilter to filter x with the FIR filter.
X_pb= lfilter(taps, 1.0, X)
# X_pb= lfilter(taps_2, 1.0, X_l)
return X_pb[N-1:]
def FpassBand_2(X, f1, f2, Fs):
Y = FhighPass(X, Fs, f1)
Y = FlowPass(X, Fs, f2)
return Y
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:]
# -----------------------------------------------------------------------------
# Segmentation By Running Sum Algorithm & Filters
# -----------------------------------------------------------------------------
def seg_RSA1(x, fs):
# # Appliying Running Sum Algorithm of PCG filtered from 0.01Hz to 1kHz
F_x = running_sum(vec_nor(butter_bp_fil(x, 0.01, 50, fs)))
# Smoothing the signal by filtering from 0.01Hz to 5Hz
F_x = butter_bp_fil(F_x,0.01, 2, fs)
# Appliying 1st derivative to indentify slope sign changes
xx = derivate_1(F_x)
#
xx = butter_bp_fil(xx, 0.01,2, fs)
# Transforming positives to 1 & negatives to -1
xxS = np.sign(xx)
return xxS, F_x
# -----------------------------------------------------------------------------
# Segmentation By Derivatives
# -----------------------------------------------------------------------------
def seg_Der1(x, fs):
# Segmenting in an specific frequency band
pcgPeaks, peaks, allpeaks = PDP(x, fs)
# -----------------------------------------------------------------------------
# Time processing
dT = 0.4 # [1] mean S1 duration 122ms, mean S2 duration 92ms
timeV = []
timeV.append(0)
pointV = []
pointV.append(0)
segmV = np.zeros(len(x))
for i in range(len(pcgPeaks)-1):
if pcgPeaks[i]>0.5:
timeV.append(i/fs) # Gives the time when a peak get found
pointV.append(i) # Gives the time when a peak get found
if (pointV[-1]/fs)-(pointV[-2]/fs)> dT:
segmV[pointV[-2]:pointV[-1]] = 0.4 # Marc a diastolic segment
else:
segmV[pointV[-2]:pointV[-1]] = 0.6 # Marc a systolic segment
return segmV, pcgPeaks
def seg_Der2(x, fs):
# Appliying Running Sum Algorithm of PCG filtered from 0.01Hz to 1kHz
F_x = running_sum(vec_nor(butter_bp_fil(x, 0.01, 999, fs)))
# Smoothing the signal by filtering from 0.01Hz to 5Hz
F_x = butter_bp_fil(F_x,0.01, 5, fs)
# Time to be represented in samples
time_samples = 0.5
# Number of samples to move over the signal
mC = int(time_samples * fs)
peaks = np.zeros(len(F_x))
p = sg.find_peaks(F_x, distance=mC)
# Defining peaks as +1
for i in range (len(p[0][:])):
peaks[p[0][i]] = 1
return peaks
# -----------------------------------------------------------------------------
# Peak Detection
# -----------------------------------------------------------------------------
def PDP(Xf, samplerate):
"""
Peak Detection Process
"""
timeCut = samplerate*0.25 # Time to count another pulse
vCorte = 0.6 # Amplitude threshold
Xf = vec_nor(Xf) # Normalize signal
dX = derivate_positive(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
'''
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 Xf[i] > 0:
if dX[i]>0:
positive[i] = Xf[i];
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]==Xf[i] 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 (Xf[i] > vCorte and p==0):
p = i
points1[i] = Xf[i]
else:
points1[i] = 0
if (p+timeCut < i):
p = 0
return points1, points, positive[0:(len(positive)-1)]
# -----------------------------------------------------------------------------
# Peak Detection 2
# -----------------------------------------------------------------------------
def PDP_2(Xf, samplerate):
"""
Peak Detection Process
"""
timeCut = samplerate*0.25 # Time to count another pulse
vCorte = 0.6 # Amplitude threshold
#Xf = vec_nor(Xf) # Normalize signal
dX = np.diff(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
'''
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 Xf[i] > 0:
if dX[i]>0:
positive[i] = Xf[i];
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]==Xf[i] 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 (Xf[i] > vCorte and p==0):
p = i
points1[i] = Xf[i]
else:
points1[i] = 0
if (p+timeCut < i):
p = 0
return points1, points, positive[0:(len(positive)-1)]
# -----------------------------------------------------------------------------
# Peak Detection 3 - Janko Slavic
# -----------------------------------------------------------------------------
def findpeaks(data, spacing=1, limit=None):
"""
Janko Slavic peak detection algorithm and implementation.
https://github.com/jankoslavic/py-tools/tree/master/findpeaks
Finds peaks in `data` which are of `spacing` width and >=`limit`.
:param ndarray data: data
:param float spacing: minimum spacing to the next peak (should be 1 or more)
:param float limit: peaks should have value greater or equal
:return array: detected peaks indexes array
"""
len = data.size
x = np.zeros(len + 2 * spacing)
x[:spacing] = data[0] - 1.e-6
x[-spacing:] = data[-1] - 1.e-6
x[spacing:spacing + len] = data
peak_candidate = np.zeros(len)
peak_candidate[:] = True
for s in range(spacing):
start = spacing - s - 1
h_b = x[start: start + len] # before
start = spacing
h_c = x[start: start + len] # central
start = spacing + s + 1
h_a = x[start: start + len] # after
peak_candidate = np.logical_and(peak_candidate, np.logical_and(h_c > h_b, h_c > h_a))
ind = np.argwhere(peak_candidate)
ind = ind.reshape(ind.size)
if limit is not None:
ind = ind[data[ind] > limit]
return ind
# -----------------------------------------------------------------------------
# Fast Fourier Transform
# -----------------------------------------------------------------------------
def fft_k(data, samplerate, showFrequency):
'''
Fast Fourier Transform using fftpack from scipy
Ref: https://docs.scipy.org/doc/scipy/reference/tutorial/fftpack.html
'''
# FFT Full Vector 'k' coefficients
pcgFFT = fft(data)
# FFT positives values
short_pcgFFT = 2.0/np.size(data) * np.abs(pcgFFT[0:(np.size(data)//2)])
#short_pcgFFT = 2.0/np.size(data) * np.abs(pcgFFT[(np.size(data)//2):None]) vector selected from the middle to the end
# Vector of frequencies (X-axes)
vTfft = np.linspace(0.0, 1.0/(2.0*(1/samplerate)), np.size(data)//2)
# find the value closest to a value
idx = (np.abs(vTfft-showFrequency)).argmin()
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
'''
# FFT Full Vector 'k' coefficients
pcgFFT = fft(data)
# FFT positives values from the middle to the end (to evoid the interference at beginning)
short_pcgFFT = 2.0/np.size(data) * np.abs(pcgFFT[0:(np.size(data)//2)])
#short_pcgFFT = 2.0/np.size(data) * np.abs(pcgFFT[(np.size(data)//2):None]) vector selected from the middle to the end
# Vector of frequencies (X-axes)
vTfft = np.linspace(0.0, 1.0/(2.0*(1/samplerate)), np.size(data)//2)
# find the value closest to a value
idx = (np.abs(vTfft-showFrequency)).argmin()
return vec_nor(short_pcgFFT[0:idx]), vTfft[0:idx]
# -----------------------------------------------------------------------------
# PCG Audio Pre-Processing
# -----------------------------------------------------------------------------
def pre_pro_audio_PCG(x, fs):
# Ensure having a Mono sound
if len(x.shape)>1:
# select the left size
x = x[:,0]
# Resampling Audio PCG to 2k Hz
Frs = 2000
Nrs = int(Frs*(len(x)/fs))
if fs > Frs:
x = sg.resample(x, Nrs)
return vec_nor(x), Frs
# -----------------------------------------------------------------------------
# Pre-Process: Signal basic information (SBI)
# -----------------------------------------------------------------------------
def pre_pro_basicInfo_PCG(x, fs):
# find the time duration of the sound
t_sound = 1/fs*len(x)
# make a vector time for the sound
vec_time = np.linspace(0, t_sound, len(x))
return t_sound, vec_time
Datasets
Models
