function [harmonic_ratios_vector, harmonic_frequency_vector, fs] = FindHarmonics(note_fundamental_frequency, filename)

% Load the input note wave file
[note_signal, fs, N] = wavread(filename);

% Calculate the total number of samples
npoints = length(note_signal);

% Create the appropriate frequency resolution vector up to fs/2
f = fs * [0:1:npoints/2-1]' / npoints;

% Calculate the input note spectrum (linear magnitude)
Plinear = abs(fft(note_signal));
% Keep only the first symmetric half
Plinear = Plinear(1:npoints/2);

% Get the frequency coefficient that corresponds to the note fundamental frequency
i = find(Plinear == max(Plinear));
% Calcte and display the measured fundamental frequency
found_note_fundamental_frequency = i * (f(2) - f(1));
disp('Harmonic ratios calculator monitor');
disp(sprintf('%s%3.2f','Original note fundamental frequency = ',note_fundamental_frequency));
disp(sprintf('%s%3.2f','Measured note fundamental frequency = ',found_note_fundamental_frequency));

% Now get the even and odd harmonics' ratios starting from the first even
harmonic_ratios_vector = [1];           % Initiated to 1 for the note fundamental frequency
harmonic_frequency_vector = [found_note_fundamental_frequency];
harmonic = 2;

while (harmonic * f(i) <= fs/2),
    current_ratio = Plinear(harmonic * i) / Plinear(i);
    harmonic_ratios_vector = [harmonic_ratios_vector; current_ratio];
    harmonic_frequency_vector = [harmonic_frequency_vector; harmonic * i * (f(2) - f(1))];
    harmonic = harmonic + 1;
end