demo.py

import os
import json
import h5py
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
import scipy.io

from tensorflow.keras import Sequential
from tensorflow.keras.layers import Conv1D, MaxPooling1D, Dropout, Flatten, Dense, InputLayer


def grabLoudestCh(song16ch, method="mean"):
    if method == "mean":
        chIdx = np.argmax(np.mean(np.abs(song16ch), axis=1))
    elif method == "max":
        chIdx = np.argmax(np.max(np.abs(song16ch), axis=1))
    return song16ch[chIdx, :], chIdx


def smooth(x, window_len=15):
    s = np.r_[x[window_len-1:0:-1], x, x[-1:-window_len:-1]]
    w = np.ones(window_len, 'd')
    y = np.convolve(w/w.sum(), s, mode='valid')
    z = y[(window_len-1)//2 + np.arange(len(x))]
    return z


def normalizePulseAndCut(signal, pulseWidth=165, signWindow=np.array([10,0]), smoothWindow=15):
    signalDur = len(signal)
    if signalDur < pulseWidth:
        raise ValueError('Signal should be longer than pulseWidth')
    pulsePad = (pulseWidth - 1)//2
    smoothedPower = smooth(signal**2, smoothWindow)
    mIdx = np.argmax(smoothedPower[pulsePad:(signalDur-pulsePad)])
    oldIm = mIdx + pulsePad
    oldI0 = oldIm-pulsePad
    oldI1 = oldIm+pulsePad+1
    newPulse = signal[oldIm-pulsePad : oldIm+pulsePad+1]
    newPulse = newPulse / np.linalg.norm(newPulse)
    if np.mean(newPulse[pulsePad-signWindow[0]:pulsePad-signWindow[1]]) < 0:
        newPulse = -newPulse
    return newPulse, oldIm, oldI0, oldI1


def load_pulse_classifier(out_path):
    if os.path.exists(out_path):
        print('Loading converted model from', out_path)
        model = tf.keras.models.load_model(out_path)
        return model
    else:
        print('model not found at', out_path)


def segPulseAux(pulseWindw, framei1, lastPulseLocation, p,clf):
	# Move samples around to get the window of desired length with overlapping samples from the previous frame
    windowi0 = framei1 - p['pulseWindowSize'] # first index of the window

    # take the channel with the largest mean amplitude
    window1ch,loudestMic = grabLoudestCh(pulseWindw, p['pulse_chSelection'])

    # Process the signal, such that the pulse candidate is centered, cropped, and normalized
    [normPcand, mIdx, i0, i1] = normalizePulseAndCut(window1ch, p['pulseWidth'])
    # Compute the index of the pulse peak in the unit of samples, from the beginning of the experiment
    pulseIdx = windowi0 + mIdx
    if lastPulseLocation is None:
        IPI = np.inf
    else:
        IPI = pulseIdx - lastPulseLocation
    if IPI < p['minIPI']:
        p_thresh = 1 #Cannot be a pulse
    elif IPI >= p['maxIPI']:
        p_thresh = p['pulse_p_thresh']
    else: #IPI is just right!
        p_thresh = p['pulse_consec_p_thresh']

    # Use the CNN to detect whether the pulse candidate is indeed a pulse.
    normPcand = normPcand.reshape(1, p['pulseWidth'], 1)
    pulseProb = clf.predict(normPcand)
    
    isPulse = pulseProb > p_thresh
    if isPulse:        
        lastPulseLocation = pulseIdx
        print("Pulse!")

    return isPulse, lastPulseLocation, pulseProb, pulseIdx

def getParams():
    p = {}

    p['Fs'] = 10000  # sampling frequency
    p['classifierWindow'] = int(384)  # int(0.1*Fs) is default, but depending on what to trigger this will change (TO DO)
    p['chunksize_samples'] = int(160)
    p['pulse_chSelection'] = "max" #or "mean" -- required
    p['pulse_model_tf_out'] = 'models/dm_pnp_cnn_170411_3_maxCh_weights_tfkeras.h5'# converted model output path
    
    p['pulse_p_thresh'] = 0.99
    p['pulse_consec_p_thresh'] = p['pulse_p_thresh'] # FRED: was 0.5 until 5/2/18 req, then changed to 0.99

    p['pulseWindowSize'] = p['classifierWindow']#245 ??
    p['pulseWidth'] = 165 # -- required

    p['minIPI'] = 150 #10 ms
    p['maxIPI'] = 600 #100 ms
    return p

def demo():

    # load parameters
    p = getParams()

    # Load example data
    #example_data = 'CLdata/20161213-155154bin_example_1565965_1587775.mat'
    example_data = 'CLdata/20161213-155154bin_example_4314723_4347189.mat'
    data16_example = scipy.io.loadmat(example_data)['data_example'].T.astype(np.float32)
    Nmics, Nsamples = data16_example.shape

    # Load the pulse classifier
    clf = load_pulse_classifier(p['pulse_model_tf_out'])

    pulseprobs = np.nan*np.ones(Nsamples)
    pulsePreds = []
    lastPulseLocation = None
    
    # Iterate through the data in chunks, simulating real-time processing
    for ii in range(p['pulseWindowSize'],Nsamples,p['chunksize_samples']):#range(0,Nsamples-pulseWindowSize,50):#range(Nsamples-pulseWidth):
        pulseWindow = data16_example[:,(ii-p['pulseWindowSize']):ii]
        
        isPulse, lastPulseLocation, pulseProb, pulseIdx = segPulseAux(pulseWindow, ii, lastPulseLocation, p=p, clf=clf)
        if isPulse:
            pulsePreds.append(pulseIdx)
        pulseprobs[pulseIdx] = pulseProb
    
    # Plotting
    nrows = 2
    ncols = 1

    t_s = np.arange(Nsamples) / p['Fs']

    plt.figure(figsize=(10,6))
    ax1 = plt.subplot(nrows,ncols,1)
    plt.plot(t_s,data16_example.T)
    [plt.axvline(x=t_s[pp], color='r', linestyle='--', alpha=0.5) for pp in pulsePreds]

    plt.subplot(nrows,ncols,2,sharex=ax1)
    plt.plot(t_s,pulseprobs,'ok', label='Pulse probability (model output)')
    #plt.plot(np.mean(np.abs(data16_example), axis=0)/np.max(np.mean(np.abs(data16_example), axis=0)), label='Normalized signal power (mean across mics)')
    plt.axhline(y=p['pulse_p_thresh'], color='r', linestyle='--', label='Pulse detection threshold')
    plt.ylabel('Pulse probability')
    plt.xlabel('time (s)')
    plt.legend()

    plt.tight_layout()
    plt.savefig('figures/pulse_detection_demo.png', dpi=300)
    plt.show()
    

if __name__ == '__main__':
    demo()