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

Back up before ISMIR
author Yading Song <yading.song@eecs.qmul.ac.uk>
date Thu, 31 Oct 2013 13:17:06 +0000
parents 6840f77b83aa
children
rev   line source
yading@10 1 function y=stpt2(x,w,N,H,t)
yading@10 2 %Analysis/synthesis of a sound using the peaks
yading@10 3 %of the short-time Fourier transform
yading@10 4 %x: input sound, w: analysis window (odd size), N: FFT size, H: hop
yading@10 5 %size, t: threshold in negative dB, y: output sound
yading@10 6
yading@10 7 %Internal variables
yading@10 8 M=length(w);
yading@10 9 N2=N/2+1;
yading@10 10 soundlength=length(x);
yading@10 11 hM=(M-1)/2;
yading@10 12
yading@10 13 %Pointers
yading@10 14 pin=1+hM;
yading@10 15 pend=soundlength-hM;
yading@10 16
yading@10 17 %Initialize buffers
yading@10 18 fftbuffer=zeros(N,1);
yading@10 19 yw=zeros(M,1);
yading@10 20 y=zeros(soundlength,1);
yading@10 21
yading@10 22 %Analysis/Synthesis windows
yading@10 23 sw=hanning(M);
yading@10 24 sw=sw./sum(sw);
yading@10 25 w=w./sum(w);
yading@10 26
yading@10 27 while pin<pend
yading@10 28 %Analysis
yading@10 29 xw=x(pin-hM:pin+hM).*w(1:M); %Frame selector
yading@10 30 fftbuffer=zeros(N,1); %Zero-phase windowing
yading@10 31 fftbuffer(1:(M+1)/2)=xw((M+1)/2:M);
yading@10 32 fftbuffer(N-(M-1)/2+1:N)=xw(1:(M-1)/2);
yading@10 33
yading@10 34 X=fft(fftbuffer); %FFT performance
yading@10 35 mX=20*log10(abs(X(1:N2))); %Compute the FFT Magnitude (dB)
yading@10 36 pX=angle(X(1:N2)); %Compute the FFT phase
yading@10 37
yading@10 38 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)
yading@10 39 pmag=mX(ploc); %Peak magnitude recovery
yading@10 40 pphase=pX(ploc); %Peak phase recovery
yading@10 41
yading@10 42 %Synthesis
yading@10 43 Y=zeros(N,1); %Memory reserve
yading@10 44 Y(ploc)=10.^(pmag/2).*exp(1i.*pphase); %Spectrum reconstruction from peak info
yading@10 45 Y(N+2-ploc)=10.^(pmag/20).*exp(-1i.*pphase); %Spectrum reconstruction negative freq.
yading@10 46 fftbuffer=real(ifft(Y)); %Compute inverse FFT
yading@10 47 yw((M+1)/2:M)=fftbuffer(1:(M+1)/2); %Zero-phase windowing undo
yading@10 48 yw(1:(M-1)/2)=fftbuffer(N-(M-1)/2+1:N);
yading@10 49 y(pin-hM:pin+hM)=y(pin-hM:pin+hM)+H*N*sw.*yw(1:M); %overlap-add
yading@10 50 pin=pin+H; %Pointer advance
yading@10 51 end
yading@10 52
yading@10 53