Chris@0: function [shiftFB,centerA4,tuningSemitones,sideinfo] = estimateTuning(f_input,parameter,sideinfo)
Chris@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Chris@0: % Name: estimateTuning
Chris@0: % Date of Revision: 2011-03
Chris@0: % Programmer: Sebastian Ewert
Chris@0: %
Chris@0: % Description:
Chris@0: % - input is either a mono audio signal or a real valued spectrogram (magnitude or power)
Chris@0: % - if input is a spectrogram you have to set parameter.timeinfo and
Chris@0: % parameter.freqinfo (vectors defining the time and freq centers for each spectrogram bin)
Chris@0: % - if input is an audio signal, a sampling freq of 22050 Hz is assumed
Chris@0: % - guesses the tuning according to a simple energy maximizing criterion
Chris@0: % - output is either: what shiftFB is best to use (shiftFB \in [0:5]).
Chris@0: %   Alternatively, the center freq for A4 is given which can be used to
Chris@0: %   specify a filterbank on your own. The second option is more fine
Chris@0: %   grained.
Chris@0: %   Alternatively, it gives a tuning in semitones, which can
Chris@0: %   easily be shifted cyclicly. For example: a tuning of -19/20 is more likely to be
Chris@0: %   +1/20 Tuning difference.
Chris@0: % - parameter.numAdditionalTunings: how many tunings besides the fixed shiftFB ones
Chris@0: % to test. For example: If set to 3, than three additional tuning settings are
Chris@0: % tested for, located at 1/4, 2/4 and 3/4 semitones below the reference
Chris@0: % tuning. If set to 5, then at 1/6, 2/6,..., 5/6 semitones.
Chris@0: % - parameter.pitchRange specifies which pitches are considered for the
Chris@0: % tuning estimation. 
Chris@0: % - parameter.pitchWeights: each pitch is considered according to a weight
Chris@0: % - Middle pitches are considered as being more important per default because
Chris@0: % here the frequency resolution is high enough. Additionally the piano has
Chris@0: % a consistent tuning only for middle pitches.
Chris@0: %
Chris@0: % Input:
Chris@0: %        f_input
Chris@0: %        parameter.numAdditionalTunings = 0;
Chris@0: %        parameter.pitchRange = [21:108];
Chris@0: %        parameter.pitchWeights = gausswin(length(parameter.pitchRange)).^2;  
Chris@0: %        parameter.fftWindowLength = 8192;
Chris@0: %        parameter.windowFunction = @hanning;
Chris@0: %        sideinfo
Chris@0: %
Chris@0: % Output:
Chris@0: %        shiftFB
Chris@0: %        centerA4
Chris@0: %        tuningSemitones
Chris@0: %        sideinfo
Chris@0: %
Chris@0: % License:
Chris@0: %     This file is part of 'Chroma Toolbox'.
Chris@0: % 
Chris@0: %     'Chroma Toolbox' is free software: you can redistribute it and/or modify
Chris@0: %     it under the terms of the GNU General Public License as published by
Chris@0: %     the Free Software Foundation, either version 2 of the License, or
Chris@0: %     (at your option) any later version.
Chris@0: % 
Chris@0: %     'Chroma Toolbox' is distributed in the hope that it will be useful,
Chris@0: %     but WITHOUT ANY WARRANTY; without even the implied warranty of
Chris@0: %     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Chris@0: %     GNU General Public License for more details.
Chris@0: % 
Chris@0: %     You should have received a copy of the GNU General Public License
Chris@0: %     along with 'Chroma Toolbox'. If not, see <http://www.gnu.org/licenses/>.
Chris@0: % 
Chris@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Chris@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Chris@0: % Check parameters
Chris@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Chris@0: 
Chris@0: if nargin<2
Chris@0:     parameter=[];
Chris@0: end
Chris@0: if nargin<1
Chris@0:     error('Please specify input data');
Chris@0: end
Chris@0: 
Chris@0: if isfield(parameter,'numAdditionalTunings')==0
Chris@0:     parameter.numAdditionalTunings = 0;
Chris@0: end
Chris@0: if isfield(parameter,'pitchRange')==0
Chris@0:     % Which pitches to consider during the estimation
Chris@0:     parameter.pitchRange = [21:108];
Chris@0: end
Chris@0: if isfield(parameter,'pitchWeights')==0
Chris@0:     % assign a weight to each pitch specified in parameter.pitchRange to
Chris@0:     % specify it's importance
Chris@0:     parameter.pitchWeights = gausswin(length(parameter.pitchRange)).^2;  
Chris@0: end
Chris@0: 
Chris@0: % the following parameters are only for audio signal input
Chris@0: if isfield(parameter,'fftWindowLength')==0
Chris@0:     parameter.fftWindowLength = 8192;
Chris@0: end
Chris@0: if isfield(parameter,'windowFunction')==0
Chris@0:     parameter.windowFunction = @hanning;  % only tested with hanning. 
Chris@0: end
Chris@0: 
Chris@0: if min(size(f_input)) == 1
Chris@0:     inputIsAudioSignal = 1;
Chris@0: else
Chris@0:     inputIsAudioSignal = 0;
Chris@0: end
Chris@0: 
Chris@0: if ~inputIsAudioSignal
Chris@0:     if isfield(parameter,'freqinfo')==0
Chris@0:         error('When using a spectrogram input you have to set parameter.freqinfo');
Chris@0:     end
Chris@0: end
Chris@0: 
Chris@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Chris@0: % Main program
Chris@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Chris@0: 
Chris@0: numTunings = 6 + parameter.numAdditionalTunings;
Chris@0: referenceFreqsA4 = zeros(numTunings,1);
Chris@0: tunings = zeros(numTunings,1);
Chris@0: 
Chris@0: tunings(1) = 0;
Chris@0: tunings(2) = -1/4;
Chris@0: tunings(3) = -1/3;
Chris@0: tunings(4) = -1/2;
Chris@0: tunings(5) = -2/3;
Chris@0: tunings(6) = -3/4;
Chris@0: for k=1:parameter.numAdditionalTunings
Chris@0:     tunings(k+6) = -k/(parameter.numAdditionalTunings+1);
Chris@0: end
Chris@0: referenceFreqsA4 = 2.^((69-69+tunings)/12) * 440;
Chris@0: 
Chris@0: if inputIsAudioSignal
Chris@0:     [s,f,t] = spectrogram(f_input, parameter.windowFunction(parameter.fftWindowLength), parameter.fftWindowLength/2, parameter.fftWindowLength, 22050);
Chris@0: else
Chris@0:     s = f_input;
Chris@0:     f = parameter.freqinfo;
Chris@0: end
Chris@0: s = abs(s);
Chris@0: 
Chris@0: directFreqBinSearch = 0;
Chris@0: if all( (f(2:end)-f(1:end-1)) - (f(2)-f(1)) < eps )
Chris@0:     directFreqBinSearch = 1;
Chris@0: end
Chris@0: 
Chris@0: averagedPowerSpectrogram = sum(s.^2,2);
Chris@0: totalPitchEnergyViaSpec = zeros(numTunings,1);
Chris@0: for tu=1:numTunings
Chris@0:     centerfreqs = 2.^((parameter.pitchRange-69)/12) * referenceFreqsA4(tu);
Chris@0:     upperborderfreqs = 2.^((parameter.pitchRange-68.5)/12) * referenceFreqsA4(tu);
Chris@0:     lowerborderfreqs = 2.^((parameter.pitchRange-69.5)/12) * referenceFreqsA4(tu);
Chris@0:     
Chris@0:     % build triangular filterbank for magnitude spectrogram
Chris@0:     spectrogramFilter = zeros(length(f),1);
Chris@0:     for k=1:length(parameter.pitchRange)        
Chris@0:         c = getCorrespondingBin(f,centerfreqs(k),directFreqBinSearch);
Chris@0:         u = getCorrespondingBin(f,upperborderfreqs(k),directFreqBinSearch);
Chris@0:         l = getCorrespondingBin(f,lowerborderfreqs(k),directFreqBinSearch);
Chris@0:         
Chris@0:         % order is important here. If third parameter is < 2, then linspace
Chris@0:         % returns the second parameter
Chris@0:         spectrogramFilter(c:u) = parameter.pitchWeights(k) * linspace(1,0,u-c+1);
Chris@0:         spectrogramFilter(l:c) = parameter.pitchWeights(k) * linspace(0,1,c-l+1);
Chris@0:     end
Chris@0: 
Chris@0:     totalPitchEnergyViaSpec(tu) = sum(spectrogramFilter.^2 .* averagedPowerSpectrogram);    
Chris@0: end
Chris@0: 
Chris@0: [ignoreMe, maxIndex] = max(totalPitchEnergyViaSpec(1:6));
Chris@0: shiftFB = maxIndex-1;
Chris@0: 
Chris@0: [ignoreMe, maxIndex] = max(totalPitchEnergyViaSpec);
Chris@0: centerA4 = referenceFreqsA4(maxIndex);
Chris@0: tuningSemitones = tunings(maxIndex);
Chris@0: 
Chris@0: sideinfo.tuning.shiftFB = shiftFB;
Chris@0: sideinfo.tuning.centerA4 = centerA4;
Chris@0: sideinfo.tuning.tuningSemitones = tuningSemitones;
Chris@0: sideinfo.tuning.method = 'estimateTuningV1';
Chris@0: sideinfo.tuning.numAdditionalTunings = parameter.numAdditionalTunings;
Chris@0: sideinfo.tuning.pitchRange = parameter.pitchRange;
Chris@0: sideinfo.tuning.pitchWeights = parameter.pitchWeights;
Chris@0: sideinfo.tuning.fftWindowLength = parameter.fftWindowLength;
Chris@0: sideinfo.tuning.windowFunction = parameter.windowFunction;
Chris@0: sideinfo.tuning.inputWasAudioSignal = inputIsAudioSignal;
Chris@0: 
Chris@0: end
Chris@0: 
Chris@0: function index = getCorrespondingBin(x,sval,directSearch)
Chris@0: % - Finds the entry in x with the smallest absolute distance to sval.
Chris@0: % - x is assumed to be sorted (ascending)
Chris@0: % - 'directSearch' means that all values in x are equally spaced
Chris@0: % - x is assumed to be at least of length 2.
Chris@0: % - If directSearch==0 then we use binary seach to find the entry
Chris@0: %
Chris@0: % You can test the correctness of this procedure by comparing it against
Chris@0: % the result of:  [ignoreMe index] = min(abs(x-sval))
Chris@0: %
Chris@0: % Author: Sebastian Ewert
Chris@0: 
Chris@0: if sval >= x(end)
Chris@0:     index = length(x);
Chris@0:     return;
Chris@0: elseif sval <= x(1)
Chris@0:     index = 1;
Chris@0:     return;
Chris@0: end
Chris@0: 
Chris@0: 
Chris@0: if directSearch
Chris@0:     index = round( (sval-x(1)) / (x(2)-x(1))) + 1;
Chris@0: else
Chris@0:     from=1;
Chris@0:     to=length(x);
Chris@0:     
Chris@0:     while from<=to
Chris@0:         mid = round((from + to)/2);
Chris@0:         diff = x(mid)-sval;
Chris@0:         if diff<0   % x(mid) < sval
Chris@0:             from=mid;
Chris@0:         else              % x(mid) => sval
Chris@0:             to=mid;
Chris@0:         end
Chris@0:         
Chris@0:         if to-from==1
Chris@0:             break;
Chris@0:         end
Chris@0:     end
Chris@0:     if abs(x(from)-sval) < abs(x(to)-sval)
Chris@0:         index = from;
Chris@0:     else
Chris@0:         index = to;
Chris@0:     end
Chris@0: end
Chris@0: 
Chris@0: end
Chris@0: