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

from signalprocess import SignalProcess

# Read data from the video processing pipeline
sig = SignalProcess(path="out.tsv")
# Uncomment to see fft plot
# sig.fft()

# Applies the butter filter from the Signal Process class
data = sig.butter_filter()


plt.plot(sig.times, data)
plt.plot(sig.times, sig.ave_counts)
# Uncomment to see raw waveform
# plt.plot(sig.times,sig.counts)

# Peak finding algo provided by scipy package
peaks, descript = signal.find_peaks((-1) * data, prominence=0.2, distance=30)

# Plot the calculated peaks
plt.scatter(
    sig.times[peaks],
    data[peaks],
)

# Array of peak to peak distances
distances = []
# Array of the peak heights as calculated in the plot
heights = []
for i in range(len(peaks)):
    heights.append(data[peaks[i]])

for i in range(1, len(peaks) - 1):
    distance = (peaks[i] - peaks[i - 1]) / 30.0
    distances.append(distance)

# Basic stats
print("distances: median={}, std={}".format(np.median(distances), np.std(distances)))
print("heights: median={}, std={}".format(np.median(heights), np.std(heights)))

# Standard deviation rejection of outliers
d = np.abs(heights - np.median(heights))
mdev = np.median(d)
s = d / mdev if mdev else 0.0
peaks = peaks[s < 3.0]

# plot the non-rejected data on top of the old peak data
plt.scatter(
    sig.times[peaks],
    data[peaks],
)

# Draw that bitch
plt.show()