Harmonic Model

function y = harmonicmodel(x, w, N, t, fs)
% Authors: J. Bonada, X. Serra, X. Amatriain, A. Loscos
% Analysis/synthesis of a sound using the sinusoidal model
% x: input sound, w: analysis window (odd size), N: FFT size,
% t: threshold in negative dB, y: output sound
%
%--------------------------------------------------------------------------
% This source code is provided without any warranties as published in
% DAFX book 2nd edition, copyright Wiley & Sons 2011, available at
% http://www.dafx.de. It may be used for educational purposes and not
% for commercial applications without further permission.
%--------------------------------------------------------------------------

M = length(w);                                % analysis window size
Ns= 1024;                                     % FFT size for synthesis (even)
H = 256;                                      % analysis/synthesishop size
N2= N/2+1;                                    % size of positive spectrum
soundlength = length(x);                      % length of input sound array
hNs = Ns/2;                                   % half synthesis window size
hM = (M-1)/2;                                 % half analysis window size
pin = max(H+1,1+hM); % initialize sound pointer to middle of analysis window
pend = soundlength-max(H,hM);                 % last sample to start a frame
fftbuffer = zeros(N,1);                       % initialize buffer for FFT
y = zeros(soundlength,1);	                    % initialize output array
w = w/sum(w);                                 % normalize analysis window
sw = zeros(Ns,1);
ow = triang(2*H-1);                           % overlapping window
ovidx = Ns/2+1-H+1:Ns/2+H;                    % overlap indexes
sw(ovidx) = ow(1:2*H-1);
bh = blackmanharris(Ns);                      % synthesis window
bh = bh ./ sum(bh);                           % normalize synthesis window
sw(ovidx) = sw(ovidx) ./ bh(ovidx);
k = [0:1:N];
freq = fs*k/N;
minf0 = 100;
maxf0 = 20000;
ef0max = 100;
while pin<pend
  %-----analysis-----%
  xw = x(pin-hM:pin+hM).*w(1:M);              % window the input sound
  fftbuffer(:) = 0;                           % reset buffer
  fftbuffer(1:(M+1)/2) = xw((M+1)/2:M);       % zero-phase window in fftbuffer
  fftbuffer(N-(M-1)/2+1:N) = xw(1:(M-1)/2);
  X = fft(fftbuffer);                              % compute the FFT
  mX = 20*log10(abs(X(1:N2)));    % magnitude spectrum of positive frequencies
  pX = unwrap(angle(X(1:N/2+1)));                  % unwrapped phase spectrum
  ploc = 1 + find((mX(2:N2-1)>t) .* (mX(2:N2-1)>mX(3:N2)) ...
                  .* (mX(2:N2-1)>mX(1:N2-2)));     % find peaks
  [ploc,pmag,pphase] = peakinterp(mX,pX,ploc);     % refine peak values
%   f0 = f0detection(mX, fs, ploc, pmag, ef0max, minf0, maxf0);
  f0 = f0detectionyin(x, fs, 1200, minf0, maxf0) ;
  nH=20;
  maxhd=.2;
  hf = (f0>0).*(f0.*(1:nH));                 % initialize harmonic frequencies
  hi = 1;                                    % initialize harmonic index
  npeaks = length(ploc);                     % number of peaks found
  while (f0>0 && hi<=nH && hf(hi)<fs/2)      % find harmonic peaks
      [dev,pei] = min(abs((ploc(1:npeaks)-1)/N*fs-hf(hi)));       % closest peak
      if ((hi==1 || ~any(hloc(1:hi-1)==ploc(pei))) && dev<maxhd*hf(hi))
        hloc(hi) = ploc(pei);                  % harmonic locations
        hmag(hi) = pmag(pei);                  % harmonic magnitudes
        hphase(hi) = pphase(pei);              % harmonic phases
      end
      hi = hi+1;                               %increase harmonic index
  end
%   f0
%   f0yin
%   subplot(2,1,1);
%   plot(freq(1:N2),mX);
%   hold on;
%   plot(f0,mX(round(f0*N/fs)),'g*');
%   plot(fs*hloc/N,hmag,'rX');

%   plot(f0yin,mX(round(f0yin*N/fs)),'g*');
%   hold off;
%   pause;
%   subplot(2,1,2);
%   plot(f0,mX(round(f0*N/fs)),'X');
  %-----synthesis-----%
  plocs = (ploc-1)*Ns/N;             % adapt peak locations to synthesis FFT
  Y = genspecsines(plocs,pmag,pphase,Ns);          % generate spec sines
  yw = fftshift(real(ifft(Y)));                    % time domain of sinusoids
  %(pin-hNs:pin+hNs-1)
  y(pin-hNs:pin+hNs-1) = y(pin-hNs:pin+hNs-1) + sw.*yw(1:Ns);    % overlap-add
  pin  = pin+H;                                    % advance the sound pointer
end