function [d] = Calculation (f1, f2, g, N, T, df1, df2, P)
  alpha1 = ((pi * df1).^2) / (4*log(2));
  alpha2 = ((pi * df2).^2) / (4*log(2));


  for(i = 1:N)
    for(j = 1:N)
      K(i, j) = eps.^(-1*alpha1*((i-j)*T).^2);
      Z(i, j) = (eps.^(j*2*pi*f1*T*i)) / sqrt(N);
      S(i, j) = (eps.^(j*2*pi*f2*T*i)) / sqrt(N);
      R(i, j) = P*eps.^(-1*alpha2*((i-j)*T).^2);
    endfor
  endfor

  [Ps, mu] = eig(K);
  mu = diag(mu);

  RA = R';
  SA = S';
  A = Z * Ps';
  BA = SA*Ps;
  B = BA';

  for(l = 1:N)
    for(k = 1:N)
      C(l, k) = 0;
      for(i = 1:N)
        C(l, k) += RA(l, i)*Ps(i, k);
      endfor
    endfor
  endfor
  D = C.'*Ps';


  if(g == false)
    window2 = figure('Name', 'Graphic');
    grid on;
    hold on;
    [mu_sorted, origpos] = sort(mu, 'descend');
    plot(1:N, A(origpos(1:4), :));
    plot(1:N, K(N/2, :));
    legend('1', '2', '3', '4', 'N/2 row of K');
  endif


  d = 0; den = 0;
  for(k = 1:N)
    d += (A(f1, k)*BA(f2, k)) / mu(k);
  endfor
  d *= d;

  for(k = 1:N)
    for(m = 1:N)
      den += (D(k, m)*B(f2, k)*BA(f2, m)) / (mu(k)*mu(m));
    endfor
  endfor
  d /= den;
end


Calculation(1,1,false);
window1 = figure('Name', 'Graphic result');
grid on;
hold on;
H = 1:16;
for(f2p = 1:7)
  for(f1p = 1:16)
    G(f1p) = Calculation(f1p, f2p, true, 16, 0.002, 1, 3, 1000);
  endfor
  plot(H, G);
endfor