katkost@1: function [y,yh,ys,fr0] = hpsmodel(x,fs,w,N,t,nH,minf0,maxf0,f0et,maxhd,stocf) 
katkost@1: %=> analysis/synthesis of a sound using the sinusoidal harmonic model 
katkost@1: % x: input sound, fs: sampling rate, w: analysis window (odd size),  
katkost@1: % N: FFT size (minimum 512), t: threshold in negative dB,  
katkost@1: % nH: maximum number of harmonics, minf0: minimum f0 frequency in Hz,  
katkost@1: % maxf0: maximim f0 frequency in Hz,  
katkost@1: % f0et: error threshold in the f0 detection (ex: 5), 
katkost@1: % maxhd: max. relative deviation in harmonic detection (ex: .2) 
katkost@1: % stocf: decimation factor of mag spectrum for stochastic analysis 
katkost@1: % y: output sound, yh: harmonic component, ys: stochastic component 
katkost@1: M = length(w);   % analysis window size 
katkost@1: Ns = 1024;                               % FFT size for synthesis 
katkost@1: H = 256;                                 % hop size for analysis and synthesis 
katkost@1: N2 = N/2+1;                              % half-size of spectrum 
katkost@1: soundlength = length(x);                 % length of input sound array 
katkost@1: hNs = Ns/2;                              % half synthesis window size 
katkost@1: hM = (M-1)/2;                            % half analysis window size 
katkost@1: pin = max(hNs+1,1+hM);   % initialize sound pointer to middle of analysis window 
katkost@1: pend = soundlength-max(hM,hNs);            % last sample to start a frame 
katkost@1: fftbuffer = zeros(N,1);                  % initialize buffer for FFT 
katkost@1: yh = zeros(soundlength+Ns/2,1);          % output sine component 
katkost@1: ys = zeros(soundlength+Ns/2,1);          % output residual component 
katkost@1: w = w/sum(w);                            % normalize analysis window 
katkost@1: sw = zeros(Ns,1); 
katkost@1: ow = triang(2*H-1);                      % overlapping window 
katkost@1: ovidx = Ns/2+1-H+1:Ns/2+H;               % overlap indexes 
katkost@1: sw(ovidx) = ow(1:2*H-1); 
katkost@1: bh = blackmanharris(Ns);                 % synthesis window 
katkost@1: bh = bh ./ sum(bh);                      % normalize synthesis window 
katkost@1: wr = bh;                                 % window for residual  
katkost@1: sw(ovidx) = sw(ovidx) ./ bh(ovidx); 
katkost@1: sws = H*hanning(Ns);                     % synthesis window for stochastic 
katkost@1: 
katkost@1: i = 0;
katkost@1: while pin<pend 
katkost@1:   i = i+1;
katkost@1:   %-----analysis-----% 
katkost@1:   xw = x(pin-hM:pin+hM).*w(1:M);         % window the input sound 
katkost@1:   fftbuffer(1:(M+1)/2) = xw((M+1)/2:M);  % zero-phase window in fftbuffer 
katkost@1:   fftbuffer(N-(M-1)/2+1:N) = xw(1:(M-1)/2); 
katkost@1:   X = fft(fftbuffer);                    % compute the FFT 
katkost@1:   mX = 20*log10(abs(X(1:N2)));           % magnitude spectrum  
katkost@1:   pX = unwrap(angle(X(1:N/2+1)));        % unwrapped phase spectrum  
katkost@1:   ploc = 1 + find((mX(2:N2-1)>t) .* (mX(2:N2-1)>mX(3:N2)) ... 
katkost@1:                   .* (mX(2:N2-1)>mX(1:N2-2)));          % find peaks 
katkost@1:   [ploc,pmag,pphase] = peakinterp(mX,pX,ploc);          % refine peak values 
katkost@1:   f0 = f0detection(mX,fs,ploc,pmag,f0et,minf0,maxf0);   % find f0 
katkost@1:   fr0(i)=f0;
katkost@1:   hloc = zeros(nH,1);                    % initialize harmonic locations 
katkost@1:   hmag = zeros(nH,1)-100;                % initialize harmonic magnitudes 
katkost@1:   hphase = zeros(nH,1);                  % initialize harmonic phases 
katkost@1:   hf = (f0>0).*(f0.*(1:nH));             % initialize harmonic frequencies 
katkost@1:   hi = 1;                                % initialize harmonic index 
katkost@1:   npeaks = length(ploc);                 % number of peaks found 
katkost@1:   while (f0>0 && hi<=nH && hf(hi)<fs/2)  % find harmonic peaks 
katkost@1:     [dev,pei] = min(abs((ploc(1:npeaks)-1)/N*fs-hf(hi)));    % closest peak 
katkost@1:     if ((hi==1 || ~any(hloc(1:hi-1)==ploc(pei))) && dev<maxhd*hf(hi)) 
katkost@1:       hloc(hi) = ploc(pei);              % harmonic locations 
katkost@1:       hmag(hi) = pmag(pei);              % harmonic magnitudes 
katkost@1:       hphase(hi) = pphase(pei);          % harmonic phases 
katkost@1:     end 
katkost@1:     hi = hi+1;                           % increase harmonic index 
katkost@1:   end 
katkost@1:   hloc(1:hi-1) = (hloc(1:hi-1)~=0).*((hloc(1:hi-1)-1)*Ns/N+1);  % synth. locs 
katkost@1:   ri= pin-hNs;                     % input sound pointer for residual analysis 
katkost@1:   xr = x(ri:ri+Ns-1).*wr(1:Ns);          % window the input sound 
katkost@1:   Xr = fft(fftshift(xr));                % compute FFT for residual analysis 
katkost@1:   Yh = genspecsines(hloc(1:hi-1),hmag,hphase,Ns);             % generate sines 
katkost@1:   Yr = Xr-Yh;                            % get the residual complex spectrum 
katkost@1:   mYr = abs(Yr(1:Ns/2+1));               % magnitude spectrum of residual 
katkost@1:   %mYs = stochenvelope(mYr,stocf);
katkost@1:   %-----transformations-----% 
katkost@1:   mYsenv = decimate(mYr,stocf,1); % decimate the magnitude spectrum
katkost@1:   %-----synthesis-----%
katkost@1:   mYs = interp(mYsenv,stocf,1); % interpolate to original size
katkost@1:   
katkost@1: % n=1:N/2+1; 
katkost@1: % plot(n/N*Ns,mX); %plotting the original spectrum
katkost@1: % hold on;
katkost@1: % plot(20*log10(abs(mYs)), 'r'); %plotting the approximation done by the decimate function
katkost@1: % 
katkost@1: % hold on;
katkost@1: % plot(hloc, hmag, 'g*');
katkost@1: % hold off;
katkost@1: % pause
katkost@1:   
katkost@1:   roffset = ceil(stocf/2)-1; % interpolated array offset
katkost@1:   mYs = [ mYs(1)*ones(roffset,1); mYs(1:Ns/2+1-roffset) ];
katkost@1:   pYs = 2*pi*rand(Ns/2+1,1);                 % generate phase random values 
katkost@1:   mYs1 = [mYs(1:Ns/2+1); mYs(Ns/2:-1:2)];    % create magnitude spectrum 
katkost@1:   pYs1 = [pYs(1:Ns/2+1); -1*pYs(Ns/2:-1:2)]; % create phase spectrum 
katkost@1:   Ys = mYs1.*cos(pYs1)+1i*mYs1.*sin(pYs1);   % compute complex spectrum 
katkost@1:   yhw = fftshift(real(ifft(Yh)));            % sines in time domain using IFFT 
katkost@1:   ysw = fftshift(real(ifft(Ys)));            % stoc. in time domain using IFFT 
katkost@1:   yh(ri:ri+Ns-1) = yh(ri:ri+Ns-1)+yhw(1:Ns).*sw;   % overlap-add for sines 
katkost@1:   ys(ri:ri+Ns-1) = ys(ri:ri+Ns-1)+ysw(1:Ns).*sws;  % overlap-add for stoch. 
katkost@1:   pin = pin+H;                                     % advance the sound pointer 
katkost@1: end 
katkost@1: 
katkost@1: %ys=tanh(10*ys);
katkost@1: y= yh+ys; % sum sines and stochastic