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@4: import librosa
Maria@4: import scipy.signal
Maria@4: # import matplotlib.pyplot as plt
Maria@4: import numpy
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@4:         
Maria@4:     def load_audiofile(self, filename='test.wav', sr=None, segment=True):
Maria@4:         self.y, self.sr = librosa.load(filename, sr=sr)
Maria@4:         #self.y, self.sr = librosa.load(filename) # default sr=22050
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, filename='1_21.wav'):
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@4:         nfft1 = int(2**numpy.ceil(numpy.log2(win1)))
Maria@4:         nmels = 40
Maria@4:         # y, sr = librosa.load(filename, sr=fs)        
Maria@4:         # melspec = librosa.feature.melspectrogram(y=y, sr=sr, n_fft=nfft, win_length=win1, window=scipy.signal.hamming, hop_length=hop1, n_mels=nmels)
Maria@4:         # melspec = librosa.feature.melspectrogram(y=self.y, sr=self.sr, n_fft=nfft, hop_length=hop1, n_mels=nmels)
Maria@4:         D = numpy.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)
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:         # librosa.display.specshow(melspec)
Maria@4:         # return melspec, melsr
Maria@4:         self.melspec = melspec
Maria@4:         self.melsr = melsr
Maria@4:  
Maria@4:     # post process spectrogram
Maria@4:     def post_process_spec(self, melspec=None, log=True, medianfilt=False, sqrt=False, smooth=False, 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 smooth:
Maria@4:             step = 50.0
Maria@4:             nmels = self.melspec.shape[0]
Maria@4:             mm = self.melspec
Maria@4:             for i in range(nmels):
Maria@4:                 mm[i,:] = numpy.convolve(mm[i,:], numpy.ones(step)/step, mode='same')
Maria@4:             self.melspec = mm
Maria@4:         if diff:
Maria@4:             self.melspec = numpy.concatenate((self.melspec,self.melspec[:,-1,None]),axis=1)  # append one frame before diff to keep number of frames the same
Maria@4:             self.melspec = numpy.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@4:             self.melspec[numpy.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:         # return melspec
Maria@4: 
Maria@4: 
Maria@4:     def onset_patterns(self, melspec=None, melsr=None, fft=True, logfilter=False, center=True, 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@4:         win2 = int(round(self.win2sec*self.melsr))
Maria@4:         hop2 = int(round(0.5*self.melsr))
Maria@4:         nmels, nmelframes = self.melspec.shape
Maria@4:         if fft:  # periodicity via fft
Maria@4:             #nfft2 = int(2**numpy.ceil(numpy.log2(win2)))
Maria@4:             nfft2 = 2048  # nfft2 does not depend on win2??
Maria@4:             melspectemp = self.melspec
Maria@4:             if ((nfft2 > win2) and (center is False)):
Maria@4:                 # pad the signal by nfft2-win2 so that frame decomposition (n_frames) 
Maria@4:                 # is still based on win2 and not nfft2
Maria@4:                 melspectemp = numpy.concatenate([numpy.zeros((nmels,int((nfft2 - win2) // 2))),self.melspec, numpy.zeros((nmels,int((nfft2 - win2) // 2)))],axis=1)
Maria@4:             if melspectemp.shape[1]<nfft2:
Maria@4:                 # pad with zeros to have at least one 8-sec window
Maria@4:                 nzeros = nfft2 - melspectemp.shape[1]
Maria@4:                 melspectemp = numpy.concatenate([numpy.zeros((nmels,int(numpy.ceil(nzeros / 2.)))),melspectemp, numpy.zeros((nmels,int(numpy.ceil(nzeros / 2.))))],axis=1)
Maria@4:                 #melspectemp = numpy.concatenate([melspectemp, numpy.zeros((nmels,nzeros))],axis=1)
Maria@4:             # nframes = int(round(numpy.ceil(round(nmelframes/hop2))))
Maria@4:             ff0 = numpy.abs(librosa.stft(y=melspectemp[0, :], win_length=win2, hop_length=hop2, n_fft=nfft2, window=scipy.signal.hamming, center=center))
Maria@4:             nframes = ff0.shape[1]
Maria@4:             # nmags, nframes = ff0.shape
Maria@4:             if logfilter:  # log filterbank
Maria@4:                 nbins = 25
Maria@4:                 binsperoct = 5
Maria@4:                 fft2 = numpy.zeros((nmels, nbins, nframes))     
Maria@4:                 logfb = librosa.filters.logfrequency(self.melsr, nfft2, n_bins=nbins, bins_per_octave=binsperoct, fmin=0.5)
Maria@4:                 # logfb = librosa.filters.logfrequency(self.melsr, nfft2, n_bins=25, bins_per_octave=5, fmin=0.5)
Maria@4:                 for i in range(nmels):
Maria@4:                     fftmags = numpy.abs(librosa.stft(y=melspectemp[i, :], win_length=win2, hop_length=hop2, n_fft=nfft2, window=scipy.signal.hamming, center=center))
Maria@4:                     fftmags = numpy.dot(logfb, fftmags)  # log
Maria@4:                     fft2[i, :, :] = fftmags
Maria@4:             else:
Maria@4:                 # nmags = min(round(nfft2/2), 200)
Maria@4:                 if bpmrange:
Maria@4:                     # freqresinbpm = float(self.melsr)/float(nfft2)*60.
Maria@4:                     freqresinbpm = float(self.melsr)/float(nfft2/2.)*60.
Maria@4:                     minmag = int(numpy.floor(30./freqresinbpm))  # min tempo 30bpm
Maria@4:                     maxmag = int(numpy.ceil(960./freqresinbpm))  # max tempo 960 bpm
Maria@4:                     magsinds = range(minmag, maxmag)
Maria@4:                 else:
Maria@4:                     magsinds = range(0, 200)
Maria@4:                 nmags = len(magsinds)
Maria@4:                 fft2 = numpy.zeros((nmels, nmags, nframes))
Maria@4:                 for i in range(nmels):
Maria@4:                     fftmags = numpy.abs(librosa.stft(y=melspectemp[i, :], win_length=win2, hop_length=hop2, n_fft=nfft2, window=scipy.signal.hamming, center=center))
Maria@4:                     fftmags = fftmags[magsinds, :]
Maria@4:                     fft2[i, :, :] = fftmags
Maria@4:             op = fft2
Maria@4:         else:  # periodicity via autocorrelation
Maria@4:             nmags = min(win2, 1000)
Maria@4:             if nmelframes < win2:
Maria@4:                 nframes = 1
Maria@4:             else:
Maria@4:                 nframes = int(nmelframes/float(win2))
Maria@4:             autocor = numpy.zeros((nmels, nmags, nframes))
Maria@4:             for i in range(nmels):
Maria@4:                 ac = numpy.zeros((nmags, nframes))
Maria@4:                 for j in range(nframes):
Maria@4:                     ac[:, j] = librosa.autocorrelate(self.melspec[i, (j*win2):min((j+1)*win2, nmelframes)], max_size=nmags)
Maria@4:                 autocor[i, :, :] = ac
Maria@4:             op = autocor
Maria@4:         self.op = op
Maria@4: 
Maria@4:     def post_process_op(self, medianfiltOP=False):
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@4:                     self.op[i,:,j] = numpy.convolve(self.op[i,:,j], numpy.ones(ks)/ks, mode='same')
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@4:         u_max = numpy.log(nmags)
Maria@4:         delta_c = numpy.pi / u_max
Maria@4:         c_max = nmagsout
Maria@4:         c = numpy.arange(delta_c, c_max, delta_c)
Maria@4:         k = range(1, nmags)
Maria@4:         exponent = 0.5 - c*1j
Maria@4: 
Maria@4:         normMat = 1./(exponent * numpy.sqrt(2*numpy.pi))
Maria@4:         normMat.shape = (normMat.shape[0], 1)
Maria@4:         normMat = numpy.repeat(normMat.T, nmels, axis=0)
Maria@4:         kernelMat = numpy.asarray([numpy.power(ki, exponent) for ki in k])
Maria@4:         opmellin = numpy.zeros((nmels, nmagsout, nframes))
Maria@4:         for i in range(nframes):
Maria@4:             self.op[:, -1, i] = 0
Maria@4:             deltaMat = - numpy.diff(self.op[:, :, i])
Maria@4:             mellin = numpy.abs(numpy.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@4:             self.opmellin = numpy.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@4:             min_opmellin = numpy.amin(self.opmellin, axis=0, keepdims=True)
Maria@4:             max_opmellin = numpy.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(medianfilt=medianfiltspec)
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@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
Maria@4:         self.onset_patterns(logfilter=False, 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@4: if __name__ == '__main__':
Maria@4:     op = OPMellin()
Maria@4:     op.get_opmellin()