function [coefs] = rand_coefs(f,N)
    % Crea un vector de N elementos con coeficientes a_i aleatorios entre
    % f(i) y -f(i).
    coefs = zeros(N,1);
    for i = 1:N
        coefs(i) = rand*2*f(i)-f(i);
    end
end

function [F] = generar_F(a,b)
    % Genera una función F a partir de sus coeficientes de Fourier.
    
    % a son los N+1 coeficientes para los elementos pares de la base
    % trigonométrica en [-pi,pi], y b son los N coeficientes de los 
    % elementos impares.

    N = length(b);
    n = (1:N)';  % Vector columna
    a = a(:);    % Asegurar que a es columna
    b = b(:);    % Asegurar que b es columna
    
    F = @(x) (a(1)/sqrt(2*pi)) + sum((a(n+1)/sqrt(pi)).*cos(n*x) + (b(n)/sqrt(pi)).*sin(n*x), 1);
end

function [coefs, funciones] = familia_funciones(f,N,Num)
    % f,N son las entradas para rand_coefs
    % Num es el número de funciones de la familia

    % Devuelve:

    % Una lista de parejas de vectores de coeficientes a y b. Para
    % obtener los coeficientes impares de la funcion 3 (por ejemplo),
    % llamas coefs(3).b

    % Una lista de funciones. Para obtener la función 3, llamas
    % funciones{3}

    coefs(Num) = struct('a',[],'b',[]);
    funciones = cell(1,Num);
    for i = 1:Num
        coefs(i).a = rand_coefs(f,N+1);
        coefs(i).b = rand_coefs(f,N);
        funciones{i} = generar_F(coefs(i).a,coefs(i).b); 
    end
end

f1 = @(x) 1./(x.^2); f2 = @(x) 1./x; f3 = @(x) 1/sqrt(x);
N = 500;
Num = 1e5;
[coefs,funciones] = familia_funciones(f2,N,Num); % Elegir la f_i que se quiera
x1 = 1; x2 = 0.2; % Distintos puntos
Sum_1 = 0;
Sum_2 = 0;
for i =1:Num
    Sum_1 = Sum_1 + funciones{i}(x1);
    Sum_2 = Sum_2 + funciones{i}(x2);
end
mean1 = (1/Num)*Sum_1
mean2 = (1/Num)*Sum_2
media_metodo_A = mean1;

S1 = 0;
S2 = 0;
C = 0;
for i =1:Num
    S1 = S1 + (funciones{i}(x1)-mean1)^2;
    S2 = S2 + (funciones{i}(x2)-mean2)^2;
    C = C + (funciones{i}(x1)-mean1)*(funciones{i}(x2)-mean2);
end
VAR1 = (1/(Num-1))*S1
VAR2 = (1/(Num-1))*S2
COV = (1/(Num-1))*C

f1 = @(x) 1./(x.^(2)); f2 = @(x) 1./x; f3 = @(x) 1/sqrt(x);
N = 100;
a1 = rand_coefs(f1,N+1); b1 = rand_coefs(f1,N);
a2 = rand_coefs(f2,N+1); b2 = rand_coefs(f2,N);
a3 = rand_coefs(f3,N+1); b3 = rand_coefs(f3,N);
F1 = generar_F(a1,b1); F2 = generar_F(a2,b2); F3 = generar_F(a3,b3);
x = linspace(-pi,pi,600);
P = linspace(-pi,pi,300);

V=zeros(1000,3);

for i=1:1000

    for k=1:length(P)-1
        V(i,1)=V(i,1)+abs(F1(P(k))-F1(P(k+1)));
    end

    for k=1:length(P)-1
        V(i,2)=V(i,2)+abs(F2(P(k))-F2(P(k+1)));
    end

    for k=1:length(P)-1
        V(i,3)=V(i,3)+abs(F3(P(k))-F3(P(k+1)));
    end
    a1 = rand_coefs(f1,N+1); b1 = rand_coefs(f1,N);
    a2 = rand_coefs(f2,N+1); b2 = rand_coefs(f2,N);
    a3 = rand_coefs(f3,N+1); b3 = rand_coefs(f3,N);
    F1 = generar_F(a1,b1); F2 = generar_F(a2,b2); F3 = generar_F(a3,b3);
end

tablaResultados = table(mean(V)', min(V)', max(V)', std(V)', ...
    'VariableNames', {'Media', 'Minimo', 'Maximo', 'Desviacion Tipica' }, ...
    'RowNames', {'Funcion 1', 'Funcion 2', 'Funcion 3'});

 plot(x,F1(x));
 hold on
 plot(x,F2(x));
 plot(x,F3(x));
 legend