yading@10: function y = harmonicmodel(x, fs, w, N, t, nH, minf0, maxf0, f0et, maxhd)
yading@10: % Analysis/synthesis of a sound using the sinusoidal harmonic model 
yading@10: % x: input sound, fs: sampling rate, w: analysis window (odd size),  
yading@10: % N: FFT size (minimum 512), t: threshold in negative dB,  
yading@10: % nH: maximum number of harmonics, minf0: minimum f0 frequency in Hz,  
yading@10: % maxf0: maximim f0 frequency in Hz,  
yading@10: % f0et: error threshold in the f0 detection (ex: 5), % maxhd: max. relative deviation in harmonic detection (ex: .2) 
yading@10: % y: output sound 
yading@10: M = length(w);                           % analysis window size 
yading@10: Ns= 1024;                                % FFT size for synthesis 
yading@10: H = 256;                                 % hop size for analysis and synthesis 
yading@10: N2 = N/2+1;                              % size postive spectrum 
yading@10: soundlength = length(x);                 % length of input sound array 
yading@10: hNs = Ns/2;                              % half synthesis window size 
yading@10: hM = (M-1)/2;                            % half analysis window size 
yading@10: pin = max(hNs+1,1+hM);   % initialize sound pointer to middle of analysis window 
yading@10: pend = soundlength-hM;                   % last sample to start a frame 
yading@10: fftbuffer = zeros(N,1);                  % initialize buffer for FFT 
yading@10: y = zeros(soundlength+Ns/2,1);           % output sound 
yading@10: w = w/sum(w);                            % normalize analysis window 
yading@10: sw = zeros(Ns,1); 
yading@10: ow = triang(2*H-1);                      % overlapping window 
yading@10: ovidx = Ns/2+1-H+1:Ns/2+H;               % overlap indexes 
yading@10: sw(ovidx) = ow(1:2*H-1); 
yading@10: bh = blackmanharris(Ns);                 % synthesis window 
yading@10: bh = bh ./ sum(bh);                      % normalize synthesis window 
yading@10: sw(ovidx) = sw(ovidx) ./ bh(ovidx); 
yading@10: while pin<pend 
yading@10:   %-----analysis-----% 
yading@10:   xw = x(pin-hM:pin+hM).*w(1:M);         % window the input sound 
yading@10:   fftbuffer(:) = 0;                      % reset buffer 
yading@10:   fftbuffer(1:(M+1)/2) = xw((M+1)/2:M);  % zero-phase window in fftbuffer 
yading@10:   fftbuffer(N-(M-1)/2+1:N) = xw(1:(M-1)/2); 
yading@10:   X = fft(fftbuffer);                    % compute the FFT 
yading@10:   mX = 20*log10(abs(X(1:N2)));           % magnitude spectrum  
yading@10:   pX = unwrap(angle(X(1:N/2+1)));        % unwrapped phase spectrum  
yading@10:   ploc = 1 + find((mX(2:N2-1)>t) .* (mX(2:N2-1)>mX(3:N2)) ... 
yading@10:                   .* (mX(2:N2-1)>mX(1:N2-2)));            % find peaks 
yading@10:   [ploc,pmag,pphase] = peakinterp(mX,pX,ploc);            % refine peak values 
yading@10:   f0 = f0detection(mX,fs,ploc,pmag,f0et,minf0,maxf0);     % find f0 
yading@10:   hloc = zeros(nH,1);                        % initialize harmonic locations 
yading@10:   hmag = zeros(nH,1)-100;                    % initialize harmonic magnitudes 
yading@10:   hphase = zeros(nH,1);                      % initialize harmonic phases 
yading@10:   hf = (f0>0).*(f0.*(1:nH));                 % initialize harmonic frequencies 
yading@10:   hi = 1;                                    % initialize harmonic index 
yading@10:   npeaks = length(ploc);                     % number of peaks found 
yading@10:   while (f0>0 && hi<=nH && hf(hi)<fs/2)      % find harmonic peaks 
yading@10:     [dev,pei] = min(abs((ploc(1:npeaks)-1)/N*fs-hf(hi)));       % closest peak 
yading@10:     if ((hi==1 || ~any(hloc(1:hi-1)==ploc(pei))) && dev<maxhd*hf(hi)) 
yading@10:       hloc(hi) = ploc(pei);                  % harmonic locations 
yading@10:       hmag(hi) = pmag(pei);                  % harmonic magnitudes 
yading@10:       hphase(hi) = pphase(pei);              % harmonic phases 
yading@10:     end 
yading@10:     hi = hi+1;                               %increase harmonic index 
yading@10:   end 
yading@10:   hloc(1:hi-1) = (hloc(1:hi-1)~=0).*((hloc(1:hi-1)-1)*Ns/N+1); % synth. locs 
yading@10:   %-----synthesis-----% 
yading@10:   Yh = genspecsines(hloc(1:hi-1),hmag,hphase,Ns);       % generate sines 
yading@10:   yh = fftshift(real(ifft(Yh)));                        % sines in time domain  
yading@10:   y(pin-hNs:pin+hNs-1) = y(pin-hNs:pin+hNs-1) + sw.*yh(1:Ns);    % overlap-add 
yading@10:   pin = pin+H;                               % advance the input sound pointer 
yading@10: end