normalizer.py

import heartbeat as hb
import matplotlib.pyplot as plt 
import numpy as np 
from scipy import fftpack
from scipy.signal import butter, lfilter, freqz
import math

def z_norm_01(result):
    """
    Normalize Data. This creats a mean of 0, and normalizes 
    all values between 0 and 1 by shifting negative values 
    upward
    """
    result_mean = result.mean()
    result -= result_mean
    if min(result)<0.0:
        result-=min(result)
    result /= max(abs(result))
    return result

def z_norm_b(result):
    """
    Normalize Data. This creats a mean of 0, and a
    std of 1
    """
    result_mean = result.mean()
    result_std = result.std()
    result -= result_mean
    result /= result_std
    return result

def z_norm(result):
    """
    Normalize Data. This fits
    all values between 0 and 1. 
    """
    result = (result-min(result))/(max(result)-min(result))
    return result

def z_norm2(result):
    """
    Normalize Data. This creats a mean of 0, and fits
    all values between -1 and 1. 
    """
    result_mean = result.mean()
    result_std = result.std()
    result -= result_mean
    result /= max(abs(result))
    return result

def dynamic_threshold(signal, fs):
    mvg_avg_percent= np.arange(5,100,5).tolist() #[0 100]
    valid_percent = []
    interval_SD=[]
  
    for percent in mvg_avg_percent: #Detect peaks with all percentages, append results to list 'rrsd'
        mov_avg, peaklist= peak_finder(signal, percent, fs)
        RR_list, bpm = get_HR(peaklist)
        print(percent,bpm)
        RR_list_STD= np.std(RR_list)
        
        interval_SD.append([RR_list_STD, bpm, percent])
    for interval,HR,perc in interval_SD: #Test list entries and select valid measures
        if ((interval > 1) and ((HR > 30) and (HR < 130))):
            valid_percent.append([interval, perc])

    mov_avg, peaklist= peak_finder(signal, min(valid_percent, key = lambda t: t[0])[1], fs) #Detect peaks with 'ma_perc' that goes with lowest rrsd
    return mov_avg, peaklist

def peak_finder(signal,mvg_perc, fs):
    """
    1. Get moving average of the data
    2. Replace NaN's with moving average
    3. 5% Amplification to prevent heart rate interference
    4. for each point in signal
        current_mean = point in moving avg
        if point < current mean and len(window) =0 #NO R Complex Activity 
            get new mean avg point 
        else if point > current mean: #region of interest  
            include point in window 
            get new mean avg point
        else 
            peak = position of the point on the X-axis
            get new mean avg point

    """
    width = fs*.1
    mvg_perc= (1+mvg_perc/100)
    mov_avg = signal.rolling(int(width)).mean() #Calculate moving average
    
    #Impute where moving average function returns NaN, which is the beginning of the signal where x hrw
    avg_hr = np.mean(signal)
    mov_avg = [avg_hr if math.isnan(x) else x for x in mov_avg]
    mov_avg = [x*mvg_perc for x in mov_avg] #For now we raise the average by 20% to prevent the secondary heart contraction from interfering, in part 2 we will do this dynamically
    
    window = []
    peaklist = []
    listpos = 0 #We use a counter to move over the different data columns
    for datapoint in signal:
        rollingmean = mov_avg[listpos] #Get local mean
        if (datapoint <= rollingmean) and (len(window) <= 1): #If no detectable R-complex activity -> do nothing
            listpos += 1
        elif (datapoint > rollingmean): #If signal comes above local mean, mark ROI
            window.append(datapoint)
            listpos += 1
        else: #If signal drops below local mean -> determine highest point
            
            beatposition = listpos - len(window) + (window.index(max(window))) #Notate the position of the point on the X-axis
            peaklist.append(beatposition) #Add detected peak to list
            window = [] #Clear marked ROI
            listpos += 1
    
    return mov_avg, peaklist 
  
class filt: 

    def fft_plot(self,sig,label,color):
        """
        Generates a FFT or fast fourier transform 
        plot for the input signal.

        Input: Signal, Label (filtered/not), Color of Figure
        Output: Histogram of Frequencies 

        """
        f_s=360
        X = fftpack.fft(sig)
        freqs = fftpack.fftfreq(len(sig)) * f_s
        plt.plot(freqs, np.abs(X))
        plt.xlabel('Frequency in Hertz [Hz]')
        plt.ylabel('Frequency Domain (Spectrum) Magnitude')
        plt.xlim(0, f_s / 2)
        plt.ylim(-5, 110)

    def butter_lowpass(self,cutoff, fs, order=5):
        """
        Butterworth lowpass filter:: Allows only signal below 
        certain cutoff to pass through. Order is a constant value
        and is part of the filter arguments. 

        Return b,a inputs for filter 

        """
        nyq = 0.5 * fs
        normal_cutoff = cutoff / nyq
        b, a = butter(order, normal_cutoff, btype='low', analog=False)
        return b, a

    def butter_lowpass_filter(self,data, cutoff, fs=360, order=5):
        """
        Low pass filter for the signal data w/r to a particular
        cut off frequency. 

        """

        b, a = self.butter_lowpass(cutoff, fs, order=order)
        y = lfilter(b, a, data)
        return y

    def low_pass_filter_plot(self,patient,cutoff,fs=360,order=5):
        """
        Generates a 3 figure subplot
        FFT, Illustration of Frequency Response, 
        and a plot of signal vs unfilt signal 
        for that patient. 

        Inputs are patient, cuttoff freq

        Outputs are the plots. 


        """

        plt.figure(figsize=(20,10))
        sig,notes=hb.get_patient_data(patient,norm=True)

        #FFT PLOT
        plt.subplot(221)
        self.fft_plot(sig,label='Orig Sig',color='b')
        self.fft_plot(self.butter_lowpass_filter(sig,cutoff),label='Filt Sig',color='r')
        plt.title('FFT Patient: {}'.format(patient))
        plt.xlabel('Frequency [Hz]')
        plt.xlim(0,cutoff+30)

        #Low pass Filter Frequency Response 
        b, a = self.butter_lowpass(cutoff, fs, order=order)
        y = lfilter(b, a, sig)
        w, h = freqz(b, a, worN=8000)
        plt.subplot(222)
        plt.plot(0.5*fs*w/np.pi, np.abs(h), 'b')
        plt.plot(cutoff, 0.5*np.sqrt(2), 'ko')
        plt.axvline(cutoff, color='k')
        plt.xlim(0, .1*fs)
        plt.title("Lowpass Filter Frequency Response")
        plt.xlabel('Frequency [Hz]')
        plt.grid()

        #Signal Plot 
        plt.subplot(212)
        plt.plot(sig, 'b-', label='data')
        plt.plot(self.butter_lowpass_filter(sig,cutoff), 'g-', linewidth=2, label='filtered data')
        plt.xlabel('Time [sec]')
        plt.xlim(0,1000)
        plt.ylim(-.5,1.1)
        plt.grid()
        plt.legend()