danstowell@1: '''MFCC.py
danstowell@1: 
danstowell@1: Calculation of MFCC coefficients from frequency-domain data
danstowell@1: 
danstowell@1: Adapted from the Vampy example plugin "PyMFCC" by Gyorgy Fazekas
danstowell@1: http://code.soundsoftware.ac.uk/projects/vampy/repository/entry/Example%20VamPy%20plugins/PyMFCC.py
danstowell@1: 
danstowell@1: Centre for Digital Music, Queen Mary University of London.
danstowell@1: Copyright (C) 2009 Gyorgy Fazekas, QMUL.
danstowell@1: '''
danstowell@1: 
danstowell@1: import sys,numpy
danstowell@1: from numpy import abs,log,exp,floor,sum,sqrt,cos,hstack
danstowell@1: from numpy.fft import *
danstowell@1: 
danstowell@1: class melScaling(object):
danstowell@1: 
danstowell@1: 	def __init__(self,sampleRate,inputSize,numBands,minHz = 0,maxHz = None):
danstowell@1: 		'''Initialise frequency warping and DCT matrix. 
danstowell@1: 		Parameters:
danstowell@1: 		sampleRate: audio sample rate
danstowell@1: 		inputSize: length of magnitude spectrum (half of FFT size assumed)
danstowell@1: 		numBands: number of mel Bands (MFCCs)
danstowell@1: 		minHz: lower bound of warping  (default = DC)
danstowell@1: 		maxHz: higher bound of warping (default = Nyquist frequency)
danstowell@1: 		'''
danstowell@1: 		self.sampleRate = sampleRate
danstowell@1: 		self.NqHz = sampleRate / 2.0
danstowell@1: 		self.minHz = minHz
danstowell@1: 		if maxHz is None : maxHz = self.NqHz
danstowell@1: 		self.maxHz = maxHz
danstowell@1: 		self.inputSize = inputSize
danstowell@1: 		self.numBands = numBands
danstowell@1: 		self.valid = False
danstowell@1: 		self.updated = False
danstowell@1: 		
danstowell@1: 	def update(self): 
danstowell@1: 		# make sure this will run only once 
danstowell@1: 		# if called from a vamp process
danstowell@1: 		if self.updated: return self.valid
danstowell@1: 		self.updated = True
danstowell@1: 		self.valid = False
danstowell@8: 		print('Updating parameters and recalculating filters: ')
danstowell@8: 		print('Nyquist: ',self.NqHz)
danstowell@1: 		
danstowell@1: 		if self.maxHz > self.NqHz : 
danstowell@1: 			raise Exception('Maximum frequency must be smaller than the Nyquist frequency')
danstowell@1: 		
danstowell@1: 		self.maxMel = 1000*log(1+self.maxHz/700.0)/log(1+1000.0/700.0)
danstowell@1: 		self.minMel = 1000*log(1+self.minHz/700.0)/log(1+1000.0/700.0)
danstowell@8: 		print('minHz:%s\nmaxHz:%s\nminMel:%s\nmaxMel:%s\n' \
danstowell@8: 		%(self.minHz,self.maxHz,self.minMel,self.maxMel))
danstowell@1: 		self.filterMatrix = self.getFilterMatrix(self.inputSize,self.numBands)
danstowell@1: 		self.DCTMatrix = self.getDCTMatrix(self.numBands)
danstowell@1: 		self.valid = True
danstowell@1: 		return self.valid
danstowell@1: 				
danstowell@1: 	def getFilterCentres(self,inputSize,numBands):
danstowell@1: 		'''Calculate Mel filter centres around FFT bins.
danstowell@1: 		This function calculates two extra bands at the edges for
danstowell@1: 		finding the starting and end point of the first and last 
danstowell@1: 		actual filters.'''
danstowell@8: 		centresMel = numpy.array(range(numBands+2)) * (self.maxMel-self.minMel)/(numBands+1) + self.minMel
danstowell@1: 		centresBin = numpy.floor(0.5 + 700.0*inputSize*(exp(centresMel*log(1+1000.0/700.0)/1000.0)-1)/self.NqHz)
danstowell@1: 		return numpy.array(centresBin,int)
danstowell@1: 		
danstowell@1: 	def getFilterMatrix(self,inputSize,numBands):
danstowell@1: 		'''Compose the Mel scaling matrix.'''
danstowell@1: 		filterMatrix = numpy.zeros((numBands,inputSize))
danstowell@1: 		self.filterCentres = self.getFilterCentres(inputSize,numBands)
danstowell@8: 		for i in range(numBands) :
danstowell@1: 			start,centre,end = self.filterCentres[i:i+3]
danstowell@1: 			self.setFilter(filterMatrix[i],start,centre,end)
danstowell@1: 		return filterMatrix.transpose()
danstowell@1: 
danstowell@1: 	def setFilter(self,filt,filterStart,filterCentre,filterEnd):
danstowell@1: 		'''Calculate a single Mel filter.'''
danstowell@1: 		k1 = numpy.float32(filterCentre-filterStart)
danstowell@1: 		k2 = numpy.float32(filterEnd-filterCentre)
danstowell@8: 		up = (numpy.array(range(filterStart,filterCentre))-filterStart)/k1
danstowell@8: 		dn = (filterEnd-numpy.array(range(filterCentre,filterEnd)))/k2
danstowell@1: 		filt[filterStart:filterCentre] = up
danstowell@1: 		filt[filterCentre:filterEnd] = dn
danstowell@3: 
danstowell@1: 	def warpSpectrum(self,magnitudeSpectrum):
danstowell@1: 		'''Compute the Mel scaled spectrum.'''
danstowell@1: 		return numpy.dot(magnitudeSpectrum,self.filterMatrix)
danstowell@1: 		
danstowell@1: 	def getDCTMatrix(self,size):
danstowell@1: 		'''Calculate the square DCT transform matrix. Results are 
danstowell@1: 		equivalent to Matlab dctmtx(n) with 64 bit precision.'''
danstowell@8: 		DCTmx = numpy.array(range(size),numpy.float64).repeat(size).reshape(size,size)
danstowell@1: 		DCTmxT = numpy.pi * (DCTmx.transpose()+0.5) / size
danstowell@1: 		DCTmxT = (1.0/sqrt( size / 2.0)) * cos(DCTmx * DCTmxT)
danstowell@1: 		DCTmxT[0] = DCTmxT[0] * (sqrt(2.0)/2.0)
danstowell@1: 		return DCTmxT
danstowell@1: 		
danstowell@1: 	def dct(self,data_matrix):
danstowell@1: 		'''Compute DCT of input matrix.'''
danstowell@1: 		return numpy.dot(self.DCTMatrix,data_matrix)
danstowell@1: 		
danstowell@1: 	def getMFCCs(self,warpedSpectrum,cn=True):
danstowell@1: 		'''Compute MFCC coefficients from Mel warped magnitude spectrum.'''
danstowell@24: 		mfccs=self.dct(numpy.log(numpy.clip(warpedSpectrum, 1e-9, numpy.inf)))
danstowell@1: 		if cn is False : mfccs[0] = 0.0
danstowell@1: 		return mfccs
danstowell@1: