view util/smoothiecore.py @ 5:e8ebba3d6294

add smoothie

"""Smoothie

A module aimed at melody salience and melody features. It contains

* a frequency transform inspired by constant-Q transforms
* an NMF-based melodic salience function
* methods that transform this salience into melody features (not implemented)

This program is free software; you can redistribute it and/or modify it under 
the terms of the GNU General Public License as published by the Free Software 
Foundation; either version 2 of the License, or (at your option) any later 
version. See the file COPYING included with this distribution for more information.

author: Matthias Mauch 
(addition of wrap_to_octave and get_chroma functions by Maria Panteli)

"""

import sys
import numpy as np
from scipy.signal import hann
from scipy.sparse import csr_matrix
import librosa


p = {# smoothie q mapping parameters
     'bpo'           :     36,
     'break1'        :    500,
     'break2'        :   4000,
     'f_max'         :   8000,
     'fs'            :  16000,
     # spectrogram parameters
     'step_size'     :    160,
     'use_A_wtg'     :  False,
     'harm_remove'   :  False,
     # NMF and NMF dictionary parameters
     's'             :      0.85, # decay parameter
     'n_harm'        :     50, # I'm gonna do you no harm
     'min_midi'      :     25,
     'max_midi'      :     90.8,
     'bps'           :      5,
     'n_iter'        :     30
     }


def update_nmf_simple(H, W, V):
	"""Update the gain matrix H using the multiplicative NMF update equation.

	Keyword arguments:
	H -- gain matrix
	W -- template matrix
	V -- target data matrix
	"""
	WHV = np.dot(W, H) ** (-1) * V
	H = H * np.dot(W.transpose(), WHV)
	return H

def make_a_weight(f):
	"""Return the values of A-weighting at the input frequencies f."""
	f = np.array(f)
	zaehler = 12200 * 12200* (f**4)
	nenner = (f*f + 20.6*20.6) * np.sqrt((f*f + 107.7*107.7) * (f*f + 737.9*737.9)) * (f*f + 12200*12200)
	return zaehler / nenner


def nextpow2(x):
    """Return the smallest integer n such that 2 ** n > x."""
    return np.ceil(np.log2(x))


def get_smoothie_frequencies(p):
	"""Calculate filter centres and widths for the smooth Q transform."""
	# I think this calculation should be reconsidered in order to have
	# a more principled transition between linear and exponential
	n = 1.0 / (2.0**(1.0/p['bpo']) - 1)
	x = p['break1'] * 1.0 / n

	# first linear stretch
	# f = (np.arange(1, int(np.round(n)) + 1) * x).round().tolist()
	f = (np.arange(1, int(np.round(n)) + 1) * x).tolist()

	# exponential stretch
	while max(f) < p['break2']:
		f.append(max(f) * 2**(1.0/p['bpo']))

	# final linear stretch
	lastdiff = f[-1] - f[-2]
	while max(f) < p['f_max']:
		f.append(max(f) + lastdiff)

	deltaf = np.diff(np.array(f))
	f = f[:-1]

	return f, deltaf


def create_smoothie_kernel(f, deltaf, fs):
	"""Create a sparse matrix that maps the complex DFT to the complex 
	smoothie representation.
	"""
	print >>sys.stdout, "[ SMOOTHIE Q kernel calculation ... ]"
	n_filter = len(f)
	n_fft = 2**nextpow2(np.ceil(fs/min(deltaf)))
	
	thresh = 0.0054
	smoothie_kernel = np.zeros([n_fft, n_filter], np.complex64)

	for i_filter in range(n_filter-1, -1, -1): # descending
		# print i_filter
		Q = f[i_filter] * 1.0 / deltaf[i_filter] # local Q for this filter
		lgth = int(np.ceil(fs * 1.0 / deltaf[i_filter]))
		lgthinv = 1.0 / lgth
		offset = int(n_fft/2 - np.ceil(lgth * 0.5))
		temp = hann(lgth) * lgthinv * \
			np.exp(2j * np.pi * Q * (np.arange(lgth) - offset) * lgthinv)
		# print(sum(hann(lgth)), Q, lgth, offset)
		temp_kernel = np.zeros(n_fft, dtype = np.complex64)
		temp_kernel[np.arange(lgth) + offset] = temp
		spec_kernel = np.fft.fft(temp_kernel)
		spec_kernel[abs(spec_kernel) <= thresh] = 0
		smoothie_kernel[:,i_filter] = spec_kernel
	return csr_matrix(smoothie_kernel.conj()/n_fft)


def smoothie_q_spectrogram(x, p):
	"""Calculate the actual spectrogram with smooth Q frequencies"""
	print >>sys.stdout, "[ SMOOTHIE Q spectrogram ... ]"

	# precalculate smoothie kernel
	f, deltaf = get_smoothie_frequencies(p)
	smoothie_kernel = create_smoothie_kernel(f, deltaf, p['fs'])
	n_fft, n_filter = smoothie_kernel.shape
	
	# some preparations
	n_sample = len(x)
	# n_frame = int(np.floor(n_sample / p['step_size']))
	n_frame = int(np.ceil(n_sample / float(p['step_size']))) # added mp
	t = (np.arange(n_frame) * p['step_size']) * 1.0 / p['fs']
	smoothie_kernelT = smoothie_kernel.T

	# allocate
	s = np.zeros((n_filter, n_frame), dtype = np.complex64)

	# pad (if wanted)
	x = np.concatenate((np.zeros(n_fft/2), x, np.zeros(n_fft/2)))
	
	for i_frame in range(n_frame):
		smpl = p['step_size'] * i_frame
		block = x[smpl + np.arange(n_fft)]
		s[:, i_frame] = smoothie_kernelT.dot(np.fft.fft(block))

	if p['use_A_wtg']:
		a_wtg = make_a_weight(f)
		s = s * a_wtg[:, np.newaxis]

	return s, t


def mel_triangles(input_f):
	"""Make matrix with mel filters at the given frequencies.
	Warning: this is a very coarse mel filterbank.
	"""
	n_linearfilters = 3
	n_logfilters0 = 30 # just something high, will be pruned later
	linear_spacing = 100
	log_spacing = 6.0/4
	n_filter0 = n_linearfilters + n_logfilters0

	freqs = np.zeros(n_filter0+2) # includes one more on either side, hence +2
	freqs[range(n_linearfilters+1)] = \
		np.arange(-1,n_linearfilters) * linear_spacing
	freqs[range(n_linearfilters+1, n_filter0+2)] = \
		freqs[n_linearfilters] * log_spacing ** np.arange(1, n_logfilters0+2)

	centre_freqs = freqs[1:-1]
	lower_freqs  = freqs[0:-2]
	upper_freqs = freqs[2:]

	n_filter = list(np.nonzero(lower_freqs < max(input_f)))[0][-1] + 1

	n_input_f = len(input_f)
	mtr = np.zeros((n_input_f, n_filter))

	for i_filter in range(n_filter):
		for i_freq in range(n_input_f):
			if input_f[i_freq] > lower_freqs[i_filter] \
				and input_f[i_freq] <= upper_freqs[i_filter]:

				if input_f[i_freq] <= centre_freqs[i_filter]:
					mtr[i_freq, i_filter] = \
						(input_f[i_freq] - lower_freqs[i_filter]) * 1.0 / \
						(centre_freqs[i_filter] - lower_freqs[i_filter])
				else:
					mtr[i_freq, i_filter] = \
						1 - (input_f[i_freq] - centre_freqs[i_filter]) / \
						(upper_freqs[i_filter] - centre_freqs[i_filter])
	return mtr


def create_smoothie_nfm_dict(smoothie_kernel, filterf, p):
	"""Create dictionary matrix with note templates."""
	n_note = int((p['max_midi'] - p['min_midi']) * p['bps'] + 1)
	n_fft, n_filter = smoothie_kernel.shape
	t = ((np.arange(n_fft) + 1) - ((n_fft + 1)*0.5))/p['fs']

	mtr = mel_triangles(filterf)
	n_template = n_note + mtr.shape[1]

	w = np.zeros((n_filter, n_template), dtype = np.complex64)
	w[:,n_note:] = mtr
	f0s = []

	smoothie_kernelT = smoothie_kernel.T

	for i_note in range(n_note):
		midi = p['min_midi'] + i_note * 1.0 / p['bps']
		f0 = 440 * 2 ** ((midi-69)*1.0/12)
		f0s.append(f0)
		sig = np.zeros(len(t))
		for i_harm in range(p['n_harm']):
			f = f0 * (i_harm + 1)
			if f > p['fs'] * 0.5:
				continue
			x = np.sin(2*np.pi*f*t) * p['s']**(i_harm)
			sig += x
		w[:, i_note] = smoothie_kernelT.dot(np.fft.fft(sig))

	for i in range(mtr.shape[1]):
		f0s.append(np.nan)

	w = abs(w)
	col_sums = w.sum(axis = 0)
	w = w / col_sums[np.newaxis, :] # normalisation
	return w, np.array(f0s)


def smoothie_salience(x, p, do_tune = False):
	"""Calculate melodic salience."""
	print >>sys.stdout, "[ SMOOTHIE Q salience ... ]"

	# precalculate nmf kernel
	f, deltaf = get_smoothie_frequencies(p)
	smoothie_kernel = create_smoothie_kernel(f, deltaf, p['fs'])
	w, f0s = create_smoothie_nfm_dict(smoothie_kernel, f, p)

	# calculate smoothiegram
	s, t = smoothie_q_spectrogram(x, p)
	s[s==0] = np.spacing(1) # very small number
	s = abs(s)

	# run NMF
	n_frame = len(t)
	print >>sys.stdout, "[ SMOOTHIE Q : NMF updates ... ]"
	sal = np.ones((w.shape[1], n_frame))
	for i_iter in range(p['n_iter']):
		sal = update_nmf_simple(sal, w, s)

	#if do_tune:
	#	error("Tuning isn't yet implemented in the Python version")

	return sal, t, f0s


def qnormalise(x, q, dim):
    """Normalise by sum."""
    nrm = np.sum(x**q, axis=dim, keepdims=True)**(1./float(q))
    nrmmatrix = np.repeat(nrm, x.shape[0], axis=0)
    x = x / nrmmatrix
    return x


def wrap_to_octave(sal, param):
    """Wrap pitch to octave."""
    epsilon = 0.00001
    nsal = qnormalise(sal, 2, 0)
    
    step = 1. / float(param['bps'])
    NMFpitch = np.arange(param['min_midi'], param['max_midi'] + step, step)
    nNote = len(NMFpitch)
    nsal = nsal[0:nNote, :]

    allsal = nsal
    allsal[nsal < epsilon] = epsilon
    
    # chroma mapping
    n = param['bps'] * 12
    mmap = np.tile(np.eye(n), (1, 5))
    chroma = mmap.dot(allsal[0:(n * 5), :])
    return chroma


def get_chroma(filename = 'test.wav'):
    """Get chroma from audio file."""
    y, sr = librosa.load(filename, sr = None)
    p['fs'] = sr
    p['step_size'] = int(round(0.005 * p['fs']))  # default hop size 0.005 sec.
    sal, t, f0s = smoothie_salience(y, p)
    sal = sal[-np.isnan(f0s), :]
    chroma = wrap_to_octave(sal, p)
    chromasr = p['fs'] / float(p['step_size'])
    return chroma, chromasr