import numpy as np, matplotlib.pyplot as plt, PhyPraKit as ppk
from scipy import interpolate, signal


fname = "signal.csv"

hlines, data = ppk.readCSV(fname, nlhead=1)
print (hlines)
print (" --> Anzahl Spalten", data.shape[1])
print (" --> Datenzeilen", data.shape[1])


t = data[0]
a = data[1]

glatt = int(input("Glättungsparameter:"))
a_smooth = ppk.meanFilter(a, glatt)

#bilde Differenz
differenz = a - a_smooth

# make plots
fig=plt.figure(1, figsize=(7.5, 7.5))
fig.suptitle('Signal mit Glättung und Differenz')
fig.subplots_adjust(left=0.14, bottom=0.1, right=0.97, top=0.93,
                wspace=None, hspace=.75)#
ax1=fig.add_subplot(3, 1, 1)
ax1.plot(t, a)
ax1.set_title('Rohdaten')
ax1.set_xlabel('Zeit  (ms)', size='large')
ax1.set_ylabel('Amplitude (mV)', size='large')
ax1.grid()

ax2=fig.add_subplot(3, 1, 2)
ax2.plot(t, a_smooth)
ax1.set_title('Geglättetes Signal')
ax2.set_xlabel('Zeit (ms)', size='large')
ax2.set_ylabel(u'geglättete Amplitude (mV)', size='large')
ax2.grid()

ax3=fig.add_subplot(3, 1, 3)
ax3.plot(t, differenz)
ax3.set_title('Differenz')
ax3.set_xlabel('Zeit (ms)', size='large')
ax3.set_ylabel(u'Differenz der Amplitude (mV)', size='large')
ax3.grid()

plt.show()

#
ac_a = ppk.autocorrelate(a_smooth)
ac_t = t

# find maxima and minima using convolution peak finder
width= 5
pidxac =  ppk.convolutionPeakfinder(ac_a, width, th=0.66) #ursprünglich 0.66 Anzahl kreuze (am ende überdeckt)
didxac =  ppk.convolutionPeakfinder(-ac_a, width, th=0.66) #ursprünglich 0.66
if len(pidxac) > 1:
    print(" --> %i auto-correlation peaks found"%(len(pidxac)))
    pidxac[0]=0 # first peak is at 0 by construction
    ac_tp, ac_ap= np.array(ac_t[pidxac]), np.array(ac_a[pidxac])
    ac_td, ac_ad= np.array(ac_t[didxac]), np.array(ac_a[didxac])
else:
    print("*!!* not enough correlation peaks found")

# make plots
# 1. signal
fig = plt.figure(1, figsize=(7.5, 9.))
fig.suptitle('Periodendauer/ Verteilung Zeitabstände')
fig.subplots_adjust(left=0.1, bottom=0.1, right=0.98, top=0.98,
                  wspace=0.3, hspace=0.5)

ax1 = fig.add_subplot(3, 1, 1)
ax1.plot(t, a_smooth)
ax1.set_title('geglättetes Signal')
ax1.set_xlabel('Zeit (ms)', size='large')
ax1.set_ylabel('$amplitude$ (a.u.)', size='large')
ax1.grid()

# 2. auto-correlation
ax2=fig.add_subplot(3,1,2)
ax2.plot(ac_tp, ac_ap, 'bx', alpha=0.9, label='peaks')
ax2.plot(ac_td, ac_ad, 'rx', alpha=0.9, label='dips')
ax2.plot(ac_t, ac_a)
ax2.set_title('Korrelation')
ax2.set_xlabel('Zeit(ms) ', size='large')
ax2.set_ylabel('$autocorrelation$', size='large')
ax2.legend(loc='best', numpoints=1, prop={'size':10})
#  ax2.set_yscale('log')
ax2.grid()

# 3. analysis of auto-correlation function
ax3 = fig.add_subplot(3, 1, 3)
ac_dtp = ac_tp[1:] - ac_tp[:-1]
ac_dtd = ac_td[1:] - ac_td[:-1]
bins=np.linspace(min(min(ac_dtp),min(ac_dtd)), max(max(ac_dtp), max(ac_dtd)), 100)
bc, be, _ = ax3.hist([ac_dtp, ac_dtd], bins, stacked = True,
                     color=['b','r'], label=['peaks','dips'], alpha=0.5)
ax3.set_title('')
ax3.set_xlabel('Zeitdifferenz von Maxima und Minima$ (ms)', size='large')
ax3.set_ylabel('Frequenz', size='large')
ax3.legend(loc='best', numpoints=1, prop={'size':10})
ax3.grid()

plt.show()

# analyis of histogram
m_dtp, s_dtp, sm_dtp = ppk.histstat(bc[0], be, pr=False)
m_dtd, s_dtd, sm_dtd = ppk.histstat(bc[1], be, pr=False)

print("Periodendauer entspricht:")
print(" --> Zeitdifferenz Auto-Correlation: (%.5g +/- %.2g) ms"%(m_dtp, sm_dtp))
print(" -->                                 (%.5g +/- %.2g) ms"%(m_dtd, sm_dtd))

#Fourierdarstellung
print("** Fourier Spectrum")
freq, amp = ppk.FourierSpectrum(t, a_smooth, fmax=20)
# freq, amp = ppk.Fourier_fft(t, a)  # use fast algorithm
frequency = freq[np.where(amp==max(amp))]
print(" --> Frequenz mit max. Amplitude: ", frequency)

# make  plots
fig=plt.figure(1, figsize=(7.5, 7.5))
fig.suptitle('Fourier-Transformerte')
fig.subplots_adjust(left=0.14, bottom=0.1, right=0.97, top=0.93,
                wspace=None, hspace=.25)#
ax1=fig.add_subplot(2, 1, 1)
ax1.plot(t, a_smooth)
ax1.set_xlabel('Zeit in ms', size='large')
ax1.set_ylabel('Amplitude mV', size='large')
ax1.grid()

ax2=fig.add_subplot(2,1,2)
ax2.plot(freq, amp, 'b')
ax2.set_xlabel('$Frequenz$ $f$ (kHz)', size='large')
ax2.set_ylabel('$Amplitude$', size='large')
ax2.set_yscale('log')
ax2.grid()

plt.show()