% This is an example script to illustrate the use of the CAdaptInstrSpec
% class. The script loads the spectra extracted two recordings of two
% different violins. It then applied the CAdaptInstrSpec class to adapt the
% spectra of the first violin to those of the second.

clearvars;
close all;

addpath('./functions/');


%% constants

% file paths of instrument spectra
noteSpectraFilePaths = {'./instrSpectra/Violin.CQT.mat'; ...
                        './instrSpectra/Violin2.CQT.mat'};
numInstr = length(noteSpectraFilePaths);

% analysis
tuningFreqInHz = 440;
costFctName = 'LS';
beta = 2; % only needed for plotting. should match the cost function above
numIter = 5;




%% load instrument spectra
instrSpec = cell(numInstr,1);
f0Idcs    = cell(numInstr,1);

for instrIdx = 1:numInstr
    load(noteSpectraFilePaths{instrIdx}); % loads variables: noteSpectra, midiPitches, freqsInHz
    instrSpec{instrIdx} = noteSpectra;
    f0Idcs{instrIdx} = midiPitch2Shift(midiPitches, tuningFreqInHz, freqsInHz)+1;
end
numBinsPerSemitone = getNumBinsPerSemitoneFromFreqVec(freqsInHz);


%% estimate filter curve

% create instrSpecFilter object
instrFiltObj = CAdaptInstrSpec(instrSpec{1}, instrSpec{2}, f0Idcs{1}, f0Idcs{2}, numBinsPerSemitone, costFctName);

% apply update function for h
for iterIdx = 1:numIter;
    instrFiltObj = instrFiltObj.updateH;
end

h = instrFiltObj.getH;
h_smooth = instrFiltObj.getSmoothedH;


%% estimate spectra with given filter curve
estSpectra = instrFiltObj.estimateSpectra(f0Idcs{1});




%% plot results

% plot original spectra and estimated spectra
yticks = [250 500 1000 2000 4000];
yticklabel = {'250'; '500'; '1k'; '2k'; '4k'};

figure;
subplot(311);
imagesc(f0Idcs{1}, freqsInHz, instrSpec{1});
axis xy;
set(gca, 'YScale', 'log', 'YTick', yticks, 'YTickLabel', yticklabel);
xlabel('pitch index');
ylabel('frequency [Hz]');
title('spectra of instrument 1');


subplot(312);
imagesc(f0Idcs{1}, freqsInHz, instrSpec{2});
axis xy;
set(gca, 'YScale', 'log', 'YTick', yticks, 'YTickLabel', yticklabel);
xlabel('pitch index');
ylabel('frequency [Hz]');
title('spectra of instrument 2');

subplot(313);
imagesc(f0Idcs{1}, freqsInHz, estSpectra);
axis xy;
set(gca, 'YScale', 'log', 'YTick', yticks, 'YTickLabel', yticklabel);
xlabel('pitch index');
ylabel('frequency [Hz]');
title('spectra of instrument 2 adapted to instrument 1')


% plot filter curve and errors
xticks = [125 250 500 1000 2000 4000];
xticklabel = {'125'; '250'; '500'; '1k'; '2k'; '4k'};

% compute beta divergence for each element
betaDivsOriginal = betaDivergencePerElement(instrSpec{1}, instrSpec{2}, beta);
betaDivsEstimate = betaDivergencePerElement(instrSpec{1}, estSpectra,   beta);

% scale per-element beta divergence for image plot
numColors = size(colormap,1);
maxDist     = max([betaDivsOriginal(:); betaDivsEstimate(:)]);
betaDivsOriginal = betaDivsOriginal / maxDist * numColors;
betaDivsEstimate = betaDivsEstimate / maxDist * numColors;



figure;
subplot(311)
%stem(freqsInHz(h ~= 0), h(h ~= 0), 'k', 'filled');
plot(freqsInHz(h ~= 0), h(h ~= 0), 'k.');
hold on;
plot(freqsInHz, h_smooth, 'k');
hold off;
axis tight;
set(gca, 'XScale', 'log', 'XTick', xticks, 'XTickLabel', xticklabel);
title('estimated filter curve ''h''');
legend('original', 'smoothed', 'Location', 'NorthWest');

subplot(312);
image(f0Idcs{1}, freqsInHz, betaDivsOriginal);
axis xy;
set(gca, 'YScale', 'log', 'YTick', yticks, 'YTickLabel', yticklabel);
xlabel('pitch index');
ylabel('frequency [Hz]');
title('elementwise differences between original sets of spectra');

subplot(313);
image(f0Idcs{1}, freqsInHz, betaDivsEstimate);
axis xy;
set(gca, 'YScale', 'log', 'YTick', yticks, 'YTickLabel', yticklabel);
xlabel('pitch index');
ylabel('frequency [Hz]');
title('differences between original and adapted set of spectra');