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

# MAT 201A Winter 2016
# 
# Homework 4
# # Yitian Shao
# 
# yitianshao@umail.ucsb.edu
# 
# Download the three HW4 audio files below. For each clip, compute the DFT of the entire file and then 
# identify the index of the most prominent peak in its magnitude spectrum. The audio files are encoded with a sampling 
# rate of 44100 Hz -- what frequencies (in Hz) correspond to the bins with the most prominent peaks?

# In[1]:


# MAT 201A Winter 2016
# Homework 4 : Download the three HW4 audio files below. For each clip, compute the DFT of the entire file and then 
# identify the index of the most prominent peak in its magnitude spectrum. The audio files are encoded with a sampling 
# rate of 44100 Hz -- what frequencies (in Hz) correspond to the bins with the most prominent peaks?
# Hint: use the fft.rfft and argmax functions. Compute one real FFT with frame size equal to the number of samples 
#     for each file.

# ======================================================================================================================
__author__ = 'Yitian Shao'
# yitianshao@umail.ucsb.edu
# Created on 02/06/2016
# The code convert a input audio signal to a image data (spectrogram)
# ======================================================================================================================
get_ipython().run_line_magic('pylab', 'inline')
from __future__ import print_function
from __future__ import division
from scipy.io import wavfile
from IPython.display import Audio

matplotlib.rcParams.update({'font.size': 16}) # Resize the font in figures

# ----------------------------------------------------------------------------------------------------------------------
# Import audios and play them (glockenspiel.wav, piano.wav, tom.wav)
sr1,glockenspiel = wavfile.read('glockenspiel.wav') 
sr2,piano = wavfile.read('piano.wav') 
sr3,tom = wavfile.read('tom.wav') 

# ----------------------------------------------------------------------------------------------------------------------
# Plot the audio signal "glockenspiel"
figure(figsize=(25,6))
plot(glockenspiel)
xlabel('Sample number')
ylabel('Magnitude')
title('glockenspiel')

print ('Signal Length = ', len(glockenspiel), ' with sampling rate = ', sr1)


# In[2]:


# Play the audio signal "glockenspiel"
Audio(glockenspiel, rate = sr1)


# In[3]:


# Plot the audio signal "piano"
figure(figsize=(25,6))
plot(piano)
xlabel('Sample number')
ylabel('Magnitude')
title('piano')

print ('Signal Length = ', len(piano), ' with sampling rate = ', sr1)


# In[4]:


# Play the audio signal "piano"
Audio(piano, rate = sr2)


# In[5]:


# Plot the audio signal "tom"
figure(figsize=(25,6))
plot(tom)
xlabel('Sample number')
ylabel('Magnitude')
title('tom')

print ('Signal Length = ', len(tom), ' with sampling rate = ', sr1)


# In[6]:


# Play the audio signal "tom"
Audio(tom, rate = sr3)


# In[7]:


# ----------------------------------------------------------------------------------------------------------------------
# Function to find the most prominent peak 
# ----------------------------------------------------------------------------------------------------------------------
def specPeak( sig, sr ):
    # Comput DFT using rfft since the input contains real number only
    # 'sig' - The input signal, in the form of an array
    # 'sr' - The sampling rate
    
    # 'Fs' - The Nyquist frequency
    Fs = 0.5 * sr

    sampleN = len(sig)

    freWave = rfft(sig)
    specMag = abs(freWave*2/sampleN)

    freq = np.linspace(0.0, Fs , 0.5 * sampleN + 1)

    figure(figsize=(25,6))
    specMag = 10 * log(specMag);
    plot(freq, specMag)
    xlim(0.0, Fs)
    xlabel('Frequency (Hz)')
    ylabel('Magnitude (DB)')
    title('Magnitude spectrum and the peak')
    
    yMin, yMax = gca().get_ylim() # Get the y-axis min and max limitation in order to plot the vertical red line

    peakInd = argmax(specMag)
    peakFreq = freq[peakInd]
     
    vlines(x = peakFreq, ymin = yMin, ymax = yMax, color = 'red')

    annotate(("Peak spectrum = %.2f Hz" % peakFreq), xy=(peakFreq, 1), xytext=(peakFreq + 200, 60), color = 'red')


# In[8]:


specPeak( glockenspiel, sr1 )


# In[9]:


specPeak( piano, sr2 )


# In[10]:


specPeak( tom, sr3 )


# <h3> In summary, the peak spectrum of "glockenspiel", "piano" and "tom" are 1323.34 Hz, 78.06 Hz and 92.57 Hz separately. </h3>