yading@10: function y=stpt2(x,w,N,H,t) 
yading@10: %Analysis/synthesis of a sound using the peaks 
yading@10: %of the short-time Fourier transform 
yading@10: %x: input sound, w: analysis window (odd size), N: FFT size, H: hop 
yading@10: %size, t: threshold in negative dB, y: output sound 
yading@10:   
yading@10: %Internal variables 
yading@10: M=length(w); 
yading@10: N2=N/2+1; 
yading@10: soundlength=length(x); 
yading@10: hM=(M-1)/2; 
yading@10:   
yading@10: %Pointers 
yading@10: pin=1+hM; 
yading@10: pend=soundlength-hM; 
yading@10:   
yading@10: %Initialize buffers 
yading@10: fftbuffer=zeros(N,1); 
yading@10: yw=zeros(M,1); 
yading@10: y=zeros(soundlength,1); 
yading@10:   
yading@10: %Analysis/Synthesis windows 
yading@10: sw=hanning(M); 
yading@10: sw=sw./sum(sw); 
yading@10: w=w./sum(w); 
yading@10:   
yading@10: while pin<pend 
yading@10:     %Analysis 
yading@10:     xw=x(pin-hM:pin+hM).*w(1:M); %Frame selector 
yading@10:     fftbuffer=zeros(N,1);        %Zero-phase windowing 
yading@10:     fftbuffer(1:(M+1)/2)=xw((M+1)/2:M); 
yading@10:     fftbuffer(N-(M-1)/2+1:N)=xw(1:(M-1)/2); 
yading@10:      
yading@10:     X=fft(fftbuffer);   %FFT performance 
yading@10:     mX=20*log10(abs(X(1:N2))); %Compute the FFT Magnitude (dB) 
yading@10:     pX=angle(X(1:N2));  %Compute the FFT phase 
yading@10:      
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)));      %Peak finder above threshold (t)    
yading@10:     pmag=mX(ploc);     %Peak magnitude recovery 
yading@10:     pphase=pX(ploc);   %Peak phase recovery 
yading@10:  
yading@10:     %Synthesis 
yading@10:     Y=zeros(N,1);    %Memory reserve 
yading@10:     Y(ploc)=10.^(pmag/2).*exp(1i.*pphase); %Spectrum reconstruction from peak info 
yading@10:     Y(N+2-ploc)=10.^(pmag/20).*exp(-1i.*pphase); %Spectrum reconstruction negative freq. 
yading@10:     fftbuffer=real(ifft(Y));     %Compute inverse FFT 
yading@10:     yw((M+1)/2:M)=fftbuffer(1:(M+1)/2); %Zero-phase windowing undo 
yading@10:     yw(1:(M-1)/2)=fftbuffer(N-(M-1)/2+1:N); 
yading@10:     y(pin-hM:pin+hM)=y(pin-hM:pin+hM)+H*N*sw.*yw(1:M); %overlap-add 
yading@10:     pin=pin+H;                   %Pointer advance 
yading@10: end 
yading@10:      
yading@10: