yading@10: function y = stftlab4(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 size,
yading@10: % t: threshold in negative dB, y: output sound
yading@10: M = length(w); % analysis window size
yading@10: N2 = N/2+1; % size of positive spectrum
yading@10: soundlength = length(x); % length of input sound array
yading@10: hM = (M-1)/2; % half analysis window size
yading@10: pin = 1+hM; % initialize sound pointer at the middle of analysis window
yading@10: pend = soundlength-hM; % last sample to start a frame
yading@10: fftbuffer = zeros(N,1); % initialize buffer for FFT
yading@10: yw = zeros(M,1); % initialize output sound frame
yading@10: y = zeros(soundlength,1); % initialize output array
yading@10: w = w/sum(w); % normalize analysis window
yading@10: sw = hanning(M); % synthesis window
yading@10: sw = sw./sum(sw);
yading@10: while pin<pend
yading@10: %-----analysis-----%
yading@10: xw = x(pin-hM:pin+hM).*w(1:M)'; % window the input sound
yading@10: fftbuffer(:) = 0; % reset buffer
yading@10: fftbuffer(1:(M+1)/2) = xw((M+1)/2:M); % zero-phase fftbuffer
yading@10: fftbuffer(N-(M-1)/2+1:N) = xw(1:(M-1)/2);
yading@10: X = fft(fftbuffer); % compute the FFT
yading@10: mX = 20*log10(abs(X(1:N2))); % magnitude spectrum of positive frequencies
yading@10: pX = unwrap(angle(X(1:N2))); % unwrapped phase spectrum
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))); % peaks
yading@10: [nploc,npmag,npphase] = peakinterp(mX,pX,ploc); % refine peak values. i use new parameters
yading@10: %so that i can plot the new values at the figures existed
yading@10: pmag = mX(ploc); % finding the magnitude of the obtained locations (peaks)
yading@10: pphase = pX(ploc); % finding the phase of the obtained locations (peaks)
yading@10:  
yading@10: % plotting the peak values on top of the magnitude and phase spectra at the
yading@10: % magnitude and phase frame
yading@10:   
yading@10: subplot(211); 
yading@10: plot(mX);
yading@10: hold on
yading@10: plot(ploc,pmag,'m*');   %m is for magenta color
yading@10: plot(nploc,npmag,'c*'); %plot the new magnitude spectra
yading@10: title('magnitude and peak values');
yading@10: hold off
yading@10: subplot(212);
yading@10: plot(pX);
yading@10: hold on
yading@10: plot(ploc,pphase,'m*');
yading@10: plot(nploc,npphase,'c*');
yading@10: title('phase and peak values');
yading@10: hold off
yading@10:  
yading@10: end