katkost@3: function y = stpt(x, w, N, H, t)
katkost@3: % Analysis/synthesis of a sound using the peaks
katkost@3: % of the short-time fourier transform
katkost@3: % x: input sound, 
katkost@3: % w: analysis window (odd size), 
katkost@3: % N: FFT size, 
katkost@3: % H: hop size,
katkost@3: % t: threshold in negative dB, 
katkost@3: % y: output sound
katkost@3: M = length(w); % analysis window size
katkost@3: N2 = N/2+1; % size of positive spectrum
katkost@3: soundlength = length(x); % length of input sound array
katkost@3: hM = (M-1)/2; % half analysis window size
katkost@3: pin = 1+hM; % initialize sound pointer at the middle of analysis window
katkost@3: pend = soundlength-hM; % last sample to start a frame
katkost@3: fftbuffer = zeros(N,1); % initialize buffer for FFT
katkost@3: yw = zeros(M,1); % initialize output sound frame
katkost@3: y = zeros(soundlength,1); % initialize output array
katkost@3: w = w/sum(w); % normalize analysis window
katkost@3: sw = hanning(M); % synthesis window
katkost@3: sw = sw./sum(sw);
katkost@3: while pin<pend
katkost@3: %-----analysis-----%
katkost@3: xw = x(pin-hM:pin+hM).*w(1:M); % window the input sound
katkost@3: fftbuffer(:) = 0; % reset buffer
katkost@3: fftbuffer(1:(M+1)/2) = xw((M+1)/2:M); % zero-phase fftbuffer
katkost@3: fftbuffer(N-(M-1)/2+1:N) = xw(1:(M-1)/2);
katkost@3: X = fft(fftbuffer); % compute the FFT
katkost@3: mX = 20*log10(abs(X(1:N2))); % magnitude spectrum of positive frequencies
katkost@3: pX = unwrap(angle(X(1:N2))); % unwrapped phase spectrum
katkost@3: ploc = 1 + find((mX(2:N2-1)>t) .* (mX(2:N2-1)>mX(3:N2)) ...
katkost@3: .* (mX(2:N2-1)>mX(1:N2-2))); % peakss
katkost@3: pmag = mX(ploc); % magnitude of peaks
katkost@3: pphase = pX(ploc); % phase of peaks
katkost@3: %-----synthesis-----%
katkost@3: Y = zeros(N,1); % initialize output spectrum
katkost@3: Y(ploc) = 10.^(pmag/20).*exp(i.*pphase); % generate positive freq.
katkost@3: Y(N+2-ploc) = 10.^(pmag/20).*exp(-i.*pphase); % generate negative freq.
katkost@3: fftbuffer = real(ifft(Y)); % real part of the inverse FFT
katkost@3: yw((M+1)/2:M) = fftbuffer(1:(M+1)/2); % undo zero phase window
katkost@3: yw(1:(M-1)/2) = fftbuffer(N-(M-1)/2+1:N);
katkost@3: y(pin-hM:pin+hM) = y(pin-hM:pin+hM) + H*N*sw.*yw(1:M); % overlap-add
katkost@3: pin = pin+H; % advance sound pointer
katkost@3: end