%==============================================================
% Harmonic Model
%==============================================================
function y = harmonicmodel(x, fs, w, N, H, mthreshold, nHarms, minf0, maxf0, maxf0error, maxharmdev)
% x: input array sound
% fs: sampling rate
% w: window (odd size), ex: blackmanharris(1023)
% N: FFT size, ex: 1024
% H: hop size, ex: 512
% mthreshold: threshold in dB used for peak detection, ex: -20
% nHarms: maximum number of harmonics to find, ex: 40
% minf0: minimum frequency for f0, ex: 80
% maxf0: maximim frequency for f0, ex: 500
% maxf0error: maximum error allowed in the f0 detection, ex: 5
% maxharmdev: maximum deviation allowed in detecting harmonics, as a
% product, ex: .2
M = length(w);   % window size
N2 = N/2; % half of FFT size
soundlength = length(x);
pin = 0; % pointer into input sound array
l = 1; % frame index
pend = soundlength - M;
f0 = zeros(1,1+ceil(pend/H));
f0error = zeros(1,1+ceil(pend/H));
sw = 2*M/N * triang(M)./blackmanharris(M); % synthesis window
y = zeros(soundlength,1);  % output sound array
fftbuffer = zeros (N,1);
bh92lobe = genbh92lobe(); % generate a spectral lobe of a sinusoid
while pin<pend
    x1 = x(pin+1:pin+M) .* w(1:M);   % window a frame of the input sound
    fftbuffer(1:((M+1)/2)) = x1((M+1)/2:M); % zero-phase window
    fftbuffer(N-(M-1)/2+1:N) = x1(1:(M-1)/2);
    X = fft(fftbuffer); % compute the FFT
    mX = 20*log10(abs(X(1:N2)));  % magnitude spectrum in dB
    pX = angle(X(1:N2)); % phase spectrum
    specder = diff([mthreshold; mX(1:N2);mthreshold]);  % calculate magnitude differences
    ploc = find(mX(1:N2) >=mthreshold & specder(1:N2)>= 0 & specder(2:N2+1) <= 0); % find local maxima
    [iploc, ipmag, ipphase] = peakinterp (mX, pX, ploc); % refine peak values
    [f0t,f0errort] = f0detection(mX,fs,iploc,ipmag); % find fundamental frequency
    f0(l) = f0t * (f0errort < maxf0error) * (f0t > minf0) * (f0t < maxf0);
    f0error(l) = f0errort;
    ihloc = zeros(nHarms,1);
    ihmag = zeros(nHarms,1) -100;
    ihphase = zeros(nHarms,1);
    if (f0(l)>0)
        for i=1:nHarms
            trajfreq = f0(l)*i;
            if (trajfreq < fs/2)
                [closestpeakmag,closestpeakindex]=min(abs((iploc-1)/N*fs-trajfreq));
                if (abs(iploc(closestpeakindex)/N*fs - trajfreq) < maxharmdev*trajfreq)
                ihloc(i) = iploc(closestpeakindex);
                ihmag(i) = ipmag(closestpeakindex);
                ihphase(i) = ipphase(closestpeakindex);
                end
            end
        end
    end
    X = genspecsines(ihloc,ihmag,ihphase,M,N,bh92lobe); % generate spectral sines
    fftbuffer = real(ifft(X)); % inverse FFT
    x1((M+1)/2:M) = fftbuffer(1:(M+1)/2);
    x1(1:(M-1)/2) = fftbuffer(N-(M-1)/2+1:N);
    y(pin+1:pin+M) = y(pin+1:pin+M) + x1(1:M) .* sw; % overlap-add
    pin  = pin + H;
    l = l + 1;
end