yading@10: function [y]=harmonicmodel1(x,w,N,t)
yading@10: %initializing values
yading@10: M = length(w); % window size - the longer the more frequency resolution
yading@10: N2 = N/2+1; % positive part of the spectrum
yading@10: Ns= 1024; % FFT size for synthesis (even)
yading@10: H = 256; % analysis/synthesishop size
yading@10: soundlength = length(x); % length of input sound array - samples
yading@10:  
yading@10: fftbuffer = zeros(N,1); % initialize buffer for FFT
yading@10:  
yading@10: %Create a loop to step through the sound array x
yading@10: %initializing the loop
yading@10: hNs = Ns/2; % half synthesis window size
yading@10: hM = (M-1)/2; % half analysis window size used to overlap windows
yading@10:  
yading@10: pin = max(hNs+1,1+hM); % initialize sound pointer to middle of analysis window
yading@10: pend = soundlength-max(H,hM); % last sample to start a frame
yading@10:  
yading@10: y = zeros(soundlength,1); % initialize output array
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-hNs+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:  
yading@10: while pin<pend
yading@10:     xw = x(pin-hM:pin+hM).*w(1:M)'; % window the input sound - STFT definition
yading@10:   
yading@10:     %zero phased window
yading@10:     fftbuffer(:) = 0; % reset buffer
yading@10:     fftbuffer(1:(M+1)/2) = xw((M+1)/2:M); % zero-phase fftbuffer
yading@10:     fftbuffer(N-(M-1)/2+1:N) = xw(1:(M-1)/2);
yading@10:  
yading@10:     %compute FFT of the zero phased frame
yading@10:     X = fft(fftbuffer);
yading@10:  
yading@10:     %calculate  magnitude and phase spectrum of of positive frequencies
yading@10:     mX = 20*log10(abs(X(1:N2)));
yading@10:     pX = unwrap(angle(X(1:N2)));
yading@10:     
yading@10:         
yading@10:     %Find the locations, ploc, of the local maxima above a given 
yading@10:     %threshold, t, in each magnitude spectrum by finding changes of slope.
yading@10:     ploc = 1+find((mX(2:N2-1)>t).*(mX(2:N2-1)>mX(3:N2)).*(mX(2:N2-1)>mX(1:N2-2)));
yading@10:     
yading@10:     %Find the magnitudes, pmag, and phases, pphase, of the obtained
yading@10:     %locations.
yading@10:     pmag = mX(ploc);
yading@10:     %pmag = mX(ploc)*0.4; 
yading@10:     pphase = pX(ploc); 
yading@10:     
yading@10:     
yading@10:     %peak interpolation
yading@10:     [iploc, ipmag, ipphase] = peakinterp (mX, pX, ploc);
yading@10:     
yading@10:     
yading@10:     %plot for a window
yading@10:     if (pin==1+hM+20*H)
yading@10:         figure
yading@10:         subplot(2,1,1)
yading@10:         plot(mX)
yading@10:         hold on
yading@10:         plot(iploc,mX(ploc),'*')
yading@10:         title('magnitude and peak values');
yading@10:         hold off
yading@10:         subplot(2,1,2)
yading@10:         plot(pX) 
yading@10:         hold on
yading@10:         plot(iploc,pX(ploc),'*')
yading@10:         title('phase and peak values');
yading@10:         hold off
yading@10:         
yading@10:         %number of peaks
yading@10:         npeaksm = length(pmag)
yading@10:         npeaksp = length(pphase)
yading@10:     end
yading@10:     
yading@10:             
yading@10:    %-----synthesis-----%
yading@10: plocs = (ploc-1)*Ns/N+1; % adapt peak locations to synthesis FFT
yading@10: Y = genspecsines(plocs,pmag,pphase,Ns); % generate spec sines
yading@10: yw = fftshift(real(ifft(Y))); % time domain of sinusoids
yading@10: y(pin-hNs:pin+hNs-1) = y(pin-hNs:pin+hNs-1)+ sw.*yw(1:Ns); % overlap-add
yading@10: pin = pin+H; % advance the sound pointer
yading@10:  
yading@10: end
yading@10: end