Diff of /normalizer.py [000000] .. [881eaf]

Switch to side-by-side view

--- a
+++ b/normalizer.py
@@ -0,0 +1,207 @@
+
+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()