Stpt

function y = stpt(x, w, N, H, fs, t)
% Analysis/synthesis of a sound using the short-time Fourier transform
% Authors: J. Bonada, X. Serra, X. Amatriain, A. Loscos
% x: input sound, w: analysis window (odd size), N: FFT size, H: hop size
% y: output sound
%
%--------------------------------------------------------------------------
% This source code is provided without any warranties as published in
% DAFX book 2nd edition, copyright Wiley & Sons 2011, available at
% http://www.dafx.de. It may be used for educational purposes and not
% for commercial applications without further permission.
%--------------------------------------------------------------------------

M = length(w);                               % analysis window size
N2 = N/2+1;                                  % size of positive spectrum
soundlength = length(x);                     % length of input sound array
hM = (M-1)/2;                                % half analysis window size
pin = 1+hM;       % initialize sound pointer in middle of analysis window
pend = soundlength-hM;                       % last sample to start a frame
fftbuffer = zeros(N,1);                      % initialize buffer for FFT
yw = zeros(M,1);                             % initialize output sound frame
y = zeros(soundlength,1);                    % initialize output array
w = w/sum(w);                                % normalize analysis window
k = [0:1:N];
freq = fs*k/N;
%t = -80 ;                                    % threshold

while pin<pend
  %-----analysis-----%
  %clf;
  xw = x(pin-hM:pin+hM).*w(1:M);             % window the input sound

  fftbuffer(:) = 0;                          % reset buffer
  fftbuffer(1:(M+1)/2) = xw((M+1)/2:M);      % zero-phase window in fftbuffer
  fftbuffer(N-(M-1)/2+1:N) = xw(1:(M-1)/2);
  X = fft(fftbuffer);                        % compute FFT
  mX = 20*log10(abs(X(1:N2)));    % magnitude spectrum of positive frequencies
  ploc = 1 + find((mX(2:N2-1)>t).*(mX(2:N2-1)>mX(3:N2)).*(mX(2:N2-1)>mX(1:N2-2)));
  pX = unwrap(angle(X(1:N2)));    % unwrapped phase spect. of positive freq.
  pmag=mX(ploc);
  pphase=pX(ploc);
   %-----synthesis-----%
  Y = zeros(N,1);     % initialize output spectrum
  Y(ploc) = 10.^(pmag/20).*exp(i.*pphase); % generate positive freq.
  Y(N-ploc) = 10.^(pmag/20).*exp(-i.*pphase);% generate negative frequency


  fftbuffer = real(ifft(Y));                 % inverse FFT
  fftbuffer2 = fftbuffer.*hanning(N);        % smoothing with a window

  yw(1:(M-1)/2) = fftbuffer(N-(M-1)/2+1:N);  % undo zero-phase window
  yw((M+1)/2:M) = fftbuffer(1:(M+1)/2);

  y(pin-hM:pin+hM) = y(pin-hM:pin+hM) + H*yw(1:M);     % overlap-add
  pin  = pin+H;                              % advance sound pointer
  end