#!/usr/bin/env python
# coding: utf-8

# # MAT 201A - ASSIGNMENT 5
# ### AMBIKA YADAV
# 
# ## IDENTIFYING PROMINENT PEAKS IN AUDIOS  
# 
# Three wav files have been given .The autocorrelation is computed , which is a measure of how similar is a signal to itself at different points. The points where autocorrelation is the highest correspond to the points where the signal is being repeated and in other words the distance between two peaks of autocorrelation give us the period of the original signal.
# This is how I have tried to extract the frequency of the signals by extracting the indices of the prominent peaks in the audios autocorrelation. 
# After this I have used the computation done as a part of the last assignment to extablish how accurate is the prediction of the signals frequency using the autocorrelation function.

# In[1]:


get_ipython().run_line_magic('pylab', 'inline')
rcParams['figure.figsize'] = (10, 4)
from __future__ import print_function
from __future__ import division


# In[2]:


from scipy.io import wavfile
from scipy.signal import freqs
sr_a, sample_a = wavfile.read('glockenspiel.wav')
sr_b, sample_b = wavfile.read('piano.wav')
sr_c, sample_c = wavfile.read('tom.wav')


# In[3]:


def assignment4(sr_1,sample_1):
    # FINDING THE MAGNITUDE SPECTRUM USING fft.rfft and normalising the magnitude
    sample1_f=fft.rfft(sample_1,len(sample_1))
    magspec = abs(sample1_f)/(len(sample_1)/2)
    # FINDING THE FREQUENCY BINS AS HERTZ
    length= int((len(sample_1)/2)+1)
    freq=[]
    for i in range(0,length):
        freq.append((i*sr_1)/len(sample_1))
    # COMPUTING THE TOTAL ENERGY
    T_Energy=sum(magspec)
    # HARD CODING TO FIND THE CUTOFF FREQUENCY , ASSUMING 0.95*TOTALENERGY
    CO_Bin=where(cumsum(magspec)>=0.95 *T_Energy)[0][0]
    CO_frequency=CO_Bin*sr_1/len(sample_1)
    # PLOT 1 - FREQUENCY VS MAGNITUDE
    subplot(2,2,1)
    plot(freq,magspec)
    plt.axvline(CO_frequency,color='red', linestyle='--',label='{:2f}Hz'.format(CO_frequency))
    legend()
    title('FREQUENCY VS MAGNITUDE')
    xlabel('Frequency in Hz')
    ylabel('Magnitude')
    # PLOT 2 - FREQUENCY VS MAGNITUDE IN DB
    subplot(2,2,2)
    plot(freq,20*log10(magspec/T_Energy))
    xlim(0,CO_frequency)
    title('FREQUENCY VS MAGNITUDE')
    xlabel('Frequency in Hz')
    ylabel('Magnitude in db')
    gcf().set_figheight(7)
    gcf().set_figwidth(20)
    # CALCULATING PEAKBIN AND PEAKFREQ
    peakbin=argmax(magspec)
    peakfreq=peakbin*sr_1/len(sample_1)
    print('The peak frequency is {:f}'.format(peakfreq))
    print('The cut-off frequency is {:f}'.format(CO_frequency))
    


# ### GLOCKENSPIEL

# In[4]:


lags,c, lines, line = acorr(sample_a.astype(double), maxlags= None);
grid();
gcf().set_figwidth(20)
# CONSIDERING ONLY THR POSITIVE HALF AS TH SIGNAL IS SYMMETRIC AROUND 0.
c_phalf=c[sample_a.shape[0]:] 
subplot(121)
title("Postive Half Of The AutoCorrelation")
plot(c_phalf)
xlabel('LAG')
ylabel('MAGNITUDE')
# LOOKING INTO THE SIGNAL BY LIMITING THE X VALUE
subplot(122)
plot(c_phalf[0:200])  
title("Zoomed in Positive Half")
xlabel('LAG')
ylabel('MAGNITUDE')
gcf().set_figwidth(20)

# EXTRACTING THE POINTS OF HIGHEST CORRELATION WHICH MEAN THE POINTS ANALOGUS TO THE HIGHEST FREQUENCY COMPONENT
array_of_crests=[]
crest_number=0
for i in range(1,len(sample_a)):        
    if (c_phalf[i]>c_phalf[i+1] and c_phalf[i-1]<c_phalf[i]) and (crest_number<2)  :  
        print("CREST LAG INDEX ", crest_number +1, "is", i)
        array_of_crests.append(i)
        plt.axvline(i,color='r',linestyle='--',label='{:d}'.format(i))
        legend()
        crest_number =crest_number+1
        # EXTRACTING THE CORRESPONDING LAG VALUES 
        if crest_number==2:      
            break
period=array_of_crests[1]-array_of_crests[0]
fundamental_frequency=sr_a/period            
print("The fundamental frequency extracted from Autocorrelation is",fundamental_frequency, "Hz")


# In[5]:


assignment4(sr_a,sample_a)


# ### TOM

# In[6]:


lags,c, lines, line = acorr(sample_c.astype(double), maxlags= None);
grid();
gcf().set_figwidth(20)
# CONSIDERING ONLY THR POSITIVE HALF AS TH SIGNAL IS SYMMETRIC AROUND 0.
c_phalf=c[sample_c.shape[0]:] 
subplot(121)
title("Postive Half Of The AutoCorrelation")
plot(c_phalf)
xlabel('LAG')
ylabel('MAGNITUDE')
# LOOKING INTO THE SIGNAL BY LIMITING THE X VALUE
subplot(122)
plot(c_phalf[0:1200])  
title("Zoomed in Positive Half")
xlabel('LAG')
ylabel('MAGNITUDE')
gcf().set_figwidth(20)

# EXTRACTING THE POINTS OF HIGHEST CORRELATION WHICH MEAN THE POINTS ANALOGUS TO THE HIGHEST FREQUENCY COMPONENT
array_of_crests=[]
crest_number=0
for i in range(1,len(sample_c)):        
    if (c_phalf[i]>c_phalf[i+1] and c_phalf[i-1]<c_phalf[i]) and (crest_number<2)  :  
        print("CREST LAG INDEX ", crest_number +1, "is", i)
        array_of_crests.append(i)
        plt.axvline(i,color='r',linestyle='--',label='{:d}'.format(i))
        legend()
        crest_number =crest_number+1
        # EXTRACTING THE CORRESPONDING LAG VALUES 
        if crest_number==2:      
            break
period=array_of_crests[1]-array_of_crests[0]
fundamental_frequency=sr_c/period             
print("The fundamental frequency extracted from Autocorrelation is",fundamental_frequency, "Hz")


# In[7]:


assignment4(sr_c,sample_c)


# ### PIANO

# In[8]:


lags,c, lines, line = acorr(sample_b.astype(double), maxlags= None);
grid();
gcf().set_figwidth(20)
# CONSIDERING ONLY THR POSITIVE HALF AS TH SIGNAL IS SYMMETRIC AROUND 0.
c_phalf=c[sample_b.shape[0]:] 
subplot(121)
title("Postive Half Of The AutoCorrelation")
plot(c_phalf)
xlabel('LAG')
ylabel('MAGNITUDE')
# LOOKING INTO THE SIGNAL BY LIMITING THE X VALUE
subplot(122)
plot(c_phalf[0:1500])  
title("Zoomed in Positive Half")
xlabel('LAG')
ylabel('MAGNITUDE')
gcf().set_figwidth(20)

# EXTRACTING THE POINTS OF HIGHEST CORRELATION WHICH MEAN THE POINTS ANALOGUS TO THE HIGHEST FREQUENCY COMPONENT
array_of_crests=[]
crest_number=0

# HARD CODING THE INDEX HERE BECAUSE THERE ARE A LOT OF PEAKS HERE ( TOO MANY HARMONICS )
for i in range(1,len(sample_b)):        
    if (c_phalf[i]>c_phalf[i+1] and c_phalf[i-1]<c_phalf[i]) and (crest_number<16)  :  
        if crest_number == 6 or crest_number == 13 :
            print("CREST LAG INDEX ", crest_number +1, "is", i)
            array_of_crests.append(i)
            plt.axvline(i,color='r',linestyle='--',label='{:d}'.format(i))
            legend()
        crest_number =crest_number+1
        # EXTRACTING THE CORRESPONDING LAG VALUES 
        if crest_number==16:      
            break
period=array_of_crests[1]-array_of_crests[0]
fundamental_frequency=sr_b/period             
print("The fundamental frequency extracted from Autocorrelation is",fundamental_frequency, "Hz")


# In[9]:


assignment4(sr_b,sample_b)