#!/usr/bin/env python3
# generate_plots.py : use the matplotlib library to generate graphical plots of filter responses

import argparse
import matplotlib.pyplot as plt
import numpy as np
import scipy.signal

#--------------------------------------------------------------------------------

# This file is normally contained in a subdirectory of the actual filters, so
# this adds the parent directory to the load path in order to import the filter modules.

import sys, os
# print("__file__: ", __file__)
parent = os.path.dirname(os.path.abspath(os.path.dirname(__file__)))
# print("parent: ", parent)
sys.path.insert(0, parent)

import biquad
import hysteresis
import linear
import median
import smoothing
import statistics

################################################################
def _make_subplots(rows):
    # Create a plot or set of plots at a common figure size and resolution.
    fig, ax = plt.subplots(nrows=rows)
    fig.set_dpi(160)
    fig.set_size_inches((8,6))
    return fig, ax

################################################################
def make_figure(t_axis, y_axis, calib_axis, name, title, path, label='chirp'):
    fig, ax = _make_subplots(1)   # Create a figure containing a single set of axes.
    ax.plot(t_axis, calib_axis, label=label)
    ax.plot(t_axis, y_axis, label=name)
    ax.set_title(title)
    ax.set_xlabel("Time (seconds)")
    ax.legend()
    fig.savefig(path)
    
################################################################
if __name__ == "__main__":

    parser = argparse.ArgumentParser()
    args = parser.parse_args()

    # generate a common time axis
    delta_t = 0.1               # 10 Hz sampling
    samples = 300
    t_axis = np.arange(0, samples*delta_t, delta_t)
    
    # generate a data sequence as a 'chirp' which increases in frequency
    f0 = 0.0   # start at zero, same as the frequency plot
    f1 = 5.0   # half the Nyquist rate, same as the frequency plot
    chirp = scipy.signal.chirp(t_axis, f0, t_axis[-1], f1, phi=-90)

    # create arrays representing the response of each filter
    bandpass_f = biquad.BiquadFilter(biquad.band_pass_10_1)
    bandpass   = [bandpass_f.update(y) for y in chirp]

    bandstop_f = biquad.BiquadFilter(biquad.band_stop_10_1)
    bandstop   = [bandstop_f.update(y) for y in chirp]

    lowpass_f  = biquad.BiquadFilter(biquad.low_pass_10_1)
    lowpass    = [lowpass_f.update(y) for y in chirp]

    highpass_f = biquad.BiquadFilter(biquad.high_pass_10_1)
    highpass   = [highpass_f.update(y) for y in chirp]

    smoothing_f = smoothing.Smoothing()
    smoothing_t = [smoothing_f.update(y) for y in chirp]

    median_f    = median.MedianFilter()
    median_t    = [median_f.update(y) for y in chirp]

    # produce time plots to show the response
    sweep = "%f to %f Hz sweep" % (f0, f1)
    if False:
        make_figure(t_axis, bandpass, chirp, 'bandpass', f'Band-pass Response (1.0 Hz center, {sweep})', '../plots/bandpass-time.png')
        make_figure(t_axis, bandstop, chirp, 'bandstop', f'Band-stop Response (1.0 Hz center, {sweep})', '../plots/bandstop-time.png')    
        make_figure(t_axis, lowpass,  chirp, 'lowpass',  f'Low-pass Response (1.0 Hz center, {sweep})', '../plots/lowpass-time.png')    
        make_figure(t_axis, highpass, chirp, 'highpass', f'High-pass Response (1.0 Hz center, {sweep})', '../plots/highpass-time.png')

    if False:
        make_figure(t_axis, smoothing, chirp, 'smoothing', f'Smoothing Response ({sweep})', '../plots/smoothing-time.png')
        make_figure(t_axis, median_t,    chirp, 'median',    f'Median Filter Response ({sweep})', '../plots/median-time.png')        
    
    # create a smoother example for the median filter
    spiky = np.sin(t_axis)
    for i in range(0, len(spiky), 10):
        spiky[i] += 0.5
    median_f    = median.MedianFilter()
    median_t    = [median_f.update(y) for y in spiky]
    if False:
        make_figure(t_axis, median_t, spiky, 'median',  'Median Filter Response', '../plots/median-spiky.png', label='spiky sine')            

    # create a digital example for the debounce filter
    digital = 1.0 * (spiky > 0)
    debounce_f = hysteresis.Debounce()
    debounce_t  = [debounce_f.update(y) for y in digital]
    if False:
        make_figure(t_axis, debounce_t, digital, 'debounce',  'Debounce Filter Response', '../plots/debounce.png', label='glitchy')