Mercurial > hg > adaptinstrspec
changeset 0:b4e26b53072f tip
Initial commit.
| author | Holger Kirchhoff <holger.kirchhoff@eecs.qmul.ac.uk> |
|---|---|
| date | Tue, 04 Dec 2012 13:57:15 +0000 |
| parents | |
| children | |
| files | CAdaptInstrSpec.m adaptInstrSpecExample.m functions/betaDivergencePerElement.m functions/compNumPartials.m functions/dceps.m functions/freqIdx2PartialIdx.m functions/getNumBinsPerSemitoneFromFreqVec.m functions/idceps.m functions/midiPitch2Freq.m functions/midiPitch2Shift.m functions/partialIdx2FreqIdx.m instrSpectra/Violin.CQT.mat instrSpectra/Violin2.CQT.mat |
| diffstat | 13 files changed, 758 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/CAdaptInstrSpec.m Tue Dec 04 13:57:15 2012 +0000 @@ -0,0 +1,484 @@ +classdef CAdaptInstrSpec +% CAdaptInstrSpec - Estimation of a filter curve that enables the +% adaptation of instrument templates from one recording to another. The +% beta-divergence is used as a cost function between the original and the +% adapted spectra. All spectra need to be provided on a logarithmic +% frequency axis (for now, an extension to linear frequency axes should +% be straightforward). For further details on the estimation method, see +% the references below. +% +% PROPERTIES +% no public properties +% +% METHODS +% CAdaptInstrSpec - constructor for CAdaptInstrSpec object +% setH - sets filter curve 'h' +% getH - returns estimate of filter curve ''h'' +% getSmoothedH - returns smoothed and interpolated version of the +% filter curve ''h'' +% updateH - performs single update of filter curve ''h'' +% estimateSpectra - estimates spectra based on current estimate of +% the filter curve. +% compBetaDivergence - compute beta-divergence between original and +% estimated spectra +% +% For further help on the methods, type 'help CAdaptInstrSpec.[methodName]' +% +% +% References: +% +% [1] H. Kirchhoff, S. Dixon, A. Klapuri. Missing spectral templates +% estimation for user-assisted music transcription. IEEE International +% Conference on Acoustics, Speech and Signal Processing, Vancouver, +% Canada, 2013, submitted. +% [2] H. Kirchhoff, S. Dixon, and A. Klapuri. Cross-recording adaptation of +% musical instrument spectra. Technical Report C4DM-TR-11-2012, +% Queen Mary University of London, 2012. +% http://www.eecs.qmul.ac.uk/~holger/C4DM-TR-11-2012 + +% Copyright (C) 2012 Holger Kirchhoff +% +% 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. +% +% This program is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program; if not, write to the Free Software +% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + + properties (Access='private') + + h = []; % filter transfer function + maxDevFromMedianInDB = 10; % maximum deviation from mean values + numCepstralCoeffs = 20; + regularisationParam = 0.001; + + spectra_DB = []; % basis functions estimated from database (e.g. RWC) + spectra_data = []; % basis functions derived from analysis spectrogram + W_DB = []; % spectra_DB reduced to peak amplitudes only + W_data = []; % spectra_data reduced to peak amplitudes only + W_data_est = []; % estimated basis functions (W_data * h) + f0Idcs_DB = []; % pitch values of columns in W_DB + f0Idcs_data = []; % pitch values of columns in W_data + numF0Idcs_DB = 0; + numF0Idcs_data = 0; + + commonF0Idcs = []; % midi pitch values that occur both in f0Idcs_DB and in f0Idcs_data + commonF0IdcsIdcs_DB = []; % pitch indices in f0Idcs_DB that also occur in f0Idcs_data + commonF0IdcsIdcs_data = []; % pitch indices in f0Idcs_data that also occur in f0Idcs_DB + numCommonShifts = 0; + + zeroIdcsW_data = []; % indices in W_data that are zero (required for numerical reasons) + zeroIdcsW_data_est = []; % indices in W_DB that are zero + zeroIdcsH = []; % indices in h that are zero + + maxNumFreqs = 0; + + costFctName = ''; + costFctNames = {'LS', 'KL', 'IS', 'BD'}; + + beta = 0; % parameter beta for beta divergence + + end % properties + + methods + + function obj = CAdaptInstrSpec(spectra_DB, spectra_data, f0Idcs_DB, f0Idcs_data, numBinsPerSemitone, costFctName, varargin) + % CAdaptInstrSpec - constructor of CAdaptInstrSpec class + % + % myObj = CAdaptInstrSpec(spectra_DB, spectra_data, f0Idcs_DB, + % f0Idcs_data, numBinsPerSemitone, costFctName) + % constructs a CAdaptInstrSpec object. + % + % Parameters: + % spectra_DB - matrix containing in its columns the the + % database templates that are to be adapted + % spectra_data - matrix containing the spectra estimated + % from the recording to which the database + % spectra should be adapted + % f0Idcs_DB - f0-indices corresponding to the columns + % in spectra_DB + % f0Idcs_data - f0-indices corresponding to the columns + % in spectra_data + % numBinsPerSemitone - pitch resolution of the constant-Q + % spectra + % costFctName - name of cost function. available cost + % functions are: + % 'LS' - least squares error + % 'KL' - generalised Kullback-Leibler div. + % 'IS' - Itakura-Saito divergence + % 'BD' - beta divergence + % + % If 'BD' is selected as the cost function, the parameter beta + % has to be provided by myObj = CAdaptInstrSpec(..., 'beta', + % betaValue) where betaValue is a real, finite scalar. + + + assert(length(f0Idcs_DB) == size(spectra_DB,2), 'number of values in ''f0Idcs_DB'' must be the same as number of columns in W_DB'); + assert(length(f0Idcs_data) == size(spectra_data,2), 'number of values in ''f0Idcs_data'' must be the same as number of columns in W_data'); + assert(size(spectra_DB,1) == size(spectra_data,1), 'number of frequency bins (rows) in W_DB and W_data must be the same'); + + % FIXME: check that inputs are valid! + + %% set member variables + obj.f0Idcs_DB = f0Idcs_DB; + obj.f0Idcs_data = f0Idcs_data; + obj.spectra_DB = spectra_DB; + obj.spectra_data = spectra_data; + + obj.maxNumFreqs = size(spectra_DB, 1); + obj.numF0Idcs_DB = size(spectra_DB, 2); + obj.numF0Idcs_data = size(spectra_data,2); + + [obj.commonF0Idcs, obj.commonF0IdcsIdcs_DB, obj.commonF0IdcsIdcs_data] = intersect(f0Idcs_DB, f0Idcs_data); + obj.numCommonShifts = length(obj.commonF0IdcsIdcs_DB); + + + %% reduce spectra to partial amplitudes only + obj.W_DB = obj.noteSpec2partialSpec(spectra_DB, f0Idcs_DB, numBinsPerSemitone); + obj.W_data = obj.noteSpec2partialSpec(spectra_data, f0Idcs_data, numBinsPerSemitone); + + %% adjust amplitudes in database spectra at common pitches + obj.W_DB(:, obj.commonF0IdcsIdcs_DB) = obj.adjustPartialPositions(obj.W_DB(:, obj.commonF0IdcsIdcs_DB), ... + obj.W_data(:, obj.commonF0IdcsIdcs_data), ... + obj.commonF0Idcs, numBinsPerSemitone); + + + %% find zero-entries in W_data and W_data_est + obj.zeroIdcsW_data = (obj.W_data(:,obj.commonF0IdcsIdcs_data) == 0); + obj.zeroIdcsW_data_est = (obj.W_DB(:,obj.commonF0IdcsIdcs_DB) == 0); % zero where W_DB is zero (see computeW_data_est) + obj.zeroIdcsH = (sum(obj.W_data,2) == 0) | (sum(obj.W_DB(:,obj.commonF0IdcsIdcs_DB),2) == 0); + + assert(ischar(costFctName), ~any(strcmpi(costFctName, obj.costFctNames)), ... + 'Argument ''costFctName'' must be a string, and must match one of the implemented cost function names.'); + obj.costFctName = costFctName; + + + %% set beta & cost function name + switch obj.costFctName + case 'LS' + obj.beta = 2; + obj.costFctName = 'BD'; + case 'KL' + obj.beta = 1; + obj.costFctName = 'BD'; + case 'IS' + obj.beta = 0; + obj.costFctName = 'BD'; + case 'BD' + %check if beta was set + if isempty(varargin) % FIXME: if more optional arguments are added later, use MATLAB's inputparser + error('When ''BD'' is used as the cost function, beta needs to be set.'); + elseif ~strcmp(varargin{1}, 'beta') + error('Name/value pair for ''beta'' not found.') + end + beta = varargin{2}; + validateattributes(beta, {'numeric'}, {'scalar', 'real', 'finite', 'nonnan'}, 'CSourceFilter', 'beta', 5) + obj.beta = beta; + end + + %% initialise h + obj.h = zeros(obj.maxNumFreqs,1); + obj.h(~obj.zeroIdcsH) = 1; + + %% compute initial WEst + obj = computeW_data_est(obj); + + end + + + function obj = updateH(obj) + % updateH - perform single update of filter curve ''h'' + % + % myObj = myObj.updateH applies the update functions to the + % filter curve ''h''. + + %% get spectra at common f0 indices + W_data = obj.W_data(:, obj.commonF0IdcsIdcs_data); + W_DB = obj.W_DB(:, obj.commonF0IdcsIdcs_DB); + W_data_est = obj.W_data_est; + + + %% compute W_data * W_data_est^(beta-2) + nomMatrix = W_data .* W_data_est .^ (obj.beta-2); + + % fix divide by 0 + if obj.beta < 2 % if beta < 2, exponent of WEst^(beta-2) is negative -> division + maxRatio = max(max(nomMatrix( ~obj.zeroIdcsW_data_est ))); + nomMatrix(obj.zeroIdcsW_data & obj.zeroIdcsW_data_est) = 1; + nomMatrix(~obj.zeroIdcsW_data & obj.zeroIdcsW_data_est) = maxRatio; + end + + %% compute W_data^(beta-1) + denomMatrix = W_data_est .^ (obj.beta-1); + + % fix divide by 0 + if obj.beta < 1 % if beta < 1, exponent of WEst^(beta-1) is negative -> division + maxRatio = max(max(denomMatrix(~obj.zeroIdcsW_data_est))); + denomMatrix(obj.zeroIdcsW_data_est) = maxRatio; + end + + %% multiply by W_DB + nomMatrix = nomMatrix .* W_DB; + denomMatrix = denomMatrix .* W_DB; + + + %% compute nominator and denominator + nom = sum(nomMatrix, 2); + denom = sum(denomMatrix, 2); + + %% compute ratio + ratio = nom ./ denom; + ratio((nom==0) & (denom==0)) = 1; + ratio((nom~=0) & (denom==0)) = max(ratio); + + %% apply update + obj.h(~obj.zeroIdcsH) = obj.h(~obj.zeroIdcsH) .* ratio(~obj.zeroIdcsH); + + %% recompute WEst + obj = computeW_data_est(obj); + end + + + function [spectra shiftVals] = estimateSpectra(obj, shiftVals) + % estimateSpectra - computes basis functions from the current + % estimates for e and h + % + % [spectra f0Idcs] = myObj.estimateSpectra(f0Idcs) estimates the + % spectra at the f0 indices provided by ''f0Idcs'' by applying + % the current estimate of the filter curve to the spectra in + % ''spectra_DB'' (see constructor). Spectra are only returned + % for those f0Idcs that exist in ''f0Idcs_DB'' specified in the + % constructor. The output is a matrix ''spectra'' containing the + % estimated spectra and a vector ''f0Idcs'' containing the + % corresponding f0 indices. + + %FIXME: check that input variable 'shiftVals' is correct + + % find values in shiftVals that also exist in obj.f0Idcs_DB + [commonShiftVals, shiftIdcs_DB, dummy] = intersect(obj.f0Idcs_DB, shiftVals); + numCommonShiftVals = length(commonShiftVals); + + h = obj.getSmoothedH(); + %h = obj.getH(); + spectra = obj.spectra_DB(:,shiftIdcs_DB) .* repmat(h, 1, numCommonShiftVals); + end + + function h = getH(obj) + % getH - get filter curve ''h'' + % + % myH = myObj.getH() returns the member variable ''h''. + + h = obj.h; + end + + function h = getSmoothedH(obj) + % getH - get smoothed version of filter curve ''h'' + % + % myH = myObj.getSmoothedH() returns a smoothed version of the + % filter curve ''h''. Smoothing is done by applying the discrete + % cepstrum spectral envelope algorithm from Diemo Schwartz to + % the filter curve ''h''. + + nonZeroFreqIdcsH = find(~obj.zeroIdcsH); + + % select nonzero entries from h + h = obj.h; + h_nonzero = h(nonZeroFreqIdcsH); + h_nonzero_DB = 20*log10(h_nonzero); + + + % correct outliers that are more than 10 dB above or below median + medianInDB = median(h_nonzero_DB); + + idcs = h_nonzero_DB > medianInDB + obj.maxDevFromMedianInDB; + h_nonzero(idcs) = 10^( (medianInDB + obj.maxDevFromMedianInDB)/20 ); + + idcs = h_nonzero_DB < medianInDB - obj.maxDevFromMedianInDB; + h_nonzero(idcs) = 10^( (medianInDB - obj.maxDevFromMedianInDB)/20 ); + + + % setup vector containing frequencies for cosine approx. + w = (1:obj.maxNumFreqs)' / obj.maxNumFreqs * pi; + + % select w at nonzero entries of h + w_nonzero = w(nonZeroFreqIdcsH); + + % copy first and last nonzero entry to boundaries + if nonZeroFreqIdcsH(1) ~= 1 + h_nonzero = [h_nonzero(1); h_nonzero]; + w_nonzero = [w(1); w_nonzero]; + end + if nonZeroFreqIdcsH(end) ~= obj.maxNumFreqs + h_nonzero = [h_nonzero; h_nonzero(end)]; + w_nonzero = [w_nonzero; w(end)]; + end + + + % apply cosine approximation to h (discrete cepstrum) + coeffs = dceps(h_nonzero, w_nonzero, obj.numCepstralCoeffs, obj.regularisationParam); + h = idceps(coeffs, w); + + end + + function obj = setH(obj, h) + % setH - set member variable h + % + % myObj = myObj.setH(myH) sets the member variable h to myH. + % myH must be a non-negative column vector of length [number of + % frequencies]. + + obj.h = h; + obj = computeW_data_est(obj); + end + + function betaDiv = compBetaDivergence(obj) + % compBetaDivergence - computes beta divergence based on the + % current estiates + % + % betaDiv = myObj.compBetaDivergence() returns the + % beta-divergence between the instrument spectra from the + % recording and the adapted database spectra based on the value + % for beta specified by the cost function in the constructor. + + W_data = obj.W_data(~obj.zeroIdcsW_data_est); + W_data_est = obj.W_data_est(~obj.zeroIdcsW_data_est); + + switch obj.beta + case 0 + betaDivMat = W_data ./ W_data_est - log(W_data ./ W_data_est) - 1; + + case 1 + betaDivMat = W_data .* log(W_data ./ W_data_est) + W_data - W_data_est; + + otherwise + betaDivMat = (W_data .^ obj.beta) / (obj.beta * (obj.beta-1)) ... + + (W_data_est .^ obj.beta) / obj.beta ... + - (W_data .* (W_data_est .^ (obj.beta-1))) / (obj.beta-1); + end + + betaDiv = sum(betaDivMat(:)); + end + + end % methods + + + methods (Access = private) + + function obj = computeW_data_est(obj) + % computes basis functions from the current estimates for s, e and h + + if ~isempty(obj.h) + obj.W_data_est = obj.W_DB(:,obj.commonF0IdcsIdcs_DB) .* repmat(obj.h, 1, obj.numCommonShifts); + end + end + + end % methods (Access = private) + + + methods (Access = private, Static) + + function partialSpectra = noteSpec2partialSpec(noteSpectra, f0Idcs, numBinsPerSemitone) + % goes through all note spectra, extracts the partial amplitudes + % and writes them to their absolute frequency positions + + % initialize matrix for result + [numFreqs numPitches] = size(noteSpectra); + partialSpectra = zeros(numFreqs, numPitches); + + % get (ideal) relative partial positions + maxNumPartials = floor(freqIdx2PartialIdx(numFreqs, numBinsPerSemitone)); + relF0IdcsOfPartials = partialIdx2FreqIdx((1:maxNumPartials)', numBinsPerSemitone); + meansF0Idcs = geomean( [relF0IdcsOfPartials(1:end-1)'; relF0IdcsOfPartials(2:end)'] )'; + relLowerBoundOfPartials = [1; floor(meansF0Idcs)+1]; + relUpperBoundOfPartials = [floor(meansF0Idcs); numFreqs]; + + % go through all spectra + for pitchIdx = 1:numPitches + + currF0Idx = f0Idcs(pitchIdx); + currNumPartials = floor(freqIdx2PartialIdx(numFreqs-currF0Idx+1, numBinsPerSemitone)); + + % go through partials + for partialIdx = 1:currNumPartials + + % find maximum with partial range + lowerBound = currF0Idx-1 + relLowerBoundOfPartials(partialIdx); + upperBound = min(currF0Idx-1 + relUpperBoundOfPartials(partialIdx), numFreqs); + [maxAmpl maxIdx] = max(noteSpectra(lowerBound:upperBound, pitchIdx)); + + % write to result matrix + partialSpectra(lowerBound-1+maxIdx, pitchIdx) = maxAmpl; + + end + end + end % noteSpec2partialSpec + + + function W_DB = adjustPartialPositions(W_DB, W_data, f0Idcs, numBinsPerSemitone) + % adjust the positions of the partials in W_DB to those in W_data + + assert(size(W_DB,2) == size(W_data,2), '''W_DB'' and ''W_data'' must contain the same number of columns'); + assert(length(f0Idcs) == size(W_DB,2), 'Number of elements in ''f0Idcs'' must be equal to number of columns in ''W_DB'''); + + [numFreqs numPitches] = size(W_DB); + + % get (ideal) relative partial positions + maxNumPartials = floor(freqIdx2PartialIdx(numFreqs, numBinsPerSemitone)); + relF0IdcsOfPartials = partialIdx2FreqIdx((1:maxNumPartials)', numBinsPerSemitone); + meansF0Idcs = geomean( [relF0IdcsOfPartials(1:end-1)'; relF0IdcsOfPartials(2:end)'] )'; + relLowerBoundOfPartials = [-ceil(numBinsPerSemitone/2); floor(meansF0Idcs)+1]; % make 1st bound -numBinsPerSemitone/2 to allow 1st partial to deviate below ideal position + relUpperBoundOfPartials = [floor(meansF0Idcs); numFreqs]; + + % go through all spectra + for pitchIdx = 1:numPitches + + currF0Idx = f0Idcs(pitchIdx); + currNumPartials = compNumPartials(numFreqs-currF0Idx+1, numBinsPerSemitone); + + % go through partials + for partialIdx = 1:currNumPartials + + % compute bounds for partial range + lowerBound = max(currF0Idx-1 + relLowerBoundOfPartials(partialIdx), 1); + upperBound = min(currF0Idx-1 + relUpperBoundOfPartials(partialIdx), numFreqs); + + % get indices of partial in both W_DB and W_data + idx_DB = find(W_DB(lowerBound:upperBound, pitchIdx)); + idx_data = find(W_data(lowerBound:upperBound, pitchIdx)); + + % set partial amplitude in W_DB to frequency index of W_data + partAmpl = W_DB(lowerBound-1+idx_DB, pitchIdx); + W_DB(lowerBound-1+idx_DB, pitchIdx) = 0; + W_DB(lowerBound-1+idx_data, pitchIdx) = partAmpl; + end + end + end % adjustPartialPositions + end + +% methods (Static) +% +% function Xshift = shiftSpectra(X, shiftVals) +% % shifts each spectrum (column in X) down by the amount specified in +% % shiftVals +% +% assert(size(X,2) == length(shiftVals), 'number values in ''shiftVals'' must be the same as number of columns in ''X'''); +% +% [maxNumFreqs numShifts] = size(X); +% Xshift = zeros(maxNumFreqs, numShifts); +% +% for pitchIdx = 1:numShifts +% phi = shiftVals(pitchIdx); +% Xshift(1:maxNumFreqs-phi, pitchIdx) = X(phi+1:maxNumFreqs, pitchIdx); +% end +% end +% end % methods (Static) + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/adaptInstrSpecExample.m Tue Dec 04 13:57:15 2012 +0000 @@ -0,0 +1,135 @@ +% This is an example script to illustrate the use of the CAdaptInstrSpec +% class. The script loads the spectra extracted two recordings of two +% different violins. It then applied the CAdaptInstrSpec class to adapt the +% spectra of the first violin to those of the second. + +clearvars; +close all; + +addpath('./functions/'); + + +%% constants + +% file paths of instrument spectra +noteSpectraFilePaths = {'./instrSpectra/Violin.CQT.mat'; ... + './instrSpectra/Violin2.CQT.mat'}; +numInstr = length(noteSpectraFilePaths); + +% analysis +tuningFreqInHz = 440; +costFctName = 'LS'; +beta = 2; % only needed for plotting. should match the cost function above +numIter = 5; + + + + +%% load instrument spectra +instrSpec = cell(numInstr,1); +f0Idcs = cell(numInstr,1); + +for instrIdx = 1:numInstr + load(noteSpectraFilePaths{instrIdx}); % loads variables: noteSpectra, midiPitches, freqsInHz + instrSpec{instrIdx} = noteSpectra; + f0Idcs{instrIdx} = midiPitch2Shift(midiPitches, tuningFreqInHz, freqsInHz)+1; +end +numBinsPerSemitone = getNumBinsPerSemitoneFromFreqVec(freqsInHz); + + +%% estimate filter curve + +% create instrSpecFilter object +instrFiltObj = CAdaptInstrSpec(instrSpec{1}, instrSpec{2}, f0Idcs{1}, f0Idcs{2}, numBinsPerSemitone, costFctName); + +% apply update function for h +for iterIdx = 1:numIter; + instrFiltObj = instrFiltObj.updateH; +end + +h = instrFiltObj.getH; +h_smooth = instrFiltObj.getSmoothedH; + + +%% estimate spectra with given filter curve +estSpectra = instrFiltObj.estimateSpectra(f0Idcs{1}); + + + + +%% plot results + +% plot original spectra and estimated spectra +yticks = [250 500 1000 2000 4000]; +yticklabel = {'250'; '500'; '1k'; '2k'; '4k'}; + +figure; +subplot(311); +imagesc(f0Idcs{1}, freqsInHz, instrSpec{1}); +axis xy; +set(gca, 'YScale', 'log', 'YTick', yticks, 'YTickLabel', yticklabel); +xlabel('pitch index'); +ylabel('frequency [Hz]'); +title('spectra of instrument 1'); + + +subplot(312); +imagesc(f0Idcs{1}, freqsInHz, instrSpec{2}); +axis xy; +set(gca, 'YScale', 'log', 'YTick', yticks, 'YTickLabel', yticklabel); +xlabel('pitch index'); +ylabel('frequency [Hz]'); +title('spectra of instrument 2'); + +subplot(313); +imagesc(f0Idcs{1}, freqsInHz, estSpectra); +axis xy; +set(gca, 'YScale', 'log', 'YTick', yticks, 'YTickLabel', yticklabel); +xlabel('pitch index'); +ylabel('frequency [Hz]'); +title('spectra of instrument 2 adapted to instrument 1') + + +% plot filter curve and errors +xticks = [125 250 500 1000 2000 4000]; +xticklabel = {'125'; '250'; '500'; '1k'; '2k'; '4k'}; + +% compute beta divergence for each element +betaDivsOriginal = betaDivergencePerElement(instrSpec{1}, instrSpec{2}, beta); +betaDivsEstimate = betaDivergencePerElement(instrSpec{1}, estSpectra, beta); + +% scale per-element beta divergence for image plot +numColors = size(colormap,1); +maxDist = max([betaDivsOriginal(:); betaDivsEstimate(:)]); +betaDivsOriginal = betaDivsOriginal / maxDist * numColors; +betaDivsEstimate = betaDivsEstimate / maxDist * numColors; + + + +figure; +subplot(311) +%stem(freqsInHz(h ~= 0), h(h ~= 0), 'k', 'filled'); +plot(freqsInHz(h ~= 0), h(h ~= 0), 'k.'); +hold on; +plot(freqsInHz, h_smooth, 'k'); +hold off; +axis tight; +set(gca, 'XScale', 'log', 'XTick', xticks, 'XTickLabel', xticklabel); +title('estimated filter curve ''h'''); +legend('original', 'smoothed', 'Location', 'NorthWest'); + +subplot(312); +image(f0Idcs{1}, freqsInHz, betaDivsOriginal); +axis xy; +set(gca, 'YScale', 'log', 'YTick', yticks, 'YTickLabel', yticklabel); +xlabel('pitch index'); +ylabel('frequency [Hz]'); +title('elementwise differences between original sets of spectra'); + +subplot(313); +image(f0Idcs{1}, freqsInHz, betaDivsEstimate); +axis xy; +set(gca, 'YScale', 'log', 'YTick', yticks, 'YTickLabel', yticklabel); +xlabel('pitch index'); +ylabel('frequency [Hz]'); +title('differences between original and adapted set of spectra'); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/functions/betaDivergencePerElement.m Tue Dec 04 13:57:15 2012 +0000 @@ -0,0 +1,14 @@ +function betaDivPerElem = betaDivergencePerElement(x, y, beta) + +switch beta + case 0 + betaDivPerElem = x./y - log(x./y) - 1; + + case 1 + betaDivPerElem = x .* log(x./y) + x - y; + + otherwise + betaDivPerElem = (x.^beta) / (beta*(beta-1)) ... + + (y.^beta) / beta ... + - (x .* (y.^(beta-1))) / (beta-1); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/functions/compNumPartials.m Tue Dec 04 13:57:15 2012 +0000 @@ -0,0 +1,6 @@ +function numPartials = compNumPartials(numFreqs, numBinsPerSemitone) +% numPartials = compNumPartials(numFreqs, numBinsPerSemitone) +% +% compute the number of partials that fall within the spectrogram range + +numPartials = floor( 2^(numFreqs / (12 * numBinsPerSemitone)) ); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/functions/dceps.m Tue Dec 04 13:57:15 2012 +0000 @@ -0,0 +1,47 @@ +function c = dceps(X, w, numCoeffs, lambda) +% computes the discrete cepstral coefficients of a sequence of partials X. +% +% c = dceps(X, w, numCoeffs) +% +% w contains the normalized frequencies of the partials. Normalized +% frequencies should be in the range [0 .. pi] +% X and w must be column vectors. +% numCoeffs specifies the number of cepstral coefficients. +% +% Further details can be found in: +% [1] Diemo Schwarz. Spectral envelopes in sound analysis and synthesis. +% Master's thesis, Universitaet Stuttgart, 1998, pp. 36-38. +% [2] T. Galas and X. Rodet. An improved cepstral method for deconvolution +% of source filter systems with discrete spectra: Application to musical +% sound signals. In International Computer Music Conference, 1990. + + +p = numCoeffs-1; +n = length(X); % number of partials +S = ones(n,1); % amplitude of partials of source spectrum +%lambda = 0.0005; % regularization term + +if sum(X == 0) %if one or more element is zero + error('All elements of X must be nonzero'); +end + +% compute all r_i (i from 0 to 2p) +r = .5 * sum(cos(w * (0:2*p))); + +% compute matrix A +A = zeros(p+1); +B = zeros(p+1); +for i = 0:p + for j = 0:p + A(i+1,j+1) = r(i+j+1) + r(abs(i-j)+1); % +1 because 1st idx in r is 0! + end + B(i+1,i+1) = 8 * (pi*(i-1))^2; +end + +A = A + lambda * B; + +% compute b +b = sum( repmat(log(X./S),1,p+1) .* cos(w *(0:p)) )'; + +% compute c +c = A \ b; % equals inv(A) * b \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/functions/freqIdx2PartialIdx.m Tue Dec 04 13:57:15 2012 +0000 @@ -0,0 +1,5 @@ +function partialIdcs = freqIdx2PartialIdx(freqIdcs, numBinsPerSemitone) +% converts frequency indices of the basis functions of the svNMD framework +% into indices of harmonic partials (! not tested !) + +partialIdcs = 2.^( (freqIdcs-1)/(12*numBinsPerSemitone) ); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/functions/getNumBinsPerSemitoneFromFreqVec.m Tue Dec 04 13:57:15 2012 +0000 @@ -0,0 +1,9 @@ +function numBinsPerSemitone = getNumBinsPerSemitoneFromFreqVec(freqsInHz) +% numBinsPerSemitone = getNumBinsPerSemitoneFromFreqVec(freqsInHz) +% +% infers the number of bins per semitone from a given constant-Q frequency +% vector in Hz + +ratios = freqsInHz(2:end) ./ freqsInHz(1:end-1); +avgRatio = mean(ratios); +numBinsPerSemitone = round( 1/(12*log2(avgRatio)) ); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/functions/idceps.m Tue Dec 04 13:57:15 2012 +0000 @@ -0,0 +1,22 @@ +function X = idceps(c, w) +% X = idceps(c, w) +% +% computes the inverse discrete cepstrum, i.e. given a number of discrete +% cepstral coefficients, it returns the spectrum. +% +% c contains the discrete cepstral coefficients +% w contains the frequency at which the spectrum is evaluated. +% c and w must be column vectors. +% +% Further details can be found in: +% [1] Diemo Schwarz. Spectral envelopes in sound analysis and synthesis. +% Master's thesis, Universitaet Stuttgart, 1998, pp. 36-38. +% [2] T. Galas and X. Rodet. An improved cepstral method for deconvolution +% of source filter systems with discrete spectra: Application to musical +% sound signals. In International Computer Music Conference, 1990. + +pmax = length(c)-1; +numW = length(w); + +exponent = sum( repmat(c', numW, 1) .* cos(w*(0:pmax)), 2 ); +X = exp(exponent); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/functions/midiPitch2Freq.m Tue Dec 04 13:57:15 2012 +0000 @@ -0,0 +1,6 @@ +function freq = midiPitch2Freq(midiPitch, tuningFreqInHz) +% freq = midiPitch2Freq(midiPitch, tuningFreqInHz) +% +% converts midi pitch index to frequency given a tuning frequeny + +freq = tuningFreqInHz * 2.^((midiPitch - 69)/12); % midi pitch 69 belongs to A4, the standard pitch, for which the tuning freq is given
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/functions/midiPitch2Shift.m Tue Dec 04 13:57:15 2012 +0000 @@ -0,0 +1,23 @@ +function shifts = midiPitch2Shift(midiPitches, tuningFreqInHz, binFreqsCQT) +% shifts = midiPitch2Shift(midiPitches, tuningFreqInHz, binFreqsCQT) +% +% converts midi pitch to shift values in constant-Q vector. +% shifts for pitches that below the lowest or above the highest CQT +% frequency are set to a value of -1 + +numMidiPitches = length(midiPitches); +numFreqsCQT = length(binFreqsCQT); + +pitchF0s = midiPitch2Freq(midiPitches, tuningFreqInHz); +outOfBoundsPitches = (pitchF0s < binFreqsCQT(1) | pitchF0s > binFreqsCQT(end)); +if any(outOfBoundsPitches) + warning('F0s of some midi pitches lie outside the constant-Q spectrogram frequency range'); +end + +freqRatios = repmat(binFreqsCQT, 1, numMidiPitches) ./ repmat(pitchF0s', numFreqsCQT, 1); +freqRatios(freqRatios < 1) = 1 ./ freqRatios(freqRatios < 1); + +[dummy shifts] = min(abs(freqRatios), [], 1); +shifts = shifts - 1; % 'no shift' is defined as 0 +shifts(outOfBoundsPitches) = -1; +shifts = shifts'; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/functions/partialIdx2FreqIdx.m Tue Dec 04 13:57:15 2012 +0000 @@ -0,0 +1,7 @@ +function freqIdcs = partialIdx2FreqIdx(partialIdcs, numBinsPerSemitone) +% freqIdcs = partialIdx2FreqIdx(partialIdcs, numBinsPerSemitone) +% +% converts indices of harmonic partials into frequency indices of the +% basis functions of the svNMD framework + +freqIdcs = 12 * numBinsPerSemitone * log2(partialIdcs) + 1; \ No newline at end of file