yading@10: function X = stft_peak(x, w, N, H, threshold)
yading@10: % Analysis/synthesis of a sound using the short-time fourier transform
yading@10: % x: input sound,
yading@10: % w: analysis window (odd number of samples required)
yading@10: % N: FFT size,
yading@10: % H: hop size
yading@10: % y: output sound
yading@10: lengthx = length(x);        % length of sound signal
yading@10: w = w/sum(w);               % normalize window
yading@10: W = length(w);              % window size (odd)
yading@10: W2 = (W-1)/2;               % half window
yading@10: N2 = N/2+1;                 % half FFT size 
yading@10: 
yading@10: p_start = 1+W2;             % first frame pointer
yading@10: p_end = lengthx-W2;         % last frame pointer
yading@10: fft_buffer = zeros(N,1);    % FFT buffer
yading@10: for g=p_start:H:p_end      % advance hop size each loop
yading@10:     % ANALYSIS
yading@10:     xw = w(1:W)'.*x(g-W2:g+W2);         % multiply window by input frame
yading@10:     fft_buffer(1:W2+1) = xw(W2+1:W);    % adjust buffer for zero-phase computation
yading@10:     fft_buffer(N-W2+1:N) = xw(1:W2);
yading@10:     X = fft(fft_buffer);                % FFT of current buffer
yading@10:     mX = 20*log10(abs(X(1:N2)));        % magnitude spectrum
yading@10:     pX = unwrap(angle(X(1:N2)));        % phase spectrum (unwrapped)
yading@10:     ploc = 1+find((mX(2:N2-1)>threshold).*(mX(2:N2-1)>mX(3:N2)).*(mX(2:N2-1)>mX(1:N2-2)));
yading@10:     pmag = mX(ploc);
yading@10:     pphase = pX(ploc);
yading@10:     plength = length(ploc);
yading@10:     if(mod((g-p_start)/H+1,(floor(p_end/H/5)))==0)    % select only 5 equally distributed frames of the sound
yading@10:         figure
yading@10:         subplot(2,1,1)
yading@10:         plot(mX)                                      % plot magnitude and phase spectrums
yading@10:         hold on
yading@10:         plot(ploc(1:plength),pmag(1:plength),'xr');   % add peaks to the plot  
yading@10:         xlim([1 N2])                                  
yading@10:         ylim([-80 0])
yading@10:         title('Magnitude spectrum of frame')
yading@10:         ylabel('Magnitude spectrum (dB)')
yading@10:         xlabel('Frequency (rad)')
yading@10:         display(['Peak positions: ' num2str(ploc')]);
yading@10:         hold off
yading@10:         subplot(2,1,2)
yading@10:         plot(pX)                                      % same with the phase
yading@10:         hold on
yading@10:         plot(ploc(1:plength),pphase(1:plength),'xr');
yading@10:         xlim([1 N2])                     
yading@10: %         ylim([-pi pi])
yading@10: %         pause()
yading@10:         hold off
yading@10:     end
yading@10: end