#Muhammad Hafiz Wan Rosli
#MAT240E-Final I [Time on Pitch & HMM]
#Raag Analysis
#Source: Master of Sitar- 01 Raag Alahya Bilawal (Early Morning Raag).mp3
#Ground Truth:
#Jhala: 00:00 - 00:42
#Gat:   00:42 - 18:46 (1080+46 = 1126 secs)
#Jhala: 18:46 - end (check spectral difference 1 minute before & after)
#Gat-Jhala transition: 18:45 - 18:46


#Code written in python        : https://www.python.org/
#Interactive computing ipython : http://ipython.org/
#Audio Analysis & MIR Library  : http://essentia.upf.edu/

%pylab inline
import matplotlib as plt
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = (16,4)

from essentia import *
from essentia.standard import *

#loader = MonoLoader(filename='/Users/muhammadhafiz/Documents/work/ComputationalEthno/samples/RaagAlahyaBilawal0_3minutes.mp3')
loader = MonoLoader(filename='/Users/muhammadhafiz/Music/iTunes/iTunes Media/Music/Ravi Shankar/Master Of Sitar/01 Raag Alahya Bilawal (Early Morning Raag).mp3')
audio = loader()

sampleRate = 44100
audioSnippet = audio[:60*44100] #45
sample_dur_secs = len(audioSnippet)/ float(sampleRate)
win_size=1024
hop = 512
window_start = arange(0, len(audioSnippet), hop)
w = Windowing(type = 'hann')

tonicNICM=TonicIndianArtMusic()
tonic = tonicNICM(audio)
#tonic = tonicNICM(audio[:60*44100])

tonic

pitchesConfidence = []
pitches= []

spectrum = Spectrum()
pitchyin = PitchYinFFT()

for frame in FrameGenerator(audioSnippet, frameSize = 1024, hopSize = 512):
    detectedPitch, pConfidence = pitchyin(spectrum(w(frame)))
    pitches.append(detectedPitch)
    pitchesConfidence.append(pConfidence)

max(pitches)

groundTruthAlap=42
plot(linspace(0, sample_dur_secs, len(pitches)), pitches)
hlines(tonic, 0, sample_dur_secs, color = 'g')
twinx()
vlines(groundTruthAlap, 0, max(pitches), color = 'g')
twinx()
plot(linspace(0, sample_dur_secs, len(pitchesConfidence)), pitchesConfidence, color='r')#pitch confidence plot

timeOnPitch = []
pc_count = zeros(int(max(pitches))+ 1) #for the maximum frequency
trackedPitch = pitches#[:1024]
for freq in trackedPitch:
      #print int(freq),
    pc_count[freq] +=1

argmax(pc_count), max(pc_count) #Hz, time

pc_count[284]

minFreqTimeThreshold = 0.01 * max(pc_count) #Time - 5% of max time
peaks = argwhere(pc_count > minFreqTimeThreshold)
#peaks, pc_count[peaks], int(tonic)
peaks.T

vlines(tonic, 0, max(pc_count), color='b', alpha = 0.5, lw= 3)#, linestyles = 'dashed')
hlines(minFreqTimeThreshold, 0, (int(max(pitches))+ 1), color='r')
bar(arange(int(max(pitches))+ 1), pc_count, 0.35)
xlim(100,800)
ylim(0,400)
stem(peaks, pc_count[peaks], linefmt='r--', markerfmt='r.')

from scipy.signal import argrelextrema
maxima = argrelextrema(peaks, np.greater)
minima = argrelextrema(peaks, np.less)

print(minima)

[a1, a2] = maxima

ratioToTonic=[]
for peak in peaks:
    ratio= peak/ float(int(tonic))
    ratioToTonic.append(ratio)

#peaks

peaks.T, ratioToTonic

def midi2Hz(midinote, tuning=440.0):
    return tuning * (2**((midinote - 69)/12.0))

def hz2Midi(frequency, tuning=440.0):
    return 69 + 12*(log2(frequency/ tuning))

for pitch in pitches[-100:]:
    for peak in peaks[1:]:
        if int(pitch) == peak: #round to nearest int = 1 hz apart
        #print argwhere(int(pitches) == peak)
            ###print pitch, pitches.index(pitch)
            #print pitches, argwhere(pitches)
            #plot(pitch.index(pitches), pitches)

distinctSvara = []
distinctSvaraMidi = []
distinctSvaraPos = []
for i, pitch in enumerate(pitches):
    for peak in peaks[1:-5]: #remove 0th, and n-1th peak (0, 2940)
        if int(pitch) == peak:
            ###print pitch, i
            svaraInMidi = hz2Midi(pitch)#the correlating MIDI note of the found frequencies
            distinctSvara.append(int(pitch))#rounded frequencies that are found using timeOnPitch
            distinctSvaraPos.append(i)#the positions where distinctSvara is found on the tracked pitch
            distinctSvaraMidi.append(int(svaraInMidi))#the notes in MIDI
            #plot(i, pitches)

plot(distinctSvaraPos, distinctSvara, color='b', alpha= 0.5)
#twinx()
#plot(distinctSvaraPos, distinctSvaraMidi, color='r', alpha= 0.5)
#ylim(0,500)
#xlim(200, 500)

#plot(linspace(0, sample_dur_secs, len(pitch)), pitch)
#plot(linspace(0, sample_dur_secs, len(distinctSvaraPos)), distinctSvara, 'rx')
plot(pitches)
plot(distinctSvaraPos, distinctSvara, 'r.')
hlines(tonic, 0, len(pitches), color = 'g', lw=3, alpha=0.75)
vlines(groundTruthAlap*44100/hop, 0, max(pitches), color = 'g', lw=3, alpha= 0.75)
ylim(0, peaks[-4]+100)
twinx()
plot(pitchesConfidence, 'k', alpha = 0.5)
#xlim(800, 1000)

changeInDistinctSvara = []
for i, svara in enumerate(distinctSvaraMidi):
    #print svara, i
    if distinctSvaraMidi[i] != distinctSvaraMidi[i-1]: #remove note sustains, only look for note changes based on timeonpitch notes
        changeInDistinctSvara.append(distinctSvaraMidi[i])
#        print 'true'
#    else:
#        print 'false'

#print(changeInDistinctSvara)
#plot(changeInDistinctSvara)
#plot(distinctSvaraPos, distinctSvara)

uniqueMidi = list(set(distinctSvaraMidi))
uniqueMidi

singleOctave=[]
for eachMidiNote in uniqueMidi:
    pitchClass = mod(eachMidiNote, 12)
    singleOctave.append(pitchClass)

singleOctave

singleOctave_changes = []
for eachMidiNote in changeInDistinctSvara:
    pitchClass = mod(eachMidiNote, 12)
    singleOctave_changes.append(pitchClass)

def transitionMatrix_counter(markovChain, transitionMatrix):
    for i in xrange(1, len(markovChain)):
        old_state = markovChain[i - 1]
        new_state = markovChain[i]
        transitionMatrix[old_state, new_state] += 1

transmat = np.zeros((12, 12))#12 by 12 for the chromatic scale
transitionMatrix_counter(singleOctave_changes, transmat)

for i in xrange(0, 12):
    #print transmat[i]
    transmat[i] = (transmat[i]/ sum(transmat[i]))
    for j in xrange(0, 12):
        if np.isnan(transmat[i][j]) == True:
            transmat[i][j] = 0
            #print transmat[i][j]
    print transmat[i]

plt.rcParams['figure.figsize'] = (8,8)
imshow(transmat, interpolation='none', cmap=cm.gray)
show()

imshow(transmat, interpolation='none', cmap=cm.gray)
for rownum, row in enumerate(transmat):
    for colnum, val in enumerate(row):
        text(colnum, rownum, round(transmat[rownum][colnum], 2), fontsize=10, color='white' if val < transmat.max()/2 else 'black', ha='center', va='center' )
#colorbar()
svaraNames=['Ni', 'Sa', 're', 'Re', 'ga', 'Ga', 'Ma', 'MA', 'Pa', 'dha', 'Dha', 'ni']
xticks(arange(12), svaraNames, fontsize=14)
yticks(arange(12), svaraNames, fontsize=14);
show()