%% Time-Varying Systems and Computations
% Klaus Diepold, Sven Gronauer, Simon Schumann
%
%% Convolution, Eigenvalue Decomposition and Fourier Transformation
%% 1. Linear Convolution
% Create column vectors u and v, both of length n. Take for now n=4.

N = 4;
u = randi([-10, 10], N, 1)
v = randi([-10, 10], N, 1)

% We consider the convolution of two vectors u and v. To do this we can use the convolution operation in Matlab.

t = conv(u,v)

% Alternatively, we can represent the convolution by a matrix-vector multiplication using Teoplitz matrix.

T = toeplitz([u;zeros(length(v)-1,1)],[u(1);zeros(length(v)-1,1)])
t_conv = T*v

% Finally, we can compare our solution with the Matlab operation.

norm(t-t_conv)^2

%% 2. Cyclic Convolution
% Cyclic convolution relates to the process where the sequences are wrapped around in a cycle. This type of convolution is used where signals are considered periodic.

t_cyclic = cconv(u, v)

% Alternatively, we can perform cyclic convolution by first converting one of the vectors into a circulant matrix and then performing matrix-vector multiplication.

n = N + N - 1;
% Zero pad u and v to the same length n
u_padded = [u; zeros(n - length(u), 1)]
v_padded = [v; zeros(n - length(v), 1)]

% Creating the circulant matrix from the zero-padded vector v

tc_padded = toeplitz([v_padded(1); flipud(v_padded(2:end))], v_padded)'
t_cyclic_manual = tc_padded * u_padded

% Verify, that the cyclic convolution result t_cyclic_conv matches t_cyclic computed with cconv.

norm(t_cyclic-t_cyclic_manual)^2

%
%% 3. Eigenvalue Decomposition of Circulant Matrices
% Exploring the eigenvalue decomposition of circulant matrices using Vandermonde matrices.
%% 3.1 Special Vandermonde Matrix
% The matrix Q represents a special Vandermonde matrix where each column is a power of the N-th root of unity.

q = exp(2 * pi * 1i / N);
Q = zeros(N);
for k = 0:N-1
    Q(:,k+1) = q.^(k*(0:N-1)) / sqrt(N);
end
Q

% This matrix Q can be used to compute the DFT of a vector v by matrix multiplication:
% V=Q\cdot v
% Here, the vector X is the DFT of the vector x, and it corresponds to the equation
% V_k =\sum_{n=0}^{N-1} v_n \cdot e^{-i\frac{2\Pi }{N}\textrm{kn}},
% where each row of the matrix Q corresponds to k=0,1,2,\ldotp \ldotp \ldotp
% Next, let's illustrate the elements as a pointer diagram.

plotRowAsArrows(Q, 1); % Plot the first row
plotRowAsArrows(Q, 2); % Plot the second row
plotRowAsArrows(Q, 3); % Plot the third row

%% 3.2 Eigenvalue Decomposition of Circulant Toeplitz Matrix
% Construct the circulant matrix for vector v.

Tc = toeplitz([v(1); flipud(v(2:end))], v)'; 

% Perform the Eigenvalue decomposition

Lambda = diag(Q' * Tc * Q)

% Check if the reconstructed matric matches the original, in particalur proof that Q is orthogonal.

Tc_reconstructed = Q * diag(Lambda) * Q';
Tc
Tc_reconstructed

% Compare the original and reconstructed matrices to verify the accuracy of the decomposition.

norm(Tc - Tc_reconstructed, 'fro')

%% 3.3 Comparsion to FFT Eigenvalue Decomposition and Build-In Eigenvalue Decomposition
% Now look at the results a FFT produces for the vector v. Looks familiar, or?

Lambda_fft = fft(v)

% It is the same result as the build-in cholesky method for eigenvalue decomposition.

Lambda_chol = eig(Tc)

% Finally compare the computational demand for both methods.

tic;
Lambda = diag(Q' * Tc * Q);
time_eig = toc;

tic;
Lambda_fft = fft(v);
time_fft = toc;

% Display results
fprintf('Time for eigenvalue decomposition: %f seconds\n', time_eig);
fprintf('Time for FFT method: %f seconds\n', time_fft);

% Why is the fft typically faster than the Toeplitz technique?
%
%% 4. Discrete Fourier Transformation and Linear Convolution
% Have a second look at the linear convolution.

y_time = conv(u, v)

% Now consider it in the fourier space.

U = fft(u, n); 
V = fft(v, n); 

% Perfom elementwise multiplication.

Y = U .* V; 

% And transform back into the time domain.

y_freq = ifft(Y)

% Lastly check if the results are the same.

norm(y_time - y_freq, 'fro')

% Exercise: Perform a computation comparion between linear convolution and the presented FFT method for multiple N.
%
%% Appendix

function plotRowAsArrows(Q, rowIndex)
    N = size(Q, 2); 
    realParts = real(Q(rowIndex, :));
    imagParts = imag(Q(rowIndex, :));

    % Define a colormap
    colors = cool(N); 

    figure;
    hold on; 
    for i = 1:N
        quiver(0, 0, realParts(i), imagParts(i), 'AutoScale', 'off', ...
            'Color', colors(i,:), 'MaxHeadSize', 0.6);
    end

    % Plot unit circle
    theta = linspace(0, 2 * pi, 100);
    unitCircleX = cos(theta);
    unitCircleY = sin(theta);
    plot(unitCircleX, unitCircleY, 'k--'); 

    hold off;

    title(sprintf('Elements of Row %d of Q with Unit Circle', rowIndex));
    xlabel('Real Part');
    ylabel('Imaginary Part');
    grid on;
    axis equal;
    maxLimit = max([abs(realParts), abs(imagParts)]);
    xlim([-maxLimit, maxLimit]);
    ylim([-maxLimit, maxLimit]);
    xticks(linspace(-maxLimit, maxLimit, 2 * maxLimit + 1));
    yticks(linspace(-maxLimit, maxLimit, 2 * maxLimit + 1));
end