view extra/stptpeaks.m @ 13:844d341cf643 tip

Back up before ISMIR

function y=stptpeaks(x, w, N, H, t)
% Analysis/synthesis of a sound using the peaks
% of the short-time fourier transform
% x: input sound, 
% w: analysis window (odd size), 
% N: FFT size, 
% H: hop size,
% t: threshold in negative dB, 
% y: output sound

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 at the 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



while pin<pend
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 fftbuffer
fftbuffer(N-(M-1)/2+1:N) = xw(1:(M-1)/2);

X = fft(fftbuffer); % compute the FFT
mX = 20*log10(abs(X(1:N2))); % magnitude spectrum of positive frequencies
pX = unwrap(angle(X(1:N2))); % unwrapped phase spectrum

ploc = 1 + find((mX(2:N2-1)>t) .* (mX(2:N2-1)>mX(3:N2)) .* (mX(2:N2-1)>mX(1:N2-2))); % peakss
pmag = mX(ploc); % magnitude of peaks
pphase = pX(ploc); % phase of peaks
posipeaks=length(ploc)/2; %Just positive peaks
    
axis=linspace(0,pi,N2); %Axis for the whole frequency plot
axispeaks=axis(ploc(1:posipeaks));%Axis for just the peaks plot

subplot(2,1,1); plot(axis,mX); hold off; %Magnitude plot
plot(axispeaks,pmag(1:posipeaks),'rx');%Magnitude peaks
subplot(2,1,2); plot(axis,pX); hold off; %Phase plot
plot(axispeaks,pphase(1:posipeaks),'rx');%Phase peaks 

pin = pin+H; % advance sound pointer
end