# -*- coding: utf-8 -*-
# Nicolas Enfon - 01/04/14 - LSIS DYNI
import readfile
import numpy
from pylab import plot, show, subplot, specgram, savefig, figure
from matplotlib.pyplot import *
from scipy.io.wavfile import read,write
from random import random

#Various parameters
nfft = 128 
overlap = 64
audio_file = 'ANTARES_66496.21_full_concatenated.wav'
dat_url = 'ANTARES_66496.21_03_09_2012_23.19.41_02.20.28_T8_Q20_J60_scalo1_L1.dat'
threshold1 = 0.15 #for the detection of peaks in the trend
threshold2 = 0.2
window_signal = 0.02 #(s)
beg = 20.6 #beginning of the interesting window (s)
end = 21.3

#Reading of the input audio data
sampling_rate, data = read(audio_file)
Pxx, freqs, bins, im = specgram( data, NFFT = nfft, Fs = sampling_rate,  noverlap = overlap)
(numBins, numSpectra) = Pxx.shape
Pxx = Pxx[::-1] #to get the low freqs on the bottom

#Converting the raw data into log data
for i in xrange(len(Pxx)):#unnecessary loop?
    for j in xrange(len(Pxx[0])):
        Pxx[i][j] = 10 * numpy.log10(Pxx[i][j])

#Reading of the .dat file
dat_file_raw1 = readfile.readdat(dat_url)
dat_file1 = numpy.array(dat_file_raw1)
dat_file1 = dat_file1[60:,:] #we keep only the fisrt 60 rows (low freq)
#dat_file1 = dat_file1[::-1] #the .dat file (obtained via octave) is inverted 

dat_file_raw2 = readfile.readdat(dat_url)
dat_file2 = numpy.array(dat_file_raw2)
dat_file2 = dat_file2[30:55,:]
#dat_file2 = dat_file2[::-1]

#Binary detector
trend1 = dat_file1.sum(0)
trend2 = dat_file2.sum(0)
smoothed_trend1 = numpy.convolve(trend1, numpy.ones(35), mode='same')
smoothed_trend2 = numpy.convolve(trend2, numpy.ones(35), mode='same')
mean1 = smoothed_trend1.mean()
mean2 = smoothed_trend2.mean()
peaks1 = []
peaks2 = []
for i in xrange(len(smoothed_trend1)):
    if smoothed_trend1[i] > mean1 * (1 + threshold1):
        peaks1.append(i) #we only keep the index of the peak, for plotting
for i in xrange(len(smoothed_trend2)):
    if smoothed_trend2[i] > mean2 * (1 + threshold2):
	peaks2.append(i)

#Scaling so that it matches data's size
size1 = 400
coeff1 = 4000

peaks1 = numpy.array(peaks1)
peaks1 = peaks1 * len(data) / float(len(dat_file1[0]))
for i in xrange(len(peaks1)):
    peaks1[i] = int( round( peaks1[i] ) )#caution, numpy array keeps float type
peaks1_gate = numpy.zeros(len(data))
for i in xrange(len(peaks1)):
    deb = int(peaks1[i]-size1)
    if deb < 0:
	deb = 0
    fin = int(peaks1[i]+size1)
    if fin >= len(data):
	fin = len(data) - 1
    peaks1_gate[ deb:fin ] = coeff1

size2 = 200
coeff2 = 7000

peaks2 = numpy.array(peaks2)
peaks2 = peaks2 * len(data) / float(len(dat_file2[0]))
for i in xrange(len(peaks2)):
    peaks2[i] = int(round(peaks2[i]))
peaks2_gate = numpy.zeros(len(data))
for i in xrange(len(peaks2)):
    deb = int(peaks2[i] - size2)
    if deb < 0:
	deb = 0
    fin = int(peaks2[i] + size2)
    if fin >= len(data):
	fin = len(data) - 1
    peaks2_gate[deb:fin] = coeff2


#Conversions
beg_ratio = float(beg) / len(data) * sampling_rate 
beg_data = int( beg_ratio * len(data) )
beg_dat = int( beg_ratio * len(dat_file1[0]) )#CAUTION: normally, dat1 & dat2 are supposed to be the same size
beg_fft = int( beg_ratio * len(Pxx[0]) )
end_ratio = float(end) / len(data) * sampling_rate
end_data = int( end_ratio * len(data) )
end_dat = int( end_ratio * len(dat_file1[0]) )
end_fft = int( end_ratio * len(Pxx[0]) )

#Plotting
fig = figure(figsize=(80,20))

subplot(5,1,1)
plot(data[beg_data:end_data])
plot(peaks1_gate[beg_data:end_data],'r')
plot(peaks2_gate[beg_data:end_data],'g')
title('Signal and L1 / L2 gates')
xlim(0,end_data - beg_data + 1)

subplot(5,1,2)
plot(smoothed_trend1[beg_dat:end_dat],'r')
title('Sperm Whale detector T8 L1 - threshold =' + str(threshold1) )
xlim(0, end_dat - beg_dat + 1)

subplot(5,1,3)
plot(smoothed_trend2[beg_dat:end_dat],'g')
title('Beaked Whale detector T8 L2 - threshold =' + str(threshold2) )
xlim(0, end_dat - beg_dat + 1)

subplot(5,1,4)
imshow(Pxx[:,beg_fft:end_fft],aspect='auto')
title('Spectrogramme - NFFT =' + str(nfft) + 'overlap =' + str(overlap) )

#subplot(7,1,5)
#imshow(dat_file1[:,beg_dat:end_dat],aspect='auto')
#title('Sperm Whale part of the scalogram (60:end)')

#subplot(7,1,6)
#imshow(dat_file2[:,beg_dat:end_dat],aspect='auto')
#title('Ziphius part of the scalogram (35:60)')

subplot(5,1,5)
dat_plot = numpy.array(dat_file_raw1)
imshow(dat_plot[:,beg_dat:end_dat],aspect='auto')
title('Total scalogram')

savefig('gates.png')