<jupyter_start><jupyter_text>Test Your Algorithm Instructions1. From the **Pulse Rate Algorithm** Notebook you can do one of the following:   - Copy over all the **Code** section to the following Code block.   - Download as a Python (`.py`) and copy the code to the following Code block.2. In the bottom right, click the Test Run button.  Didn't PassIf your code didn't pass the test, go back to the previous Concept or to your local setup and continue iterating on your algorithm and try to bring your training error down before testing again. PassIf your code passes the test, complete the following! You **must** include a screenshot of your code and the Test being **Passed**. Here is what the starter filler code looks like when the test is run and should be similar. A passed test will include in the notebook a green outline plus a box with **Test passed:** and in the Results bar at the bottom the progress bar will be at 100% plus a checkmark with **All cells passed**.![Example](example.png)1. Take a screenshot of your code passing the test, make sure it is in the format `.png`. If not a `.png` image, you will have to edit the Markdown render the image after Step 3. Here is an example of what the `passed.png` would look like 2. Upload the screenshot to the same folder or directory as this jupyter notebook.3. Rename the screenshot to `passed.png` and it should show up below.![Passed](passed.png)4. Download this jupyter notebook as a `.pdf` file. 5. Continue to Part 2 of the Project.<jupyter_code># replace the code below with your pulse rate algorithm.
# import numpy as np
# import scipy as sp
# import scipy.io

# def RunPulseRateAlgorithm(data_fl, ref_fl):
#     ref = sp.io.loadmat(ref_fl)['BPM0'].reshape(-1)
#     sample_idx = np.array([0, 250, 500, 750, 1000, 1250, 1500, 1750, 2000, 2250, 2500, 2750, 3000, 3250, 3500, 3750, 4000, 4250, 4500, 4750, 5000, 5250, 5500, 5750, 6000, 6250, 6500, 6750, 7000, 7250, 7500, 7750, 8000, 8250, 8500, 8750, 9000, 9250, 9500, 9750, 10000, 10250, 10500, 10750, 11000, 11250, 11500, 11750, 12000, 12250, 12500, 12750, 13000, 13250, 13500, 13750, 14000, 14250, 14500, 14750, 15000, 15250, 15500, 15750, 16000, 16250, 16500, 16750, 17000, 17250, 17500, 17750, 18000, 18250, 18500, 18750, 19000, 19250, 19500, 19750, 20000, 20250, 20500, 20750, 21000, 21250, 21500, 21750, 22000, 22250, 22500, 22750, 23000, 23250, 23500, 23750, 24000, 24250, 24500, 24750, 25000, 25250, 25500, 25750, 26000, 26250, 26500, 26750, 27000, 27250, 27500, 27750, 28000, 28250, 28500, 28750, 29000, 29250, 29500, 29750, 30000, 30250, 30500, 30750, 31000, 31250, 31500, 31750, 32000, 32250, 32500, 32750, 33000, 33250, 33500, 33750, 34000, 34250, 34500, 34750, 35000, 35250, 35500, 35750, 36000, 36250, 36500, 36750, 37000, 37250, 37500, 37750, 38000])
#     pr_conf = np.array([0.021337153215872807, 0.024064924996852304, 0.02463020949150536, 0.01850781929621976, 0.01798823606401656, 0.015214068166993795, 0.022201019635165256, 0.023496350505833742, 0.02284078817745135, 0.020909684967874518, 0.01563895587931265, 0.018860372259259035, 0.01962636016755157, 0.020945427799021643, 0.020210754809204513, 0.02060068846306182, 0.014944962628767091, 0.016877531273948632, 0.016231267453579862, 0.01722237292934081, 0.02254199324355656, 0.021656123596103808, 0.021185002253523294, 0.014420016080975444, 0.010299535624137768, 0.01038923070423021, 0.02343286315153606, 0.02860563654164102, 0.029497906615316982, 0.02388068014629748, 0.016667878703836216, 0.020242117638493318, 0.018640855806529574, 0.021831371828396888, 0.018254153939637585, 0.018829649283083266, 0.016142159003363295, 0.011337153547969549, 0.012296850438562329, 0.014694922793814558, 0.022694057873903946, 0.022166526448369032, 0.02289066235147056, 0.014872491261353479, 0.013945491280205995, 0.013469187903620325, 0.02133081232636454, 0.025958326164926982, 0.027304335603705383, 0.027709186267831084, 0.017495162891028045, 0.013364331102384868, 0.013070030370849745, 0.02031604844340968, 0.019293254059602423, 0.021605291432033193, 0.021023889750224935, 0.017321384494621166, 0.015149938381169913, 0.018838962183624786, 0.025925988532678426, 0.02903979368896086, 0.02868298533314268, 0.023057954230949956, 0.015228263334995955, 0.014019147093563533, 0.017022121236989998, 0.019831336765500744, 0.02061937636876943, 0.021636620552160615, 0.017595357423085337, 0.015051476624567706, 0.015426293389341369, 0.022688369156979554, 0.024392854120312033, 0.023280894840993358, 0.022340994234100887, 0.01732588839304203, 0.017110754222286795, 0.015151854898099672, 0.022079944297238606, 0.021111087340285863, 0.021399886991470905, 0.019877249666323385, 0.012442094909825943, 0.013659691010914637, 0.015758158418059307, 0.023152968433707524, 0.024554383607661916, 0.025609371650330826, 0.01648589212834377, 0.010967020013546163, 0.014360217112495348, 0.021658825752178697, 0.02581615068336972, 0.024365834821670977, 0.021359180583158867, 0.014597287141889194, 0.013871293617294026, 0.014950177066137254, 0.021663566020926135, 0.022148074222132084, 0.025118713242048184, 0.023559799018199337, 0.020626412198778414, 0.018086950917227427, 0.01461833235997816, 0.02308451506996894, 0.020805442648792628, 0.024085578390203723, 0.020377343948877333, 0.013399120284383535, 0.012353069061434472, 0.017042942709335387, 0.02318520616268311, 0.022847671258043923, 0.02409832825334578, 0.01784425180604167, 0.017275343657836877, 0.016994825583948103, 0.02317302366180944, 0.02539252664198577, 0.026498867456691444, 0.025584933217795643, 0.020291893799630097, 0.016119578944077263, 0.01561926846014978, 0.023966078931255276, 0.026596112265947987, 0.027874868480540355, 0.02562343627211549, 0.016491737162859402, 0.013560977931000669, 0.014156112036075757, 0.02081362361637795, 0.02109895800587278, 0.024756849920808358, 0.023707623471236807, 0.01828992367871122, 0.015203248827854104, 0.016907540501959875, 0.02056417746512964, 0.018946284075005378, 0.022342383921091453, 0.017858937318584536, 0.012816042168425438, 0.01289760679993461, 0.016218133075097767, 0.02006455714343904, 0.022335287222511373, 0.027133097359541277, 0.01687396198324497, 0.009829891484050932])
#     pr_est = np.array([49.93333333333333, 49.96666666666667, 49.93333333333333, 49.93333333333333, 49.983333333333334, 49.93333333333333, 50.11666666666667, 49.93333333333333, 49.95, 49.93333333333333, 49.93333333333333, 50.016666666666666, 49.93333333333333, 49.916666666666664, 49.9, 49.93333333333333, 49.93333333333333, 50.0, 50.0, 49.916666666666664, 49.916666666666664, 49.93333333333333, 49.93333333333333, 49.93333333333333, 24.966666666666665, 24.933333333333334, 50.083333333333336, 50.083333333333336, 49.983333333333334, 49.95, 49.93333333333333, 50.0, 49.93333333333333, 50.11666666666667, 24.983333333333334, 50.13333333333333, 49.93333333333333, 32.78333333333333, 48.46666666666667, 50.11666666666667, 49.93333333333333, 49.93333333333333, 49.93333333333333, 49.93333333333333, 49.96666666666667, 49.95, 49.93333333333333, 49.95, 49.983333333333334, 49.96666666666667, 49.95, 49.96666666666667, 48.61666666666667, 50.11666666666667, 50.06666666666666, 49.96666666666667, 49.93333333333333, 49.93333333333333, 49.93333333333333, 49.95, 49.95, 49.983333333333334, 49.95, 49.93333333333333, 49.93333333333333, 50.0, 50.11666666666667, 50.1, 49.95, 49.95, 49.93333333333333, 49.96666666666667, 50.06666666666666, 50.1, 50.1, 49.96666666666667, 49.95, 49.95, 49.96666666666667, 49.916666666666664, 50.1, 50.06666666666666, 49.96666666666667, 49.93333333333333, 49.95, 50.0, 49.916666666666664, 49.93333333333333, 50.06666666666666, 49.96666666666667, 50.1, 50.083333333333336, 49.96666666666667, 49.916666666666664, 49.93333333333333, 49.93333333333333, 49.93333333333333, 49.93333333333333, 49.983333333333334, 50.11666666666667, 50.1, 49.93333333333333, 49.95, 49.93333333333333, 49.916666666666664, 49.96666666666667, 49.916666666666664, 50.11666666666667, 50.083333333333336, 50.0, 49.93333333333333, 49.93333333333333, 49.96666666666667, 49.93333333333333, 49.93333333333333, 50.0, 49.95, 49.93333333333333, 49.96666666666667, 49.95, 50.1, 49.9, 49.96666666666667, 49.93333333333333, 49.93333333333333, 49.96666666666667, 50.083333333333336, 50.1, 50.06666666666666, 49.983333333333334, 49.93333333333333, 49.93333333333333, 50.03333333333333, 50.083333333333336, 50.083333333333336, 49.96666666666667, 49.95, 49.93333333333333, 49.96666666666667, 49.983333333333334, 50.1, 50.1, 49.9, 49.93333333333333, 49.93333333333333, 49.983333333333334, 50.0, 49.916666666666664, 49.916666666666664, 49.96666666666667, 50.0, 49.96666666666667, 53.233333333333334])
#     ref_idx = np.cumsum(np.ones(len(ref)) * 125 * 2) - 125 * 2
#     return pr_est[:len(ref)] * 60 - ref, pr_conf[:len(ref)]

import glob

import numpy as np
import scipy as sp
import scipy.io
import matplotlib.pyplot as plt
import scipy.signal


def LoadTroikaDataset():
    """
    Retrieve the .mat filenames for the troika dataset.

    Review the README in ./datasets/troika/ to understand the organization of the .mat files.

    Returns:
        data_fls: Names of the .mat files that contain signal data
        ref_fls: Names of the .mat files that contain reference data
        <data_fls> and <ref_fls> are ordered correspondingly, so that ref_fls[5] is the 
            reference data for data_fls[5], etc...
    """
    data_dir = "./datasets/troika/training_data"
    data_fls = sorted(glob.glob(data_dir + "/DATA_*.mat"))
    ref_fls = sorted(glob.glob(data_dir + "/REF_*.mat"))
    return data_fls, ref_fls

def LoadTroikaDataFile(data_fl):
    """
    Loads and extracts signals from a troika data file.

    Usage:
        data_fls, ref_fls = LoadTroikaDataset()
        ppg, accx, accy, accz = LoadTroikaDataFile(data_fls[0])

    Args:
        data_fl: (str) filepath to a troika .mat file.

    Returns:
        numpy arrays for ppg, accx, accy, accz signals.
    """
    data = sp.io.loadmat(data_fl)['sig']
    return data[2:]

def AggregateErrorMetric(pr_errors, confidence_est):
    """
    Computes an aggregate error metric based on confidence estimates.

    Computes the MAE at 90% availability. 

    Args:
        pr_errors: a numpy array of errors between pulse rate estimates and corresponding 
            reference heart rates.
        confidence_est: a numpy array of confidence estimates for each pulse rate
            error.

    Returns:
        the MAE at 90% availability
    """
    # Higher confidence means a better estimate. The best 90% of the estimates
    #    are above the 10th percentile confidence.
    percentile90_confidence = np.percentile(confidence_est, 10)

    # Find the errors of the best pulse rate estimates
    best_estimates = pr_errors[confidence_est >= percentile90_confidence]

    # Return the mean absolute error
    return np.mean(np.abs(best_estimates))

def Evaluate():
    """
    Top-level function evaluation function.

    Runs the pulse rate algorithm on the Troika dataset and returns an aggregate error metric.

    Returns:
        Pulse rate error on the Troika dataset. See AggregateErrorMetric.
    """
    # Retrieve dataset files
    data_fls, ref_fls = LoadTroikaDataset()
    errs, confs = [], []
    for data_fl, ref_fl in zip(data_fls, ref_fls):
        # Run the pulse rate algorithm on each trial in the dataset
        errors, confidence = RunPulseRateAlgorithm(data_fl, ref_fl)
        errs.append(errors)
        confs.append(confidence)
        # Compute aggregate error metric
    errs = np.hstack(errs)
    confs = np.hstack(confs)
    return AggregateErrorMetric(errs, confs)

def RunPulseRateAlgorithm(data_fl, ref_fl):
    # Load data using LoadTroikaDataFile
    ppg, accx, accy, accz = LoadTroikaDataFile(data_fl)
    # Load the BMP reference data
    bpm_ref = np.array(sp.io.loadmat(ref_fl)['BPM0'])
    
    # Compute pulse rate estimates and estimation confidence.
    # Extracting the dominant frequencies of signals
    f_ppg, cf = freq_analyzer(ppg)
    f_ax, _ = freq_analyzer(accx)
    f_ay, _ = freq_analyzer(accy)
    f_az, _ = freq_analyzer(accz)
    # As demostrated in readme file, if the first strongest frequency of PPG is 
    # same as the acc strongest freq, the secon strongest freq is picked
    bpm=[]
    conf=[]
    for i in range(f_ppg.shape[0]):
#         bpm.append(f_ppg[i,0])
#         conf.append(cf[i,0])
        k=0
        for j, f in enumerate(f_ppg[i,:]):
            if f in (f_ax[i,k], f_ay[i,k], f_az[i,k]):
                k+=1
                pass
            else:
                bpm.append(f)
                conf.append(cf[i,j])
                break
            if k==3:
                bpm.append(f_ppg[i,0])
                conf.append(cf[i,0])

    errors = np.abs(bpm_ref.squeeze() - np.array(bpm))
    confidence = conf
    # Return per-estimate mean absolute error and confidence as a 2-tuple of numpy arrays.
#     errors, confconfidence = confidence = np.ones(100), np.ones(100)  # Dummy placeholders. Remove
    return errors, confidence


##
# Helper functions
fs=125
# freq_analyzer use the specgram functions to explore the fft of the signal in window of 8s
# and overlap of 6s. The strongest freqs are then seleccted using sp.signal.find_peaks method
# which are later sorted accorging to their corresponding magneitude. Finally three top strong freqs
# are extracted.
def freq_analyzer(signal, fs=fs):
    dominant_fr=[]
    confidence=[]
    spec, fr, t, _ = plt.specgram(bpf(signal, fs), NFFT=fs*8 , Fs=fs, noverlap=6*fs)
    spec_energy = np.sum(spec**2, axis=0)

    for i in range(spec.shape[1]):
        index = sp.signal.find_peaks(spec[:40,i])[0]
        idx2 = sorted(index, key=lambda x: spec[x,i], reverse=True)
        selected_fr = fr[idx2]*60
        dominant_fr.append(selected_fr[:3])
        
        # uisng signal energy critiria to calculate the confidence
        # The proportion of selected freq energy within 1Hz to the total signal energy
        spec_energy = spec[:,i]**2
        energy_ratio=[]
        for j in idx2[:3]:
            # get signals 0.5 Hz lower and higher than the selected peak freq:
            lower_bound = j-4 if j>7 else 0
            upper_bound = j+5
            energy_ratio.append(np.sum(spec_energy[lower_bound:upper_bound])/np.sum(spec_energy))
        confidence.append(energy_ratio)
        
    return np.array(dominant_fr), np.array(confidence)

# Bandpass filter function to filter signal freq between 40 bpm and 240 bpm
def bpf(signal, fs):
    b, a = sp.signal.butter(3, (40/60, 240/60), btype='bandpass', fs=fs)
    return sp.signal.filtfilt(b, a, signal)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10027341079972658, 'fraction_numerical': 0.38803827751196174, 'mean_word_length': 4.386966126656849, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 114, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '482544', 'license': 'MIT', 'n_tokens_mistral': 9322, 'n_tokens_neox': 5540, 'n_words': 1330, 'path': 'unit_test.ipynb', 'repo_name': 'ashad12/PhysiologicalSignalsAI-Udacity', 'size': '14630'}