%% Ejercicios de estimación paramétrica
% 2022 G.T.A.S. (Grupo Tratamiento Avanzado de Señal). Universidad de Cantabria.
%
%
%% Ejercicio 1:
%
% Deteccion de una señal conocida en WGN
% Modelo:
% H_0: \ x[n] = r[n], \quad n=1,2,...,N-1 \\
% H_1: \ x[n] = s[n] + r[n], \quad n=1,2,...,N-1
% donde s[n] es una señal conocida y r[n] es WGN con varianza \sigma^2.
%
% Datos:

N=10; % numero de muestras de la se?al
varr=1.0; % varianza del ruido

n=(0:N-1)';
w0=pi/4;
s=cos(w0*n); % se?al s(n)
Es=sum(s.^2); % energia de la se?al

P_H0 = 1/2; % probabilidad a priori de H0
P_H1 = 1-P_H0; % probabilidad a priori de H1

% Genera realizaciones del estadistico bajo ambas hipotesis:

Nsim=3e4; % numero de realizaciones del estadistico bajo ambas hipotesis

% Ciclo a las realizaciones bajo H0
T_H0=zeros(Nsim,1);
for ns=1:Nsim 
    r = sqrt(varr) * randn(N,1);  % realizacion de ruido
    x = r;  % realizacion bajo H0
    T_H0(ns) = x'*s; % valor del estadistico para la realizacion
end

% Ciclo a las realizaciones bajo H1
T_H1=zeros(Nsim,1);
for ns=1:Nsim   
    r = sqrt(varr) * randn(N,1);  % realizacion de ruido
    x = s + r; % realizacion bajo H1
    T_H1(ns) = x'*s; % valor del estadistico para la realizacion
end

% Dibuja una realizacion bajo cada hipotesis:

r = sqrt(varr) * randn(N,1);
x_H0 = r;  % realizacion bajo H0
r = sqrt(varr) * randn(N,1);
x_H1 = r + s ; % realizacion bajo H1
figure(10);plot(n,x_H0,'bo',n,x_H1,'ro',n,s,'ko-') % dibuja las realizaciones
figure(10);legend({'x[n] | H0', 'x[n] | H1',' s[n]'},'FontSize',20)
figure(10);title('Ejemplo de realizaciones','FontSize',22)
figure(10);xlabel('n','FontSize',20)
figure(10);axis tight

% Histogramas del estadistico bajo ambas hipotesis:

figure(20);hold on
figure(20); histogram(T_H0,'Normalization','pdf') % histograma de T_H0
figure(20); histogram(T_H1,'Normalization','pdf') % histograma de T_H0
T = [T_H0;T_H1]; 
t = linspace(min(T),max(T),1000); % valores de T para dibujar las pdf analiticas
pdf_T_H0=normpdf(t,0,sqrt(varr*Es)); % pdf analitica de T_H0
pdf_T_H1=normpdf(t,Es,sqrt(varr*Es)); % pdf analitica de T_H1
figure(20);plot(t,pdf_T_H0,'k',t,pdf_T_H1,'k') % dibuja las pdf analiticas
figure(20);title('p.d.f y histogramas del estadistico T(x)','FontSize',22)
figure(20);xlabel('T','FontSize',20)
figure(20);hold off

% Dibuja la RoC y las prestaciones para cada criterio de detección:

% ROC's
labels = [zeros(1,length(T_H0)) ones(1,length(T_H1))];
T=[T_H0' T_H1'];
[Pfa_e,Pd_e,g] = perfcurve(labels,T,'1');

% Pd y Pfa analiticas
Pfa_a=qfunc(g/sqrt(varr*Es)); % Pfa analitica
Pd_a=qfunc((g-Es)/sqrt(varr*Es)); % Pd analitica

% Clasificador de minima Pe
umbral_0 = P_H0/P_H1; % valor optimo del umbral 
umbral = varr*log(umbral_0)+(Es/2);
Pfa=sum(T_H0 > umbral)/Nsim;
Pd=sum(T_H1 > umbral)/Nsim;
Pe = Pfa * P_H0 + (1-Pd) * P_H1;
disp(' ')
disp('Clasificador de minima Pe:')
disp(['     ','      umbral optimo =',num2str(umbral),'       Pe =', num2str(Pe),'       Pfa =',num2str(Pfa),'       Pd =',num2str(Pd) ])
Pfa_minPe = Pfa; Pd_minPe = Pd;

% Clasificador de minimo coste Bayesiano
C_01 = 10; C_10 = 1; C_00 = 0; C_11=0; % costes Bayesianos 
umbral_0 = (P_H0/P_H1)*(C_10-C_00)/(C_01-C_11); % valor optimo del umbral 
umbral = varr*log(umbral_0)+(Es/2);
Pfa=sum(T_H0 > umbral)/Nsim;
Pd=sum(T_H1 > umbral)/Nsim;
Pe = Pfa * P_H0 + (1-Pd) * P_H1;
disp(' ')
disp('Clasificador de minimo coste Bayesiano:')
disp(['     ','      umbral optimo =',num2str(umbral),'       Pe =', num2str(Pe),'       Pfa =',num2str(Pfa),'       Pd =',num2str(Pd) ])
Pfa_minC = Pfa; Pd_minC = Pd;

% Clasificador de Neyman-Pearson
Pfa = 0.1; % probabilidad de falsa alarma deseada
indice_g=find(Pfa_e < Pfa,1,'last');
umbral = g(indice_g);
Pd=sum(T_H1 > umbral)/Nsim;
Pe = Pfa * P_H0 + (1-Pd) * P_H1;
disp(' ')
disp('Clasificador de Neyman-Pearson:')
disp(['     ','      umbral optimo =',num2str(umbral),'       Pe =', num2str(Pe),'       Pfa =',num2str(Pfa),'       Pd =',num2str(Pd) ])
Pfa_NP = Pfa; Pd_NP = Pd;

% Dibuja las ROC y los puntos de la RoC para cada criterio
figure(30); hold on
figure(30);plot(Pfa_e,Pd_e,'b-',Pfa_a,Pd_a,'r-','LineWidth',1.5)
figure(30);xlabel('Pfa','FontSize',20)
figure(30);ylabel('Pd','FontSize',20)
figure(30);title('RoC','FontSize',24)
figure(30);grid on
figure(30);plot(Pfa_minPe,Pd_minPe,'go',Pfa_minC,Pd_minC,'gs',Pfa_NP,Pd_NP,'gd','MarkerSize',14,'LineWidth',2)
figure(30);legend({'RoC experimental','RoC', 'criterio min Pe', 'criterio min C','criterio NP'},'Location','SouthEast','FontSize',20)
figure(30); hold off

%
%% Ejercicio 3:
%
% Modelo:
% H_0: \ x[n] = s_0[n] + r[n], \quad n=1,2,...,N-1 \\
% H_1: \ x[n] = s_1[n] + r[n], \quad n=1,2,...,N-1
% donde s[n] es una señal conocida y r[n] es WGN con varianza \sigma^2.
%
% Datos:

N=6; % numero de muestras de la se?al
varr=1.0; % varianza del ruido

n=(0:N-1)';
s0=cos(pi/2*n); % señal s0(n)
s1=cos(pi/4*n); % señal s1(n)
E0=sum(s0.^2); % energia de la se?al s1[n]
E1=sum(s1.^2); % energia de la se?al s2[n]
E01=sum(s0.*s1);

P_H0 = 1/2; % probabilidad a priori de H0
P_H1 = 1-P_H0; % probabilidad a priori de H1

% Genera realizaciones del estadistico bajo ambas hipotesis:


Nsim=1e6; % numero de realizaciones del estadistico bajo ambas hipotesis

% Ciclo a las realizaciones bajo H0
T_H0=zeros(Nsim,1);
for ns=1:Nsim 
    r = sqrt(varr) * randn(N,1);  % realizacion de ruido
    x = s0 + r;  % realizacion bajo H0
    T_H0(ns) = x'*(s1-s0); % valor del estadistico para la realizacion
end

% Ciclo a las realizaciones bajo H1
T_H1=zeros(Nsim,1);
for ns=1:Nsim   
    r = sqrt(varr) * randn(N,1);  % realizacion de ruido
    x = s1 + r; % realizacion bajo H1
    T_H1(ns) = x'*(s1-s0); % valor del estadistico para la realizacion
end

% Dibuja una realizacion bajo cada hipotesis:

r = sqrt(varr) * randn(N,1);
x_H0 = s0 + r;  % realizacion bajo H0
r = sqrt(varr) * randn(N,1);
x_H1 = s1 + r ; % realizacion bajo H1
figure(10);plot(n,x_H0,'bo',n,x_H1,'ro',n,s0,'b-',n,s1,'r-') % dibuja las realizaciones
figure(10);legend({'x[n] | H0', 'x[n] | H1',' s0[n]',' s1[n]'},'FontSize',20)
figure(10);title('Ejemplo de realizaciones','FontSize',22)
figure(10);xlabel('n','FontSize',20)

% Histogramas del estadistico bajo ambas hipotesis:

figure(20);hold on
figure(20); histogram(T_H0,'Normalization','pdf') % histograma de T_H0
figure(20); histogram(T_H1,'Normalization','pdf') % histograma de T_H0
T = [T_H0;T_H1]; 
t = linspace(min(T),max(T),1000); % valores de T para dibujar las pdf analiticas
pdf_T_H0=normpdf(t,E01-E0,sqrt(varr*(E0+E1-2*E01))); % pdf analitica de T_H0
pdf_T_H1=normpdf(t,E1-E01,sqrt(varr*(E0+E1-2*E01))); % pdf analitica de T_H1
figure(20);plot(t,pdf_T_H0,'k',t,pdf_T_H1,'k','LineWidth',1) % dibuja las pdf analiticas
figure(20);title('p.d.f y histogramas del estadistico T(x)','FontSize',22)
figure(20);xlabel('T','FontSize',20)
figure(20);hold off

% Dibuja la RoC y las prestaciones para cada criterio de detección:

% ROC experimental
labels = [zeros(1,length(T_H0)) ones(1,length(T_H1))];
T=[T_H0' T_H1'];
[Pfa_e,Pd_e,g] = perfcurve(labels,T,'1');

% Pd y Pfa analiticas
Pfa_a=qfunc((g+E0-E01)/sqrt(varr*(E0+E1-2*E01))); % Pfa analitica
Pd_a=qfunc((g+E01-E1)/sqrt(varr*(E0+E1-2*E01))); % Pd analitica

% Clasificador de minima Pe
umbral_0 = P_H0/P_H1; % valor optimo del umbral 
umbral = varr*log(umbral_0)+((E1-E0)/2);
Pfa=sum(T_H0 > umbral)/Nsim;
Pd=sum(T_H1 > umbral)/Nsim;
Pe = Pfa * P_H0 + (1-Pd) * P_H1;
disp(' ')
disp('Clasificador de minima Pe:')
disp(['     ','      umbral optimo =',num2str(umbral),'       Pe =', num2str(Pe),'       Pfa =',num2str(Pfa),'       Pd =',num2str(Pd) ])
Pfa_minPe = Pfa; Pd_minPe = Pd;

% Clasificador de minimo coste Bayesiano
C_01 = 10; C_10 = 1; C_00 = 0; C_11=0; % costes Bayesianos 
umbral_0 = (P_H0/P_H1)*(C_10-C_00)/(C_01-C_11); % valor optimo del umbral 
umbral = varr*log(umbral_0)+((E1-E0)/2);
Pfa=sum(T_H0 > umbral)/Nsim;
Pd=sum(T_H1 > umbral)/Nsim;
Pe = Pfa * P_H0 + (1-Pd) * P_H1;
disp(' ')
disp('Clasificador de minimo coste Bayesiano:')
disp(['     ','      umbral optimo =',num2str(umbral),'       Pe =', num2str(Pe),'       Pfa =',num2str(Pfa),'       Pd =',num2str(Pd) ])
Pfa_minC = Pfa; Pd_minC = Pd;

% Clasificador de Neyman-Pearson
Pfa = 0.05; % probabilidad de falsa alarma deseada
indice_g=find(Pfa_e < Pfa,1,'last');
umbral = g(indice_g);
Pd=sum(T_H1 > umbral)/Nsim;
Pe = Pfa * P_H0 + (1-Pd) * P_H1;
disp(' ')
disp('Clasificador de Neyman-Pearson:')
disp(['     ','      umbral optimo =',num2str(umbral),'       Pe =', num2str(Pe),'       Pfa =',num2str(Pfa),'       Pd =',num2str(Pd) ])
Pfa_NP = Pfa; Pd_NP = Pd;

% Dibuja las ROC y los puntos de la RoC para cada criterio
figure(30); hold on
figure(30);plot(Pfa_e,Pd_e,'b-',Pfa_a,Pd_a,'r-','LineWidth',1)
figure(30);xlabel('Pfa','FontSize',20)
figure(30);ylabel('Pd','FontSize',20)
figure(30);title('RoC','FontSize',24)
figure(30);grid on
figure(30);plot(Pfa_minPe,Pd_minPe,'go',Pfa_minC,Pd_minC,'gs',Pfa_NP,Pd_NP,'gd','MarkerSize',14,'LineWidth',2)
figure(30);legend({'RoC experimental','RoC', 'criterio min Pe', 'criterio min C','criterio NP'},'Location','SouthEast','FontSize',20)
figure(30); hold off

%
%% Ejercicio 5:
%
% Detector de energía
%
% Datos:

N=10; % numero de muestras de la se?al
n=(0:N-1)';

var0=0.5; % varianza del ruido para H0
var1=1.0; % varianza del ruido para H1. Suponemos que var1 > var0

P_H0 = 3/4; % probabilidad a priori de H0
P_H1 = 1-P_H0; % probabilidad a priori de H1

% Genera realizaciones del estadistico bajo ambas hipotesis:

Nsim=3e4; % numero de realizaciones del estadistico bajo ambas hipotesis

% Ciclo a las realizaciones bajo H0
T_H0=zeros(Nsim,1);
for ns=1:Nsim 
    x = sqrt(var0) * randn(N,1);  % realizacion bajo H0
    T_H0(ns) = x'*x; % valor del estadistico para la realizacion
end

% Ciclo a las realizaciones bajo H1
T_H1=zeros(Nsim,1);
for ns=1:Nsim   
    x =  sqrt(var1) * randn(N,1); % realizacion bajo H1
    T_H1(ns) = x'*x; % valor del estadistico para la realizacion
end

% Dibuja una realizacion bajo cada hipotesis:

x_H0 = sqrt(var0) * randn(N,1);  % realizacion bajo H0
x_H1 = sqrt(var1) * randn(N,1); % realizacion bajo H1
figure(10);plot(n,x_H0,'bo-',n,x_H1,'ro-') % dibuja las realizaciones
figure(10);legend({'x[n] | H0', 'x[n] | H1'},'FontSize',20)
figure(10);title('Ejemplo de realizaciones','FontSize',22)
figure(10);xlabel('n','FontSize',20)

% Histogramas del estadistico bajo ambas hipotesis:

figure(20);hold on
figure(20); histogram(T_H0,'Normalization','pdf') % histograma de T_H0
figure(20); histogram(T_H1,'Normalization','pdf') % histograma de T_H0
figure(20);title('Histogramas de los estadisticos T(x|H0) y T(x|H1)','FontSize',20)
figure(20);xlabel('T','FontSize',20)
figure(20);hold off

% Dibuja la RoC y las prestaciones para cada criterio de clasificación:

% ROC experimental
labels = [zeros(1,Nsim) ones(1,Nsim)];
T=[T_H0' T_H1'];
[Pfa_e,Pd_e,umbrales] = perfcurve(labels,T,'1');

% Clasificador de minima Pe
umbral_0 = P_H0/P_H1; 
umbral = (2*log(umbral_0)+N*log(var1/var0))/(1/var0-1/var1); 
Pfa=sum(T_H0 > umbral)/Nsim;
Pd=sum(T_H1 > umbral)/Nsim;
Pe = Pfa * P_H0 + (1-Pd) * P_H1;
disp(' ')
disp('Clasificador de minima Pe:')
disp(['     ','      umbral optimo =',num2str(umbral),'       Pe =', num2str(Pe),'       Pfa =',num2str(Pfa),'       Pd =',num2str(Pd) ])
Pfa_minPe = Pfa; Pd_minPe = Pd;

% Clasificador de minimo coste Bayesiano
C_01 = 10; C_10 = 1; C_00 = 0; C_11=0; % costes Bayesianos 
umbral_0 = (P_H0/P_H1)*(C_10-C_00)/(C_01-C_11); 
umbral = (2*log(umbral_0)+N*log(var1/var0))/(1/var0-1/var1);
Pfa=sum(T_H0 > umbral)/Nsim;
Pd=sum(T_H1 > umbral)/Nsim;
Pe = Pfa * P_H0 + (1-Pd) * P_H1;
disp(' ')
disp('Clasificador de minimo coste Bayesiano:')
disp(['     ','      umbral optimo =',num2str(umbral),'       Pe =', num2str(Pe),'       Pfa =',num2str(Pfa),'       Pd =',num2str(Pd) ])
Pfa_minC = Pfa; Pd_minC = Pd;

% Clasificador de Neyman-Pearson
Pfa = 0.1; % probabilidad de falsa alarma deseada
indice=find(Pfa_e > Pfa,1,'first');
umbral = umbrales(indice);
Pd = Pd_e(indice);
Pe = Pfa * P_H0 + (1-Pd) * P_H1;
disp(' ')
disp('Clasificador de Neyman-Pearson:')
disp(['     ','      umbral optimo =',num2str(umbral),'       Pe =', num2str(Pe),'       Pfa =',num2str(Pfa),'       Pd =',num2str(Pd) ])
Pfa_NP = Pfa; Pd_NP = Pd;

% Dibuja las ROC y los puntos de la RoC para cada criterio
figure(30); hold on
figure(30);plot(Pfa_e,Pd_e,'b-','LineWidth',1.5)
figure(30);xlabel('$P_{FA}(\gamma)$','FontSize',20,'Interpreter','latex')
figure(30);ylabel('$P_{D}(\gamma)$','FontSize',20,'Interpreter','latex')
figure(30);title('RoC experimental','FontSize',20)
figure(30);grid on
figure(30);plot(Pfa_minPe,Pd_minPe,'go',Pfa_minC,Pd_minC,'gs',Pfa_NP,Pd_NP,'gd','MarkerSize',14,'LineWidth',2)
figure(30);legend({'RoC experimental', 'criterio min Pe', 'criterio min C','criterio NP'},'Location','SouthEast','FontSize',20)
figure(30); hold off

%
%%
%% Ejercicio 8:
%
% Modelo:  H_m : x[n]=s_m[n]+r[n], \quad n=0,...,N-1, \quad m=0,...,M-1
%
% Datos:

M=3;
N=5;

s0 = [0 1 -1 1 0]';
s1 = [1 -1 2 -1 0]';
s2 = [-1 1 -2 0 1]';

n=(0:N-1)';

varr=0.5; % varianza del ruido para H0

P_H0 = 2/M; % probabilidad a priori de H0
P_H1 = (1/2)/M; % probabilidad a priori de H1
P_H2 = (1/2)/M; % probabilidad a priori de H1

% Genera realizaciones del estadistico bajo las hipotesis:

Nr=1e5; % numero de realizaciones del estadistico bajo ambas hipotesis

clase = randsrc(Nr,1,[0 1 2; P_H0 P_H1 P_H2]);

g=zeros(M,1); % inicializa las funciones discriminantes
clase_e = zeros(Nr,1); % clase estimada

% Ciclo a las realizaciones
for nr=1:Nr 
    r = sqrt(varr) * randn(N,1);
    if clase(nr) == 0
        x = s0 + r;
    end
    if clase(nr) == 1
        x = s1 + r;
    end
    if clase(nr) == 2
        x = s2 + r;
    end
    
    g(1) = 2*varr*log(P_H0) - norm(x-s0)^2; 
    g(2) = 2*varr*log(P_H1) - norm(x-s1)^2;
    g(3) = 2*varr*log(P_H2) - norm(x-s2)^2;
    
    [maximo,clase_e(nr)] = max(g);
    clase_e(nr) = clase_e(nr)-1;
end

% Probabilidad de error:

Pe = 1-sum(clase_e==clase)/Nr;
disp(' ')
disp(['Pe: ',num2str(Pe)])

% Matriz de confusion:

MC = confusionmat(clase,clase_e);
disp(' ')
disp('Matriz de confusion:')
disp(MC)

%
%% Ejercicio 10:
%
% Datos:

M=4;
s=[-3 -1 1 3]'; % simbolos modulacion M-PAM
varr=0.25; % varianza del ruido

% Realizaciones bajo las hipótesis:

Nr=1e6; % Numero de realizaciones

% Realizaciones de hipótesis (equiprobables). Se nombran desde H_1 hasta H_M
clase=randi(4,Nr,1);

%%%% Ruido blanco Gaussiano

clase_e = zeros(Nr,1); % clase estimada

% Realizaciones de ruido blanco Gaussiano
r = sqrt(varr) * randn(Nr,1);

% Ciclo a las realizaciones
for nr=1:Nr 
    x = s(clase(nr)) + r(nr);
    g = abs(x-s);    
    [minimo,clase_e(nr)] = min(g);
end

% Probabilidad de error:

Pe = 1-sum(clase_e==clase)/Nr;
disp(' ')
disp(['Pe ruido blanco Gaussiano: ',num2str(Pe)])

% Ruido Laplaciano:

clase_e = zeros(Nr,1); % clase estimada

% Realizaciones independientes de ruido Laplaciano
a = sqrt(2/varr);
r = ruido_laplaciano(Nr,a);

% Ciclo a las realizaciones
for nr=1:Nr 
    x = s(clase(nr)) + r(nr);
    g = abs(x-s);  
    [minimo,clase_e(nr)] = min(g);
end

% Probabilidad de error
Pe = 1-sum(clase_e==clase)/Nr;
disp(['Pe ruido i.i.d. Laplaciano: ',num2str(Pe)])

%

function r=ruido_laplaciano(Nr,varr)

% Genera realizaciones independientes de ruido Laplaciano
%
% Variables de entrada
% - Nr: numeronde realizaciones de ruido
% - varr: varianza de ruido

u=rand(Nr,1);
u=u-1/2;
a=sqrt(2/varr);
r=(-1/a)*sign(u).*log(1-2*abs(u));
end