yading@10: 
yading@10: function [f0, f0error] = TWM (ploc, pmag, N, fs, minf0, maxf0) 
yading@10: % Two-way mismatch algorithm (by Beauchamp&Maher) 
yading@10: % ploc, pmag: peak locations and magnitudes, N: size of complex spectrum 
yading@10: % fs: sampling rate of sound, minf0: minimum f0, maxf0: maximum f0 
yading@10: % f0: fundamental frequency detected, f0error: error measure 
yading@10: pfreq = (ploc-1)/N*fs;                                 % frequency in Hertz of peaks 
yading@10: [zvalue,zindex] = min(pfreq); 
yading@10: if (zvalue==0)                                   % avoid zero frequency peak 
yading@10:   pfreq(zindex) = 1; 
yading@10:   pmag(zindex) = -100; 
yading@10: end 
yading@10: ival2 = pmag; 
yading@10: [Mmag1,Mloc1] = max(ival2);                      % find peak with maximum magnitude 
yading@10: ival2(Mloc1) = -100;                             % clear max peak 
yading@10: [Mmag2,Mloc2]= max(ival2);                       % find second maximum magnitude peak 
yading@10: ival2(Mloc2) = -100;                             % clear second max peak 
yading@10: [Mmag3,Mloc3]= max(ival2);                       % find third maximum magnitude peak 
yading@10: nCand = 3;                     % number of possible f0 candidates for each max peak 
yading@10: f0c = zeros(1,3*nCand);                          % initialize array of candidates 
yading@10: f0c(1:nCand)=(pfreq(Mloc1)*ones(1,nCand))./((nCand+1-(1:nCand))); % candidates 
yading@10: f0c(nCand+1:nCand*2)=(pfreq(Mloc2)*ones(1,nCand))./((nCand+1-(1:nCand))); 
yading@10: f0c(nCand*2+1:nCand*3)=(pfreq(Mloc3)*ones(1,nCand))./((nCand+1-(1:nCand))); 
yading@10: f0c = f0c((f0c<maxf0)&(f0c>minf0));              % candidates within boundaries 
yading@10: if (isempty(f0c))                                % if no candidates exit 
yading@10:   f0 = 0; f0error=100; 
yading@10:   return 
yading@10: end 
yading@10: harmonic = f0c; 
yading@10: ErrorPM = zeros(fliplr(size(harmonic)));                   % initialize PM errors 
yading@10: MaxNPM = min(10,length(ploc)); 
yading@10: for i=1:MaxNPM % predicted to measured mismatch error 
yading@10:   difmatrixPM = harmonic' * ones(size(pfreq))'; 
yading@10:   difmatrixPM = abs(difmatrixPM-ones(fliplr(size(harmonic)))*pfreq'); 
yading@10:   [FreqDistance,peakloc] = min(difmatrixPM,[],2); 
yading@10:   Ponddif = FreqDistance .* (harmonic'.^(-0.5)); 
yading@10:   PeakMag = pmag(peakloc); 
yading@10:   MagFactor = 10.^((PeakMag-Mmag1)./20); 
yading@10:   ErrorPM = ErrorPM+(Ponddif+MagFactor.*(1.4*Ponddif-0.5)); 
yading@10:   harmonic = harmonic+f0c; 
yading@10: end 
yading@10: ErrorMP = zeros(fliplr(size(harmonic)));                   % initialize MP errors 
yading@10: MaxNMP = min(10,length(pfreq)); 
yading@10: for i=1:length(f0c)                         % measured to predicted mismatch error 
yading@10:   nharm = round(pfreq(1:MaxNMP)/f0c(i)); 
yading@10:   nharm = (nharm>=1).*nharm + (nharm<1); 
yading@10:   FreqDistance = abs(pfreq(1:MaxNMP) - nharm*f0c(i)); 
yading@10:   Ponddif = FreqDistance.* (pfreq(1:MaxNMP).^(-0.5)); 
yading@10:   PeakMag = pmag(1:MaxNMP); 
yading@10:   MagFactor = 10.^((PeakMag-Mmag1)./20); 
yading@10:   ErrorMP(i) = sum(MagFactor.*(Ponddif+MagFactor.*(1.4*Ponddif-0.5))); 
yading@10: end 
yading@10: Error = (ErrorPM/MaxNPM) + (0.3*ErrorMP/MaxNMP);           % total errors 
yading@10: [f0error, f0index] = min(Error);                           % get the smallest error 
yading@10: f0 = f0c(f0index);                                         % f0 with the smallest error