% Programma per il riconoscimento di note da un racciato audio
% Corso di Abilità Informatiche - Lezione 4-5
%
% 40: DO
% 41: DO#
% 42: RE
% 43: RE#
% 44: MI
% 45: FA
% 46: FA#
% 47: SOL
% 48: SOL#
% 49: LA
% 50: LA#
% 51: SI

fid=fopen('music.inp');
a=textscan(fid,'%f');
a=a{1};

freqsamp=2^13; % Hz
sampling=1/freqsamp; % sec
N=length(a);
t=0:sampling:(N-1)*sampling;

figure(1)
set(gcf,'Name','Grafico del segnale')
plot(t,a)
xlabel('Time')
ylabel('Amplitude')
title('Plot Signal')
axis tight

% Ricerca degli intervalli delle note
aa=a(1:N-2).^2+a(2:N-1).^2+a(3:N).^2;
startnote=find(aa(1:end-1)==0 & aa(2:end)>0);
endnote=find(aa(1:end-1)>0 & aa(2:end)==0);
startnote=[1; startnote];

Notes=length(startnote);

figure(2)
clf;
set(gcf,'Name','Singole note')
for i=1:Notes
    subplot(4,2,i)
    plot(t(startnote(i):endnote(i)),a(startnote(i):endnote(i)))
    xlabel('Time')
    ylabel('Amplitude')
    title(['Nota N. ' num2str(i)])
    axis tight
end

% Intervalli con numero di dati potenza di 2
for i=1:Notes
    nn=endnote(i)-startnote(i)+1;
    fprintf('Nota %i: start: %i; end: %i; n: %i ',...
        i,startnote(i),endnote(i),nn)
    log2nn=floor(log2(nn));
    newnn=2^log2nn;
    notamedia=round((startnote(i)+endnote(i))/2);
    startnote(i)=notamedia-newnn/2;
    endnote(i)=notamedia+newnn/2-1;
    fprintf(' - new n: %i (%i-%i)\n',...
        newnn,startnote(i),endnote(i))
end
  
figure(3)
clf;
set(gcf,'Name','Singole note (potenza di 2)')
for i=1:Notes
    subplot(4,2,i)
    plot(t(startnote(i):endnote(i)),a(startnote(i):endnote(i)))
    xlabel('Time')
    ylabel('Amplitude')
    title(['Nota N. ' num2str(i)])
    axis tight
end

% Stima diretta (distanza tra picchi)
fprintf('\nStima mediante distanza tra picchi\n')
figure(4); clf;
for i=1:Notes
    i1=startnote(i); i2=endnote(i);
    aa=a(i1:i2); tt=t(i1:i2);
    Nn=i2-i1+1;
ind=find(aa(2:Nn-1)>aa(1:Nn-2) & aa(2:Nn-1)>aa(3:Nn));
ind=ind+1; % perche' nel find ind parte da 2
lind=length(ind);
freq=1./(tt(ind(2:lind))-tt(ind(1:lind-1)));
freq1=mean(freq);
err1=std(freq);
note=49+round(12*log(freq1/440)/log(2)); % 49 corrisponde a LA
fprintf('Nota %i - frequenza media: %f +/- %f Hz; Nota %i\n',i,freq1,err1,note)
subplot(Notes,1,i)
plot(tt,aa)
hold on
plot(tt(ind),aa(ind),'*r')
axis tight
end

% Stima mediante autocorrelazione
fprintf('\nAutocorrelazione\n')
figure(5);
clf;
for i=1:Notes
    i1=startnote(i); i2=endnote(i);
    aa=a(i1:i2); tt=t(i1:i2);
    Ncorr=i2-i1+1;
    xcorr=0:sampling:(Ncorr-1)*sampling; % sec
    c=zeros(Ncorr,1);
    for ic=1:Ncorr
        c(ic)=sum(aa(1:Ncorr-ic+1).*aa(ic:Ncorr))/Ncorr;
%        c(ic)=sum(aa(1:Ncorr-ic+1).*aa(ic:Ncorr))/(Ncorr-ic+1);
    end
    subplot(Notes,1,i)
    plot(xcorr,c)
    ind=find(c(2:Ncorr-1)>c(1:Ncorr-2)&c(2:Ncorr-1)>c(3:Ncorr),1);
    freq1=1/xcorr(ind+1);
    note=49+round(12*log(freq1/440)/log(2));
    fprintf('Nota %i: frequenza: %f Hz - note %i\n',i,1/xcorr(ind(1)+1),note);
end

% Stima basata sulla FFT
fprintf('\nStima di Fourier:\n')
figure(6); clf;
set(gcf,'Name','Stima di Fourier')
for i=1:Notes
    i1=startnote(i); i2=endnote(i);
    aa=a(i1:i2); tt=t(i1:i2);
    Nn=i2-i1+1;
    yy=fft(aa);
    pyy=yy.*conj(yy)/Nn;
    pyy=pyy(2:Nn/2); % Tolti primo punto e seconda metà
    fff=freqsamp*(1:(Nn/2-1))/Nn;
    subplot(Notes,1,i)
    plot(fff,pyy)
    if i==Notes, xlabel('Frequency [Hz]'); end
    ylabel('Fourier energy')
    lpyy=length(pyy);
    ind=find(pyy(2:lpyy-1)>pyy(1:lpyy-2)&pyy(2:lpyy-1)>pyy(3:lpyy));
    freq1=fff(ind(1)+1);
    note=49+round(12*log(freq1/440)/log(2));
    fprintf('Nota %i: frequenza: %f Hz - note %i\n',i,fff(ind(1)+1),note);
end

% Autocorrelation FFT
fprintf('\nAutocorrelation (FFT)\n')
figure(7);
clf;
for i=1:Notes
    i1=startnote(i); i2=endnote(i);
    aa=a(i1:i2); tt=t(i1:i2);
    Nn=i2-i1+1;
    xcorr=0:sampling:(Nn-1)*sampling; % sec
    fa=fft(aa);
    c=ifft(fa.*conj(fa))/Nn;
    subplot(Notes,1,i)
    plot(xcorr,c)
    ind=find(c(2:Nn-1)>c(1:Nn-2)&c(2:Nn-1)>c(3:Nn));
    freq1=1/xcorr(ind(1)+1);
    note=49+round(12*log(freq1/440)/log(2));
%    note=49+12*log(freq1/440)/log(2);
    fprintf('Nota %i: frequenza: %f Hz - note %i\n',i,1/xcorr(ind(1)+1),note);
end

sound(a)