clc; clear;

phi = 0.5 * pi;
N = 5000;
t = linspace(-2 * pi, 2 * pi, N)';
x = sin(2 * pi * t) + 3 * sin(3 * 2 * pi * t);
x_2 = sin(2 * pi * t + phi) + 3 * sin(3 * 2 * pi * t + phi);

f = fftshift(fft(x));
F = dftmtx(N);
f_hat = fftshift(F * x);

f_2 = fftshift(fft(x_2));

J = circshift(eye(N), 1, 2);
[J_eigenvector, J_eigenvalue] = eig(J);
J_eigenvalue = diag(J_eigenvalue);

CMatrix = toeplitz(x);
[CMatrix_eigenvector, CMatrix_eigenvalue] = eig(CMatrix);
CMatrix_eigenvalue = diag(CMatrix_eigenvalue);

CMatrix_eigenvalue_hat = zeros(size(CMatrix_eigenvalue));
for i = 1: N
    tmp = f(i) * J_eigenvalue(i) ^ (i - 1);
    CMatrix_eigenvalue_hat(i) = CMatrix_eigenvalue_hat(i) + tmp;
end


if 1 == 1
    subplot(311);
    plot(x);
    
    subplot(312);
    plot(abs(f));
    
    subplot(313);
    % plot(abs(f_hat));
    plot(abs(f_2));

    % polarplot(angle(J_eigenvalue),abs(J_eigenvalue),"o")
else
    subplot(311);
    plot(abs(f_hat));
    
    subplot(312);
    plot(sort(abs(CMatrix_eigenvalue)));

    subplot(313);
    plot(sort(abs(CMatrix_eigenvalue_hat)));
    
end