annotate extra/Lab4 stptpeaks.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=stptpeaks(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,
yading@10 5 % w: analysis window (odd size),
yading@10 6 % N: FFT size,
yading@10 7 % H: hop size,
yading@10 8 % t: threshold in negative dB,
yading@10 9 % y: output sound
yading@10 10
yading@10 11 M = length(w); % analysis window size
yading@10 12 N2 = N/2+1; % size of positive spectrum
yading@10 13 soundlength = length(x); % length of input sound array
yading@10 14 hM = (M-1)/2; % half analysis window size
yading@10 15
yading@10 16 pin = 1+hM; % initialize sound pointer at the middle of analysis window
yading@10 17 pend = soundlength-hM; % last sample to start a frame
yading@10 18
yading@10 19 fftbuffer = zeros(N,1); % initialize buffer for FFT
yading@10 20 yw = zeros(M,1); % initialize output sound frame
yading@10 21 y = zeros(soundlength,1); % initialize output array
yading@10 22
yading@10 23 w = w/sum(w); % normalize analysis window
yading@10 24 sw = hanning(M); % synthesis window
yading@10 25 sw = sw./sum(sw);
yading@10 26
yading@10 27 while pin<pend
yading@10 28 xw = x(pin-hM:pin+hM).*w(1:M); % window the input sound
yading@10 29 fftbuffer(:) = 0; % reset buffer
yading@10 30 fftbuffer(1:(M+1)/2) = xw((M+1)/2:M); % zero-phase fftbuffer
yading@10 31 fftbuffer(N-(M-1)/2+1:N) = xw(1:(M-1)/2);
yading@10 32
yading@10 33 X = fft(fftbuffer); % compute the FFT
yading@10 34 mX = 20*log10(abs(X(1:N2))); % magnitude spectrum of positive frequencies
yading@10 35 pX = unwrap(angle(X(1:N2))); % unwrapped phase spectrum
yading@10 36
yading@10 37 ploc = 1 + find((mX(2:N2-1)>t) .* (mX(2:N2-1)>mX(3:N2)) .* (mX(2:N2-1)>mX(1:N2-2))); % peakss
yading@10 38 pmag = mX(ploc); % magnitude of peaks
yading@10 39 pphase = pX(ploc); % phase of peaks
yading@10 40 num_peak=length(ploc)/2; %Just positive peaks
yading@10 41
yading@10 42 axis=linspace(0,pi,N2); %Axis for the whole frequency plot
yading@10 43 axis_peak=axis(ploc(1:num_peak));%Axis for just the peaks plot
yading@10 44
yading@10 45 subplot(2,1,1); plot(axis,mX); hold on; %Magnitude plot
yading@10 46 plot(axis_peak,pmag(1:num_peak),'rx');%Magnitude peaks
yading@10 47 subplot(2,1,2); plot(axis,pX); hold on; %Phase plot
yading@10 48 plot(axis_peak,pphase(1:num_peak),'rx');%Phase peaks
yading@10 49
yading@10 50 pin = pin+H; % advance sound pointer
yading@10 51 end