Maria@4: # -*- coding: utf-8 -*-
Maria@4: """
Maria@4: Created on Mon Aug 17 13:37:00 2015
Maria@4: 
Maria@4: @author: mariapanteli
Maria@4: """
Maria@4: 
Maria@39: import numpy as np
Maria@39: import scipy.signal
Maria@4: import librosa
Maria@39: 
Maria@4: 
Maria@4: class OPMellin:
Maria@4:     def __init__(self, win2sec=8):
Maria@4:         self.y = None
Maria@4:         self.sr = None
Maria@4:         self.melspec = None
Maria@4:         self.melsr = None
Maria@4:         self.op = None
Maria@4:         self.opmellin = None
Maria@4:         self.win2sec = win2sec
Maria@39:         self.hop2sec = 0.5
Maria@4:         
Maria@39: 
Maria@39:     def load_audiofile(self, filename, sr=None, segment=True):
Maria@4:         self.y, self.sr = librosa.load(filename, sr=sr)
Maria@4:         if segment:
Maria@4:             tracklength = self.y.shape[0]/float(self.sr)
Maria@4:             startSample = 0
Maria@4:             endSample = None
Maria@4:             if tracklength > 90:
Maria@4:                 startPointSec = (tracklength/2.)-20
Maria@4:                 startSample = round(startPointSec*self.sr)
Maria@4:                 endSample = startSample+45*self.sr
Maria@4:             self.y = self.y[startSample:endSample]
Maria@4: 
Maria@4: 
Maria@4:     def mel_spectrogram(self, y=None, sr=None):
Maria@4:         if self.y is None:
Maria@4:             self.y = y
Maria@4:         if self.sr is None:
Maria@4:             self.sr = sr
Maria@4:         win1 = int(round(0.04*self.sr))
Maria@4:         hop1 = int(round(win1/8.))
Maria@39:         nfft1 = int(2**np.ceil(np.log2(win1)))
Maria@4:         nmels = 40
Maria@39:         D = np.abs(librosa.stft(self.y, n_fft=nfft1, hop_length=hop1, win_length=win1, window=scipy.signal.hamming))**2
Maria@4:         melspec = librosa.feature.melspectrogram(S=D, sr=self.sr, n_mels=nmels, fmax=8000)
Maria@4:         melsr = self.sr/float(hop1)
Maria@4:         self.melspec = melspec
Maria@4:         self.melsr = melsr
Maria@4:  
Maria@39: 
Maria@39:     def post_process_spec(self, melspec=None, log=True, medianfilt=True, sqrt=True, diff=True, subtractmean=True, halfwave=True, maxnormal=True):
Maria@4:         if self.melspec is None:
Maria@4:             self.melspec = melspec
Maria@4:         if log:
Maria@4:             self.melspec = librosa.logamplitude(self.melspec)
Maria@4:         if medianfilt:
Maria@4:             ks = int(0.1 * self.melsr) # 100ms kernel size
Maria@4:             if ks % 2 == 0: # ks must be odd
Maria@4:                 ks += 1
Maria@4:             nmels = self.melspec.shape[0]
Maria@4:             for i in range(nmels):
Maria@4:                 self.melspec[i,:] = scipy.signal.medfilt(self.melspec[i,:],kernel_size = ks)
Maria@4:         if sqrt:
Maria@4:             self.melspec = self.melspec**.5
Maria@4:         if diff:
Maria@39:             # append one frame before diff to keep number of frames the same
Maria@39:             self.melspec = np.concatenate((self.melspec,self.melspec[:,-1,None]),axis=1)  
Maria@39:             self.melspec = np.diff(self.melspec, n=1, axis=1)
Maria@4:         if subtractmean:
Maria@4:             mean = self.melspec.mean(axis=1)
Maria@4:             mean.shape = (mean.shape[0], 1)
Maria@4:             self.melspec = self.melspec - mean
Maria@4:         if halfwave:
Maria@39:             self.melspec[np.where(self.melspec < 0)] = 0
Maria@4:         if maxnormal:
Maria@4:             if self.melspec.max() != 0:  # avoid division by 0
Maria@4:                 self.melspec = self.melspec/self.melspec.max()
Maria@4: 
Maria@4: 
Maria@39:     def onset_patterns(self, melspec=None, melsr=None, center=False, bpmrange=True):
Maria@4:         if self.melspec is None:
Maria@4:             self.melspec = melspec
Maria@4:         if self.melsr is None:
Maria@4:             self.melsr = melsr
Maria@39:         win2 = int(round(self.win2sec * self.melsr))
Maria@39:         hop2 = int(round(self.hop2sec * self.melsr))
Maria@4:         nmels, nmelframes = self.melspec.shape
Maria@39:         
Maria@39:         #nfft2 = int(2**np.ceil(np.log2(win2)))
Maria@39:         nfft2 = 2048  # nfft2 does not depend on win2??
Maria@39:         melspectemp = self.melspec
Maria@39:         if ((nfft2 > win2) and (center is False)):
Maria@39:             # in librosa version < 6.0, window is padded to the size of nfft 
Maria@39:             # so if win2<nfft2 the frames returned are less than expected
Maria@39:             # solution: pad the signal by (nfft2-win2)/2 on the edges 
Maria@39:             # then frame decomposition (n_frames) matches the one expected using win2
Maria@39:             melspectemp = np.concatenate([np.zeros((nmels,int((nfft2 - win2) // 2))),self.melspec, np.zeros((nmels,int((nfft2 - win2) // 2)))],axis=1)
Maria@39:         if melspectemp.shape[1]<nfft2:
Maria@39:             # pad with zeros to have at least one 8-sec window
Maria@39:             nzeros = nfft2 - melspectemp.shape[1]
Maria@39:             melspectemp = np.concatenate([np.zeros((nmels,int(np.ceil(nzeros / 2.)))),melspectemp, np.zeros((nmels,int(np.ceil(nzeros / 2.))))],axis=1)
Maria@39:             #melspectemp = np.concatenate([melspectemp, np.zeros((nmels,nzeros))],axis=1)
Maria@39:         # nframes = int(round(np.ceil(round(nmelframes/hop2))))
Maria@39:         ff0 = np.abs(librosa.stft(y=melspectemp[0, :], win_length=win2, hop_length=hop2, n_fft=nfft2, window=scipy.signal.hamming, center=center))
Maria@39:         nframes = ff0.shape[1]
Maria@39:         # nmags, nframes = ff0.shape
Maria@39:         if bpmrange:
Maria@39:             freqresinbpm = float(self.melsr)/float(nfft2/2.)*60.
Maria@39:             minmag = int(np.floor(30./freqresinbpm))  # min tempo 30bpm
Maria@39:             maxmag = int(np.ceil(960./freqresinbpm))  # max tempo 960 bpm
Maria@39:             magsinds = range(minmag, maxmag)
Maria@39:         else:
Maria@39:             magsinds = range(0, 200)
Maria@39:         nmags = len(magsinds)
Maria@39:         fft2 = np.zeros((nmels, nmags, nframes))
Maria@39:         for i in range(nmels):
Maria@39:             fftmags = np.abs(librosa.stft(y=melspectemp[i, :], win_length=win2, hop_length=hop2, n_fft=nfft2, window=scipy.signal.hamming, center=center))
Maria@39:             fftmags = fftmags[magsinds, :]
Maria@39:             fft2[i, :, :] = fftmags
Maria@39:         op = fft2
Maria@4:         self.op = op
Maria@4: 
Maria@39: 
Maria@39:     def post_process_op(self, medianfiltOP=True):
Maria@4:         if medianfiltOP:
Maria@4:             hop2 = int(round(0.5*self.melsr))
Maria@4:             ssr = self.melsr/float(hop2)
Maria@4:             ks = int(0.5 * ssr)  # 100ms kernel size
Maria@4:             if ks % 2 == 0:  # ks must be odd
Maria@4:                 ks += 1
Maria@4:             nmels, nmags, nframes = self.op.shape
Maria@4:             for i in range(nmels):
Maria@4:                 for j in range(nframes):
Maria@4:                     #self.op[i,:,j] = scipy.signal.medfilt(self.op[i,:,j],kernel_size = ks)
Maria@39:                     self.op[i,:,j] = np.convolve(self.op[i,:,j], np.ones(ks)/ks, mode='same')
Maria@39: 
Maria@4: 
Maria@4:     def mellin_transform(self, op=None):
Maria@4:         if self.op is None:
Maria@4:             self.op = op
Maria@4:         nmels, nmags, nframes = self.op.shape
Maria@4:         #nmagsout = min(200, nmags)
Maria@4:         nmagsout = 200
Maria@39:         u_max = np.log(nmags)
Maria@39:         delta_c = np.pi / u_max
Maria@4:         c_max = nmagsout
Maria@39:         c = np.arange(delta_c, c_max, delta_c)
Maria@4:         k = range(1, nmags)
Maria@4:         exponent = 0.5 - c*1j
Maria@4: 
Maria@39:         normMat = 1./(exponent * np.sqrt(2*np.pi))
Maria@4:         normMat.shape = (normMat.shape[0], 1)
Maria@39:         normMat = np.repeat(normMat.T, nmels, axis=0)
Maria@39:         kernelMat = np.asarray([np.power(ki, exponent) for ki in k])
Maria@39:         opmellin = np.zeros((nmels, nmagsout, nframes))
Maria@4:         for i in range(nframes):
Maria@4:             self.op[:, -1, i] = 0
Maria@39:             deltaMat = - np.diff(self.op[:, :, i])
Maria@39:             mellin = np.abs(np.dot(deltaMat, kernelMat) * normMat)
Maria@4:             opmellin[:, :, i] = mellin[:, :nmagsout]
Maria@4:         self.opmellin = opmellin
Maria@4: 
Maria@4: 
Maria@4:     def post_process_opmellin(self, opmellin=None, aveFreq=False, normFrame=True):
Maria@4:         if self.opmellin is None:
Maria@4:             self.opmellin = opmellin
Maria@4:         if aveFreq:
Maria@39:             self.opmellin = np.mean(self.opmellin, axis=0)
Maria@4:         else:  # just reshape
Maria@4:             nmels, nmags, nframes = self.opmellin.shape
Maria@4:             self.opmellin = self.opmellin.reshape((nmels*nmags, nframes))
Maria@4:         if normFrame:
Maria@39:             min_opmellin = np.amin(self.opmellin, axis=0, keepdims=True)
Maria@39:             max_opmellin = np.amax(self.opmellin, axis=0, keepdims=True)
Maria@4:             denom = max_opmellin - min_opmellin
Maria@4:             denom[denom==0] = 1 # avoid division by 0 if frame is all 0s-silent
Maria@4:             self.opmellin = (self.opmellin - min_opmellin) / denom
Maria@4: 
Maria@4: 
Maria@4:     def get_opmellin(self, filename='test.wav', logfilter=False, center=True, medianfiltspec=False, medianfiltOP=False):
Maria@4:         self.load_audiofile(filename=filename)        
Maria@4:         self.mel_spectrogram()
Maria@4:         self.post_process_spec(log=False, sqrt=True, medianfilt=medianfiltspec)  # sqrt seems to work better
Maria@4:         self.onset_patterns(logfilter=logfilter, center=center)
Maria@4:         self.post_process_op(medianfiltOP=medianfiltOP)
Maria@4:         self.mellin_transform()
Maria@4:         self.post_process_opmellin()
Maria@4:         return self.opmellin
Maria@4: 
Maria@39: 
Maria@4:     def get_opmellin_from_melspec(self, melspec=[], melsr=[]):
Maria@4:         self.melspec = melspec
Maria@4:         self.melsr = melsr
Maria@4:         self.post_process_spec(log=False, sqrt=True, medianfilt=True)  # sqrt seems to work better
m@43:         self.onset_patterns(center=False)
Maria@4:         self.post_process_op(medianfiltOP=True)
Maria@4:         self.mellin_transform()
Maria@4:         self.post_process_opmellin()
Maria@4:         return self.opmellin
Maria@4: 
Maria@39: 
Maria@4: if __name__ == '__main__':
Maria@4:     op = OPMellin()
Maria@4:     op.get_opmellin()