% Analisis de la dependencia de las probabilidades de falsa alarma y de no
% alarma en funcion de la relacion señal a ruido

% Fecha de creacion: 10 de Abril de 2007

N=4960;
num_band=33;
over_lap=31/32;
finf=300;
fsup=2000;
alpha=0.18;                  % umbral para considerar deteccion correcta

time_detect=8000;       % numero de muestras correspondientes al tiempo
                        % del trozo de audio que se coge. 24000 se 
                        % corresponden a 3segundos*8000muestras/segundo


% Calculo de los fingerprints de los anuncios para la base de datos

f_22_23_01_8000_8=fingerprint_actualizable(wavread('Anuncios\Procesado\22_23_01_8000_8'),N,over_lap,num_band,finf,fsup);
f_22_23_02_8000_8=fingerprint_actualizable(wavread('Anuncios\Procesado\22_23_02_8000_8'),N,over_lap,num_band,finf,fsup);
f_22_23_03_8000_8=fingerprint_actualizable(wavread('Anuncios\Procesado\22_23_03_8000_8'),N,over_lap,num_band,finf,fsup);
f_22_23_04_8000_8=fingerprint_actualizable(wavread('Anuncios\Procesado\22_23_04_8000_8'),N,over_lap,num_band,finf,fsup);
f_01_07_01_8000_8=fingerprint_actualizable(wavread('Anuncios\Procesado\01_07_01_8000_8'),N,over_lap,num_band,finf,fsup);
f_08_19_01_8000_8=fingerprint_actualizable(wavread('Anuncios\Procesado\08_19_01_8000_8'),N,over_lap,num_band,finf,fsup);
f_09_19_01_8000_8=fingerprint_actualizable(wavread('Anuncios\Procesado\09_19_01_8000_8'),N,over_lap,num_band,finf,fsup);
f_10_19_01_8000_8=fingerprint_actualizable(wavread('Anuncios\Procesado\10_19_01_8000_8'),N,over_lap,num_band,finf,fsup);
f_11_19_01_8000_8=fingerprint_actualizable(wavread('Anuncios\Procesado\11_19_01_8000_8'),N,over_lap,num_band,finf,fsup);
f_12_19_01_8000_8=fingerprint_actualizable(wavread('Anuncios\Procesado\12_19_01_8000_8'),N,over_lap,num_band,finf,fsup);
f_13_19_01_8000_8=fingerprint_actualizable(wavread('Anuncios\Procesado\13_19_01_8000_8'),N,over_lap,num_band,finf,fsup);
f_14_19_01_8000_8=fingerprint_actualizable(wavread('Anuncios\Procesado\14_19_01_8000_8'),N,over_lap,num_band,finf,fsup);
f_17_19_01_8000_8=fingerprint_actualizable(wavread('Anuncios\Procesado\17_19_01_8000_8'),N,over_lap,num_band,finf,fsup);
f_18_23_01_8000_8=fingerprint_actualizable(wavread('Anuncios\Procesado\18_23_01_8000_8'),N,over_lap,num_band,finf,fsup);
f_20_23_01_8000_8=fingerprint_actualizable(wavread('Anuncios\Procesado\20_23_01_8000_8'),N,over_lap,num_band,finf,fsup);


% Simulamos como si se cogiera un intervalo de tiempo de una emision de
% radio a la cual le sumamos un ruido blanco gaussiano. 
% Luego se analiza si se corresponde con algun anuncio de la base de datos 
% y las probabilidades de error (Pn y Pf)

anuncio_detectado=wavread('Anuncios\Procesado\22_23_01_8000_8');
energia=calc_energia(anuncio_detectado);
long=length(anuncio_detectado);

p=1;

for snr=5:60
    
    energia_ruido=energia*10^(-snr/10);
    amp_noise=sqrt(energia_ruido/long);
    ruido=amp_noise*randn(long,1);
    senal_captada=anuncio_detectado+ruido;

    t=1;
    s=1;
    key=1;
    for repeticion=1:1000                                        
         l=round((length(senal_captada)-time_detect)*rand(1)); 

         finger_detectado=fingerprint_actualizable(senal_captada(l:l+time_detect-1,1),N,over_lap,num_band,finf,fsup);

         if key==1
            key=0;
            [m,n]=size(finger_detectado);
            num_bits=m*n;    %numero de bits que se comparan en la busqueda
            probab_alarm=zeros(1,m*n);
            probab_no_alarm=zeros(1,m*n);
         end

         j=1;
         distancia(j)=min(busqueda(finger_detectado,f_22_23_01_8000_8));
         if distancia(j)~=0
             probab_alarm(1,distancia(j))=probab_alarm(1,distancia(j))+1;
         end
         dist_alarm(p,t)=distancia(j);      
         t=t+1;
         j=j+1;
         distancia(j)=min(busqueda(finger_detectado,f_22_23_02_8000_8));
         if distancia(j)~=0
             probab_alarm(1,distancia(j))=probab_alarm(1,distancia(j))+1;
         end
         dist_alarm(p,t)=distancia(j);
         t=t+1;
         j=j+1;
         distancia(j)=min(busqueda(finger_detectado,f_22_23_03_8000_8));
         if distancia(j)~=0
              probab_alarm(1,distancia(j))=probab_alarm(1,distancia(j))+1;
         end
         dist_alarm(p,t)=distancia(j);
         t=t+1;
         j=j+1;
         distancia(j)=min(busqueda(finger_detectado,f_22_23_04_8000_8));
         if distancia(j)~=0
              probab_alarm(1,distancia(j))=probab_alarm(1,distancia(j))+1;
         end
         dist_alarm(p,t)=distancia(j);
         t=t+1;
         j=j+1;
         distancia(j)=min(busqueda(finger_detectado,f_01_07_01_8000_8));
         probab_no_alarm(1,distancia(j))=probab_no_alarm(1,distancia(j))+1;
         dist_no_alarm(p,s)=distancia(j);     % almacena distancia en deteccion correcta para calculo de probabilidad de falsa alarma
         s=s+1;
         j=j+1;
         distancia(j)=min(busqueda(finger_detectado,f_08_19_01_8000_8));
         probab_no_alarm(1,distancia(j))=probab_no_alarm(1,distancia(j))+1;
         dist_no_alarm(p,s)=distancia(j);
         s=s+1;
         j=j+1;
         distancia(j)=min(busqueda(finger_detectado,f_09_19_01_8000_8));
         probab_no_alarm(1,distancia(j))=probab_no_alarm(1,distancia(j))+1;
         dist_no_alarm(p,s)=distancia(j);
         s=s+1;
         j=j+1;
         distancia(j)=min(busqueda(finger_detectado,f_10_19_01_8000_8));
         probab_no_alarm(1,distancia(j))=probab_no_alarm(1,distancia(j))+1;
         dist_no_alarm(p,s)=distancia(j);
         s=s+1;
         j=j+1;
         distancia(j)=min(busqueda(finger_detectado,f_11_19_01_8000_8));
         probab_no_alarm(1,distancia(j))=probab_no_alarm(1,distancia(j))+1;
         dist_no_alarm(p,s)=distancia(j);
         s=s+1;
         j=j+1;
         distancia(j)=min(busqueda(finger_detectado,f_12_19_01_8000_8));
         probab_no_alarm(1,distancia(j))=probab_no_alarm(1,distancia(j))+1;
         dist_no_alarm(p,s)=distancia(j);
         s=s+1;
         j=j+1;
         distancia(j)=min(busqueda(finger_detectado,f_13_19_01_8000_8));
         probab_no_alarm(1,distancia(j))=probab_no_alarm(1,distancia(j))+1;
         dist_no_alarm(p,s)=distancia(j);
         s=s+1;
         j=j+1;
         distancia(j)=min(busqueda(finger_detectado,f_14_19_01_8000_8));
         probab_no_alarm(1,distancia(j))=probab_no_alarm(1,distancia(j))+1;
         dist_no_alarm(p,s)=distancia(j);
         s=s+1;
         j=j+1;
         distancia(j)=min(busqueda(finger_detectado,f_17_19_01_8000_8));
         probab_no_alarm(1,distancia(j))=probab_no_alarm(1,distancia(j))+1;
         dist_no_alarm(p,s)=distancia(j);
         s=s+1;
         j=j+1;
         distancia(j)=min(busqueda(finger_detectado,f_18_23_01_8000_8));
         probab_no_alarm(1,distancia(j))=probab_no_alarm(1,distancia(j))+1;
         dist_no_alarm(p,s)=distancia(j);
         s=s+1;
         j=j+1;
         distancia(j)=min(busqueda(finger_detectado,f_20_23_01_8000_8));
         probab_no_alarm(1,distancia(j))=probab_no_alarm(1,distancia(j))+1;
         dist_no_alarm(p,s)=distancia(j);
         s=s+1;


    end

    SNR(p)=snr;
    
    % Calculo de la BER entre fingerprints y su desviacion tipica
    ber_media(p)=(mean(dist_alarm(p,:)))/(m*n);
    ber_desv(p)=(std(dist_alarm(p,:)))/(m*n);
    
    % Dibujamos algunas PDFs
    if (snr==5 || snr==10 || snr==15|| snr==20 || snr==25 || snr==30)
        figure;
        plot((1:m*n)/m/n,probab_alarm/sum(probab_alarm),'rx')
        hold on;
        plot((1:m*n)/m/n,probab_no_alarm/sum(probab_no_alarm),'yx')
        plot((1:m*n)/m/n,zeros(1,m*n),'kx')
        %plot([alpha alpha],[0 max(probab_no_alarm/sum(probab_no_alarm))])
        title(['SNR (dB) =',num2str(snr)]);
        xlabel('Tasa de error de bit entre fingerprints');
        ylabel('pdf');
        hold off;
    end
 
    % Para evitar problemas con 'wblfit'
    if(min(dist_alarm(p,:))==0)      
        dist_alarm(p,:)=dist_alarm(p,:)+10^(-6);
    end

    % Calculo de los parametros de la distribucion que mas se aproxima a 
    % cada uno de los dos casos:
    [a,b]=wblfit(dist_alarm(p,:));
    [c,d]=normfit(dist_no_alarm(p,:));


    % Probabilidad de producirse una deteccion incorrecta

    prob1(p)=1-wblcdf(alpha,a(1,1)/(m*n),a(1,2));
    prob2(p)=normcdf(alpha,c/(m*n),d/(m*n));

    p=p+1;
end

% Dibujamos la relacion Pn con SNR
figure;
hold on;
plot(SNR,prob1);

xlabel('SNR (db)')

% Dibujamos la relacion Pf con SNR
plot(SNR,prob2,'r');

legend('Pn','Pf');

% Dibujamos el limite por encima del cual consideramos probabilidad de
% error permitida


hold off

% Dibujamos la relacion BER con SNR

figure;
plot(SNR,ber_media);
hold on
plot(SNR,ber_media+(ber_desv/2),'r');
plot(SNR,ber_media-(ber_desv/2),'r');
hold off