function [C,freq] = coherencia(ARcoeff,E,f,fm,p) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Cálculo de la coherencia por métodos paramétricos:
%
%  ARcoeff : Coeficientes del modelo AR (estimados por arfit)
%  E : Matriz del error de predicción (dada por arfit)
%  f : Vector de frecuencias
%  fm : Frecuencia de muestreo
%  p : Orden máximo del modelo AR
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

m = size(ARcoeff,1);       % Número de vías
A = [eye(min(size(ARcoeff))) -ARcoeff];     % Tomamos todos los coefs AR, incluida la identidad
N = length(f);

if f(1) == -fm/2,       %%% Conservamos solo la parte positiva del vector de frecuencias
    freq = 0:fm/N:fm/2-fm/N;
else
    freq = f(f>=0);
end

for i1 = 1:length(freq)
    f1 = freq(i1);
    %%% Cálculo de la matriz de transferencia del sistema 
    H = zeros(m,m);
    for i2 = 0:p
        H = H+A(:,i2*m+1:(i2+1)*m)*exp(-i2*i*2*pi*f1/fm); 
    end
    H = H^-1;    
    %%% Cálculo de la matriz de densidad espectral S(f)    
    S(:,(i1-1)*m+1:(i1-1)*m+m) = H * E * H';
end;

C = zeros(m,m*length(freq));

%%% Cálculo de la coherencia
for i3 = 1:length(freq)
    for i4 = 1:m
        for i5 = 1:m
           C(i4,(i3-1)*m+i5) = abs(S(i4,(i3-1)*m+i5))./sqrt((abs(S(i4,(i3-1)*m+i4)) * abs(S(i5,(i3-1)*m+i5))));
        end
    end
end