comparison CAdaptInstrSpec.m @ 0:b4e26b53072f tip

Initial commit.

+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