yading@10: function [mY, pY] = stft(x, w, N, H)
yading@10: % Analysis/synthesis of a sound using the short-time fourier transform
yading@10: % x: input sound, w: analysis window (odd size), N: FFT size, H: hop size
yading@10: % 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 in 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 or reset buffer
yading@10: w = w/sum(w); % normalize analysis window
yading@10: nframes=mod(soundlength-hM,H);
yading@10: frame=1;
yading@10: %fftbuffer=zeros(N,1);
yading@10: mY=zeros(N2,nframes);
yading@10: pY=zeros(N2,nframes);
yading@10: 
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 window in fftbuffer
yading@10: fftbuffer(N-(M-1)/2+1:N) = xw(1:(M-1)/2);
yading@10: X = fft(fftbuffer); % compute 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 spect. of positive freq.
yading@10: mY(:,frame)=mX;
yading@10: pY(:,frame)=pX;
yading@10: frame=frame+1;
yading@10: pin = pin+H; % advance sound pointer
yading@10: end
yading@10: % Copied from Carles
yading@10: mY_norm=mY-(mY./mY).*max(max(mY)); 
yading@10: %Spectrogram mode respresentation
yading@10: surface(1:frame-1, linspace(0,pi,N2),mY_norm(:,1:frame-1),'LineStyle','none');
yading@10: %Axes adjustment
yading@10: xlim([1 frame-1]);
yading@10: ylim([0 pi]);
yading@10: 
yading@10: 
yading@10: 
yading@10: [x,fs,nbit]=wavread('Ocarina.wav');