y@2: %% Default Values
y@2: %Audio and score information
y@2: if(~exist('midiCell','var'))
y@2:     midiCell = cell(2,1);
y@2:     midiCell{1} = 'Lussier_bassoon.mid';
y@2:     midiCell{2} = 'Lussier_trumpet.mid';
y@2:     mixName = 'Lussier_mix.wav';
y@2:     %mixName = '311CLNOF.WAV';
y@2:     %midiCell{1} = '311CLNOF.mid';
y@2: end
y@2: if(~exist('blockSize','var'))   blockSize = 2048;   end
y@2: if(~exist('overlap','var'))     overlap = 2;        end
y@2: %Maximum frequency of the filter
y@2: if(~exist('bandPassFilterMaxFreq','var'))    bandPassFilterMaxFreq = 10000; end
y@2: %Number of band pass filters in the filter
y@2: if(~exist('nbBandPassFilters','var'))    nbBandPassFilters = 20;end
y@2: %Maximum updates per frame
y@2: if(~exist('maxUpdates','var'))    maxUpdates = 2;end
y@2: %Number of partials in the model
y@2: if(~exist('nbPartials','var'))    nbPartials = 8;end
y@2: if(~exist('plotRefinedScore','var'))    plotRefinedScore = 0;end
y@2: if(~exist('plotPartial','var'))    plotPartial = 0;end
y@2: if(~exist('plotFilter','var'))    plotFilter = 0;end
y@2: %Nummber of instrument in the mix.
y@2: nbInstrument = length(midiCell);
y@2: fftSize = blockSize;
y@2: 
y@2: %Weights of the harmonic series
y@2: partialsWeights = ones(nbPartials,nbInstrument);
y@2: 
y@2: colors = 'rwm';
y@2: %% Initialisation
y@2: hopSize =blockSize/overlap;                   %Hopsize
y@2: HFFTS = fftSize/2+1;                    %Size of the usefull part of the fft (half+1)
y@2: [audio, fs] = wavread(mixName);
y@2: 
y@2: % zeropadding the audio length
y@2: audio(end:end+hopSize-mod(length(audio),hopSize)) = 0;
y@2: audioLength = length(audio);
y@2: 
y@2: scoreLength = audioLength/hopSize - overlap+1; %Number of frames
y@2: wdw = hamming(blockSize);
y@2: 
y@2: tic
y@2: %% Create filters
y@2: filterBasis = createBPFilters(bandPassFilterMaxFreq,nbBandPassFilters,fftSize,fftSize/fs);
y@2: filterCoefficients = ones(nbBandPassFilters,nbInstrument);
y@2: %% Creating the score
y@2: originalScore = createScore(scoreLength,midiCell,fs/hopSize,fftSize/fs);
y@2: originalScore(scoreLength+1:end,:) = [];
y@2: gains = double(originalScore'>0);
y@2: %% Creating spectrogram
y@2: magnitude = zeros(fftSize, scoreLength);
y@2: phase = zeros(fftSize, scoreLength);
y@2: audioVector = zeros(fftSize,1);
y@2: for curFrame = 1:scoreLength
y@2:     audioVector(1:blockSize) = audio(hopSize*(curFrame-1)+1:hopSize*(curFrame-1)+blockSize,1).*wdw;
y@2:     spectrogramVector = fft(audioVector);
y@2:     magnitude(:,curFrame) = abs(spectrogramVector);
y@2:     phase(:,curFrame) = angle(spectrogramVector);
y@2: end
y@2: magnitude(HFFTS+1:end,:) = [];
y@2: phase(HFFTS+1:end,:) = [];
y@2: maxSizeWindows = 400;
y@2: hammingWindows = cell(1,maxSizeWindows);
y@2: for len = 1:maxSizeWindows
y@2:     hammingWindows{len} = hamming(len);
y@2: end
y@2: 
y@2: %% Refine Score and create Excitation Spectrum
y@2: refinedScore = originalScore;   %score to be refiened
y@2: basicSpec = zeros(HFFTS, scoreLength);        %Basic theoretical spectrogram
y@2: excitationSpectrums = cell(nbInstrument,scoreLength);   %cell for each frame
y@2: largerES = cell(nbInstrument,scoreLength);   %cell for each frame
y@2: 
y@2: for ins = 1:nbInstrument                %for each isntrument
y@2:     for curFrame = 1: scoreLength        %for each frame
y@2:         excitationSpectrums{ins,curFrame} = sparse(zeros(HFFTS,nbPartials));
y@2:         largerES{ins,curFrame} = sparse(zeros(HFFTS,nbPartials));
y@2:         if(originalScore(curFrame,ins))        %if there is a note
y@2:             %Refine score and creates Excitation Spectrum
y@2:             [refinedScore(curFrame,ins), width] = refineScore(magnitude(:,curFrame),originalScore(curFrame,ins));
y@2:             excitationSpectrums{ins,curFrame} = sparse(createExcitationSpectrums(refinedScore(curFrame,ins),HFFTS,nbPartials,partialsWeights(:,ins),width,hammingWindows));
y@2:             largerES{ins,curFrame} = sparse(createExcitationSpectrums(refinedScore(curFrame,ins),HFFTS,nbPartials,partialsWeights(:,ins),2*width,hammingWindows));
y@2:         end
y@2:         basicSpec(:,curFrame) = basicSpec(:,curFrame) + (excitationSpectrums{ins,curFrame}*partialsWeights(:,ins))*gains(ins,curFrame); %Spectrogram without update
y@2:     end
y@2: end
y@2: %% Plot spectrogram and add original and refined score
y@2: if (plotRefinedScore)
y@2:     figure(1),imagesc(log(10*magnitude(1:fftSize/16,:)+1)), axis xy
y@2:     plotableScore = originalScore+1;
y@2:     plotableScore(plotableScore==1) = 0;
y@2:     hold on ,plot(plotableScore,'.k'),
y@2:     
y@2:     for harm = 1 : nbPartials
y@2:         for ins = nbInstrument:-1:1
y@2:             plot(refinedScore(:,ins)*harm+1,['.' colors(ins)]), colorbar
y@2:         end
y@2:     end
y@2: end
y@2: %% Online Processing
y@2: %Values that can be changed by the model.
y@2: maskForChangeableValues = basicSpec > 10^-4;
y@2: maskedSpectrogram = magnitude.*maskForChangeableValues;
y@2: bufferSize = 100;
y@2: nextProcess = 1;
y@2: for curFrame = 1:scoreLength
y@2:     if sum(refinedScore(curFrame,:),2)
y@2:         UpdateGainCoefficientsOL2
y@2:         if (curFrame==nextProcess)
y@2:             first = max(curFrame-bufferSize+1,1);
y@2:             currentMag = magnitude(:,first:curFrame);
y@2:             tmpMaskForChangeableValues = basicSpec(:,first:curFrame) > 10^-4;
y@2:             tmpMaskedSpectrogram = currentMag.*tmpMaskForChangeableValues;
y@2:             len = curFrame-first+1;
y@2:             for  nbUpdates = 1:maxUpdates           %update loop
y@2:                 UpdatePartialCoefficientsOL
y@2:                 UpdateFilterCoefficientsOL
y@2:                 if (plotFilter)
y@2:                     %%%Filter
y@2:                     figure(4),plot(outlook(1:blockSize/8,:)), title('Filter')
y@2:                 end
y@2:                 if (plotPartial)
y@2:                     %%%Decay
y@2:                     figure(5),stem(partialsWeights), title('Partial Weights')
y@2:                 end
y@2:             end
y@2:             nextProcess = min(ceil(curFrame*1.2),scoreLength)
y@2:             %bufferSize = round(max(bufferSize,curFrame*0.5));
y@2:             if (plotFilter)
y@2:                 %%%Filter
y@2:                 figure(4),plot(outlook(1:blockSize/8,:)), title('Filter')
y@2:             end
y@2:             if (plotPartial)
y@2:                 %%%Decay
y@2:                 figure(5),stem(partialsWeights), title('Partial Weights')
y@2:                 
y@2:             end
y@2:         end
y@2:     end
y@2: end
y@2: timeV = toc;
y@2: fprintf('Time: %6.3g s, %6.2g times real time \t \n ', timeV, timeV*fs/audioLength);
y@2: 
y@2: 
y@2: % Remodel Spectrogram
y@2: ECA = cell(nbInstrument,1);
y@2: xHat = zeros(HFFTS,scoreLength);
y@2: for ins = 1:nbInstrument
y@2:     ECA{ins} = sparse(zeros(HFFTS,scoreLength));
y@2:     for frame = 1:scoreLength
y@2:         ECA{ins}(:,frame) = filterBasis * filterCoefficients(:,ins).*(largerES{ins,frame}*partialsWeights(:,ins));
y@2:     end
y@2:     xHat = xHat + ECA{ins}.*repmat(gains(ins,:),HFFTS,1);
y@2: end
y@2: %% Final plot
y@2: %%%Spectrogram (comparison with basic and final)
y@2: % % bMax = max(max(log(10*basicSpec+1)));
y@2: % % xMax = max(max(log(10*xHat+1)));
y@2: % % mMax = max(max(log(10*magnitude+1)));
y@2: % % tMax = max([bMax mMax xMax]);
y@2: % % figure(3), subplot(3,1,1), imagesc(log(10*basicSpec(1:hopSize/2,:)+1)), axis xy, linkaxes, caxis([0 tMax]), colorbar, title('Basic Model')
y@2: % % subplot(3,1,2), imagesc(log(10*xHat(1:hopSize/2,:)+1)), axis xy, linkaxes, caxis([0 tMax]), colorbar, title('Optimised Spectrogram')
y@2: % % subplot(3,1,3), imagesc(log(10*magnitude(1:hopSize/2,:)+1)), axis xy, linkaxes, caxis([0 tMax]), colorbar, title('Actual Spectrogram')
y@2: % %
y@2: % %
y@2: % % %%%Filter
y@2: % % figure(4),plot(outlook(1:blockSize/8,:)), title('Filter')
y@2: % %
y@2: % % %%%Decay
y@2: % % figure(5),stem(partialsWeights), title('Partial Weights')
y@2: 
y@2: %% Writing
y@2: %extraInfo = [sprintf('_%d_%d_%d_%d_%d_%d',blockSize,overlap,maxUpdates,nbPartials,nbBandPassFilters,bandPassFilterMaxFreq) specExp];
y@2: finalScale = zeros(audioLength,1);
y@2: freq = zeros(fftSize,1);
y@2: finalAudio = cell(nbInstrument+1,1);
y@2: finalAudioMasked = cell(nbInstrument,1);
y@2: 
y@2: wdw2 = hamming(fftSize);
y@2: nbMaskingTechniques = 4;
y@2: maskingCell = cell(1,nbMaskingTechniques);
y@2: maskingCell{1} = 'raw';
y@2: maskingCell{2} = 'per';
y@2: maskingCell{3} = 'sha';
y@2: maskingCell{4} = '10p';
y@2: finalMask = cell(nbInstrument,nbMaskingTechniques); %4 is the number of different masking techniques
y@2: for ins = 1:nbInstrument
y@2:     finalMask{ins,1} = sparse(zeros(HFFTS,scoreLength));
y@2:     finalMask{ins,4} = sparse(zeros(HFFTS,scoreLength));
y@2:     for curFrame = 1:scoreLength
y@2:         finalMask{ins,1}(:,curFrame) = filterBasis * filterCoefficients(:,ins).*(excitationSpectrums{ins,curFrame}*partialsWeights(:,ins));
y@2:         finalMask{ins,4}(:,curFrame) = filterBasis * filterCoefficients(:,ins).*(largerES{ins,curFrame}*partialsWeights(:,ins));
y@2:     end
y@2:     
y@2:     finalMask{ins,1} = finalMask{ins,1}.*repmat(gains(ins,:),HFFTS,1);
y@2:     finalMask{ins,4} = finalMask{ins,4}.*repmat(gains(ins,:),HFFTS,1);
y@2:     
y@2:     % masking strategy 1
y@2:     finalMask{ins,2} = magnitude .* (finalMask{ins,4} ./ (xHat+eps));
y@2:     
y@2:     % masking strategy 2
y@2:     finalMask{ins,3} = magnitude .* ( (finalMask{ins,4}+eps) ./ (xHat+eps));
y@2:     
y@2:     % masking strategy 3
y@2:     finalMask{ins,4} = magnitude .* ( (finalMask{ins,4}+0.1 * eps) ./ (xHat+eps));
y@2:     
y@2:     finalAudio{ins} = zeros(audioLength,1);
y@2:     finalAudioMasked{ins} = zeros(audioLength,1);
y@2: end
y@2: 
y@2: % sonifying the model
y@2: finalAudio{nbInstrument+1} = zeros(audioLength,1);
y@2: %Overlap and add
y@2: for curFrame = 1:scoreLength
y@2:     freq(1:HFFTS) = xHat(:,curFrame).*exp(1i*phase(:,curFrame));
y@2:     freq(HFFTS+1:fftSize) = freq(HFFTS-1:-1:2);
y@2:     faudio = real(ifft(freq)).*wdw2;
y@2:     finalAudio{nbInstrument+1}((curFrame-1)*hopSize+1:(curFrame-1)*hopSize+fftSize,1) = finalAudio{nbInstrument+1}((curFrame-1)*hopSize+1:(curFrame-1)*hopSize+fftSize,1) + faudio;
y@2:     for ins = 1:nbInstrument
y@2:         freq(1:HFFTS) = finalMask{ins,1}(:,curFrame).*exp(1i*phase(:,curFrame));
y@2:         freq(HFFTS+1:fftSize) = freq(HFFTS-1:-1:2);
y@2:         faudio = real(ifft(freq)).*wdw2;
y@2:         finalAudio{ins}((curFrame-1)*hopSize+1:(curFrame-1)*hopSize+fftSize,1) = finalAudio{ins}((curFrame-1)*hopSize+1:(curFrame-1)*hopSize+fftSize,1) + faudio;
y@2:     end
y@2:     finalScale((curFrame-1)*hopSize+1:(curFrame-1)*hopSize+fftSize,1) = finalScale((curFrame-1)*hopSize+1:(curFrame-1)*hopSize+fftSize,1) + wdw2;
y@2: end
y@2: 
y@2: for ins = 1:nbInstrument
y@2:     finalAudio{ins} = finalAudio{ins}./finalScale;
y@2:     %    wavwrite(finalAudio{ins},fs,16,[midiCell{ins}(1:4) '_' num2str(ins) '_' maskingCell{1} extraInfo '.wav']);
y@2: end
y@2: finalAudio{nbInstrument+1} = finalAudio{nbInstrument+1}./finalScale;
y@2: %wavwrite(finalAudio{nbInstrument+1},fs,16,[mixName(1:4) '_0_' maskingCell{1} extraInfo '.wav']);
y@2: 
y@2: % sonifying the masked version
y@2: %Overlap and add
y@2: for tech = 2:nbMaskingTechniques
y@2:     finalScale = zeros(audioLength,1);
y@2:     for curFrame = 1:scoreLength
y@2:         for ins = 1:nbInstrument
y@2:             freq(1:HFFTS) = finalMask{ins,tech}(:,curFrame).*exp(i*phase(:,curFrame));
y@2:             freq(HFFTS+1:fftSize) = freq(HFFTS-1:-1:2);
y@2:             faudio = real(ifft(freq)).*wdw2;
y@2:             finalAudioMasked{ins}((curFrame-1)*hopSize+1:(curFrame-1)*hopSize+fftSize,1) = finalAudioMasked{ins}((curFrame-1)*hopSize+1:(curFrame-1)*hopSize+fftSize,1) + faudio;
y@2:         end
y@2:         finalScale((curFrame-1)*hopSize+1:(curFrame-1)*hopSize+fftSize,1) = finalScale((curFrame-1)*hopSize+1:(curFrame-1)*hopSize+fftSize,1) + wdw2;
y@2:     end
y@2:     for ins = 1:nbInstrument
y@2:         finalAudioMasked{ins} = finalAudioMasked{ins}./finalScale;
y@2:         %    wavwrite(finalAudioMasked{ins},fs,16,[midiCell{ins}(1:4) '_' num2str(ins) '_' maskingCell{tech} extraInfo '.wav']);
y@2:     end
y@2: end
y@2: 
y@2: 
y@2: