yading@10: function [y] = extractMelody(x,w,N,melodyFile)
yading@10: %e.g. y = extractMelody(x,hamming(2025),4096,'baby.txt');
yading@10: 
yading@10: 
yading@10: % Initialize
yading@10: M = length(w);   % analysis window size
yading@10: Ns = 1024;                               % FFT size for synthesis
yading@10: H = 256;                                 % hop size for analysis and synthesis
yading@10: soundlength = length(x);                 % length of input sound array
yading@10: hNs = Ns/2;                              % half synthesis window size
yading@10: hM = (M-1)/2;                            % half analysis window size
yading@10: pin = max(hNs+1,1+hM);   % initialize sound pointer to middle of analysis window
yading@10: pend = soundlength-max(hM,hNs);          % last sample to start a frame
yading@10: fftbuffer = zeros(N,1);                  % initialize buffer for FFT
yading@10: y = zeros(soundlength+Ns/2,1);           % output sine component
yading@10: w = w/sum(w);                            % normalize analysis window
yading@10: sw = zeros(Ns,1);
yading@10: ow = triang(2*H-1);                      % overlapping window
yading@10: ovidx = Ns/2+1-H+1:Ns/2+H;               % overlap indexes
yading@10: sw(ovidx) = ow(1:2*H-1);
yading@10: bh = blackmanharris(Ns);                 % synthesis window
yading@10: bh = bh ./ sum(bh);                      % normalize synthesis window
yading@10: sw(ovidx) = sw(ovidx) ./ bh(ovidx);
yading@10: 
yading@10: 
yading@10: % Load melody file
yading@10: melody = loadMelodyFile(melodyFile);
yading@10: melody = [melody zeros(2,1000)];
yading@10: i=0;
yading@10: 
yading@10: % For each segment
yading@10: while pin<pend
yading@10:     
yading@10:     i=i+1;
yading@10:     
yading@10:     % Find predominant pitch for segment
yading@10:     melody_seg = melody(:,ceil((pin-hM)/128):round((pin+hM)/128));
yading@10:     ind = find(melody_seg(2,:)>150);       % Find pitches larger than 150Hz
yading@10:     medpitch = median(melody_seg(2,ind));  % Median pitch in segment
yading@10:     
yading@10: 
yading@10:     % Compute FFT for segment
yading@10:     xw = x(pin-hM:pin+hM).*w(1:M);         % window the input sound
yading@10:     fftbuffer(1:(M+1)/2) = xw((M+1)/2:M);  % zero-phase window in fftbuffer
yading@10:     fftbuffer(N-(M-1)/2+1:N) = xw(1:(M-1)/2);
yading@10:     X = fft(fftbuffer);                    % compute the FFT
yading@10:     
yading@10:     % Keep only the melody for each segment
yading@10:     if(isempty(ind)) % if there is no melody in segment, then silence segment
yading@10:         X = zeros(N,1);
yading@10:     else             % if there is melody, then mask everything else
yading@10: 
yading@10:         melodyBinsStart = (medpitch/10.6568) * [1:40];
yading@10:         melodyBinsEnd =  N + 2 - ((medpitch/10.6568) * [1:40]);
yading@10:         melodyBins = round([melodyBinsStart melodyBinsEnd]);
yading@10:         melodyBins = [melodyBins (melodyBins-1) (melodyBins+1) (melodyBins-2) (melodyBins+2) (melodyBins-3) (melodyBins+3) (melodyBins-4) (melodyBins+4) (melodyBins-5) (melodyBins+5) (melodyBins-6) (melodyBins+6)];
yading@10:         
yading@10:         C = setdiff(1:N,melodyBins);
yading@10:         
yading@10:         X(C) = 0;
yading@10: 
yading@10:     end;
yading@10:     
yading@10:     
yading@10:     ri= pin-hNs;                           % input sound pointer for residual analysis
yading@10:     yw = ifft(X);    
yading@10:     y(ri:ri+Ns-1) = y(ri:ri+Ns-1)+yw(1:Ns).*sw;
yading@10:     pin = pin+H;     
yading@10:     
yading@10: end
yading@10: 
yading@10: y = (max(x)/max(y))*y; % scale y to original amplitude
yading@10: 
yading@10: %wavwrite(y,44100,'test.wav');