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

Back up before ISMIR

function X = stft_peak(x, w, N, H, threshold)
% Analysis/synthesis of a sound using the short-time fourier transform
% x: input sound,
% w: analysis window (odd number of samples required)
% N: FFT size,
% H: hop size
% y: output sound
lengthx = length(x);        % length of sound signal
w = w/sum(w);               % normalize window
W = length(w);              % window size (odd)
W2 = (W-1)/2;               % half window
N2 = N/2+1;                 % half FFT size 

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