%==============================================================
% Sinusoidal Model
%==============================================================
function y = sinemodel(x, w, N, H, t)
% x: input sound array
% w: window (should be odd size)
% N: FFT size
% H: hop size
% t: threshold in dB
M = length(w);   % window size
N2 = N/2; % half of FFT size
soundlength = length(x);
pin = 0;
pend = soundlength - M;
ws = .5* N/M * 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:N)));  % magnitude spectrum in dB
   pX = angle(X(1:N)); % phase spectrum
   specder = diff([t; mX(1:N2); t]);  % calculate magnitude differences
   ploc = find(mX(1:N2) >= t & specder(1:N2)>= 0 & specder(2:N2+1) <= 0); % find local maxima
   [iploc, ipmag, ipphase] = peakinterp (mX, pX, ploc);
   X = genspecsines(iploc,ipmag,ipphase,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).*ws; % overlap-add
   pin  = pin + H;  
end