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

Back up before ISMIR

function y=stpt2(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 
  
%Internal variables 
M=length(w); 
N2=N/2+1; 
soundlength=length(x); 
hM=(M-1)/2; 
  
%Pointers 
pin=1+hM; 
pend=soundlength-hM; 
  
%Initialize buffers 
fftbuffer=zeros(N,1); 
yw=zeros(M,1); 
y=zeros(soundlength,1); 
  
%Analysis/Synthesis windows 
sw=hanning(M); 
sw=sw./sum(sw); 
w=w./sum(w); 
  
while pin<pend 
    %Analysis 
    xw=x(pin-hM:pin+hM).*w(1:M); %Frame selector 
    fftbuffer=zeros(N,1);        %Zero-phase windowing 
    fftbuffer(1:(M+1)/2)=xw((M+1)/2:M); 
    fftbuffer(N-(M-1)/2+1:N)=xw(1:(M-1)/2); 
     
    X=fft(fftbuffer);   %FFT performance 
    mX=20*log10(abs(X(1:N2))); %Compute the FFT Magnitude (dB) 
    pX=angle(X(1:N2));  %Compute the FFT phase 
     
    ploc=1+find((mX(2:N2-1)>t).*(mX(2:N2-1)>mX(3:N2)).*(mX(2:N2-1)>mX(1:N2-2)));      %Peak finder above threshold (t)    
    pmag=mX(ploc);     %Peak magnitude recovery 
    pphase=pX(ploc);   %Peak phase recovery 
 
    %Synthesis 
    Y=zeros(N,1);    %Memory reserve 
    Y(ploc)=10.^(pmag/2).*exp(1i.*pphase); %Spectrum reconstruction from peak info 
    Y(N+2-ploc)=10.^(pmag/20).*exp(-1i.*pphase); %Spectrum reconstruction negative freq. 
    fftbuffer=real(ifft(Y));     %Compute inverse FFT 
    yw((M+1)/2:M)=fftbuffer(1:(M+1)/2); %Zero-phase windowing undo 
    yw(1:(M-1)/2)=fftbuffer(N-(M-1)/2+1:N); 
    y(pin-hM:pin+hM)=y(pin-hM:pin+hM)+H*N*sw.*yw(1:M); %overlap-add 
    pin=pin+H;                   %Pointer advance 
end