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

Back up before ISMIR

function [y]=sinemodel_(x,w,N,t)
%initializing values
M = length(w); % window size - the longer the more frequency resolution
N2 = N/2+1; % positive part of the spectrum
Ns= 2048; % FFT size for synthesis (even)
H = 512; % analysis/synthesishop size
soundlength = length(x); % length of input sound array - samples
 
fftbuffer = zeros(N,1); % initialize buffer for FFT
 
%Create a loop to step through the sound array x
%initializing the loop
hNs = Ns/2; % half synthesis window size
hM = (M-1)/2; % half analysis window size used to overlap windows
 
pin = max(H+1,1+hM); % initialize sound pointer to middle of analysis window
pend = soundlength-max(H,hM); % last sample to start a frame
 
y = zeros(soundlength,1); % initialize output array
w = w/sum(w); % normalize analysis window
sw = zeros(Ns,1);
ow = triang(2*H-1); % overlapping window
ovidx = Ns/2+1-H+1:Ns/2+H; % overlap indexes
sw(ovidx) = ow(1:2*H-1);
bh = blackmanharris(Ns); % synthesis window
bh = bh ./ sum(bh); % normalize synthesis window
sw(ovidx) = sw(ovidx) ./ bh(ovidx);
 
while pin<pend
    xw = x(pin-hM:pin+hM).*w(1:M)'; % window the input sound - STFT definition
  
    %zero phased window
    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);
 
    %compute FFT of the zero phased frame
    X = fft(fftbuffer);
 
    %calculate  magnitude and phase spectrum of of positive frequencies
    mX = 20*log10(abs(X(1:N2)));
    pX = unwrap(angle(X(1:N2))); % unwrapped phase spectrum
    
        
    %Find the locations, ploc, of the local maxima above a given 
    %threshold, t, in each magnitude spectrum by finding changes of slope.
    ploc = 1+find((mX(2:N2-1)>t).*(mX(2:N2-1)>mX(3:N2)).*(mX(2:N2-1)>mX(1:N2-2))); %peaks
    
    %Find the magnitudes, pmag, and phases, pphase, of the obtained
    %locations.
    pmag = mX(ploc);
    %pmag = mX(ploc)*0.4; 
    pphase = pX(ploc); 
    
    %peak interpolation
    [iploc, ipmag, ipphase] = peakinterp (mX, pX, ploc);
    
    %plot for a window
   
        subplot(2,1,1)
        plot(mX)
        hold on
        plot(ploc,pmag,'*');
        plot(iploc,ipmag,'c');
        title('magnitude peak values');
        hold off
        subplot(2,1,2)
        plot(pX) 
        hold on
        plot(ploc,pphase,'m');
        plot(iploc,ipphase,'c*')
        title('phase peak values');
        hold off
        
        %number of peaks
        
end