function c = analizza_suono(nomefile, freq_base)
%ANALIZZASUONO Mostra le frequenze fondamentali in un file audio. 
%

figure;

if ischar(nomefile)
    [s, Fs] = audioread(nomefile);
else
    s = nomefile;
    Fs = 44100;
end

% Se viene specificata una frequenza di base la usiamo per stimare i
% coefficienti di fourier. 
if ~exist('freq_base', 'var')
    freq_base = 0;
end

% Facciamo una copia del suono per intero che ci servirà più tardi, mentre
% s verrà tagliato a solo una parte iniziale per meglio identificare le
% frequenze. 
s_intero = s;

% Teniamo solo la parte iniziale del suono -- questo ci permette di
% ignorare il fatto che durante un suono lungo i coefficienti di Fourier
% cambiano un po'. 
offset = round(0.25 * 44100); % 1/4 di secondo
s = s(offset : offset + 8192);

% Calcoliamo la trasformata di Fourier discreta dei sample. Questa è una 
% trasformata fatta utilizzando i numeri complessi, e restituisce dei 
% valori da cui si possono ricavare i coefficienti a_i, b_i della serie
% di Fourier come li abbiamo introdotti noi.
sf = fft(s);
sf = sf / max(abs(sf));

% Leggiamo le frequenze prendendo la prima metà del suono (ovvero estraiamo
% gli a_i). 
l = ceil((length(sf) + 1) / 2);
freqs = abs(sf(1 : l));

x = ( (0:l-1) ./ l * Fs / 2);

subplot(1, 2, 1);

ind = floor(Fs * 0.1) : floor(Fs * 0.12);
plot(ind ./ Fs, s_intero(ind));
xlabel('t');
ylabel('s(t)');
title('Forma d''onda');

set(gcf, 'position', [ 400 400 1400 700 ]);

subplot(1, 2, 2);
semilogy(x, freqs);

% Zoomiamo nel range interessante, dove ci sono le frequenze
% ragionevolmente udibili per noi. 
axis([ 0 8000 0 1.1 * max(freqs) ]);

title('Coefficienti di Fourier');
xlabel('Frequenza (Hz)');
ylabel('Coefficiente di Fourier');

% Identificazione della frequenza base
if freq_base ~= 0
    % Ci limitiamo a 8000 Hz, che bastano per ricostruire fedelmente la
    % maggior parte dei suoni
    ntones = floor(8000 / freq_base);
    
    c = zeros(1, ntones);
    
    for j = 1 : ntones
        jj = freq_base * l * 2.0 * j / Fs;
        j1 = round(jj * 0.98);
        j2 = round(jj * 1.02);
        c(j) = max(freqs(j1 : j2));
        
        % Azzeriamo le frequenze sotto un'opportuna scelta di tolleranza,
        % plausibilmente corrispondono a valori nulli, e quello che
        % rileviamo è solo rumore di fondo.
        if c(j) < 1e-2
            c(j) = 0.0;
        end
    end
    
    hold on;
    plot(freq_base * (1 : ntones), c, 'r*');
else
    c = [];
end