% Plot of some kind of input image via its Fourier components. See also
% https://www.dynamicmath.xyz/fourier-epicycles/
% https://mathed.miamioh.edu/index.php/ggbj/article/view/197/159
% and https://www.youtube.com/watch?v=r6sGWTCMz2k&t=1213s

clear; close all
figure('units','pixels','position',[0 0 3840 2160],'ToolBar','none');
vid = VideoWriter('discreteFT.mp4','MPEG-4');
vid.FrameRate = 30;    % Default 30
vid.Quality = 100;    % Default 75
open(vid);
set(gcf,'Color','k'); % Black background

pos1 = [0 0 1 1]; % Use up all of the screen

% Generate or read in complex numbers representing input image. Can be done
% quite nicely with GRABIT.m
% https://ch.mathworks.com/matlabcentral/fileexchange/7173-grabit

th = 0:pi/90:2*pi; % Angles defining circles drawn for each rotating Fourier component

freqMax = 350; % Maximum frequency in terms of base frequency (1 Hz). 
% Frequencies span from -freqMax to +freqMax

S = load('matterhorn001.txt'); % Input image. Can be created using GRABIT
x = S(:,1);
y = S(:,2);
N = 2*freqMax; % Number of evenly-spaced points defining 2D image (join the dots!)
% Now interpolate input image to make dots evenly spaced along the curve 
% Download https://ch.mathworks.com/matlabcentral/fileexchange/34874-interparc
evenArc = interparc(N,x,y,'linear'); % Take original input data and interpolate points evenly
xR = evenArc(:,1); % x-axis real coordinates
xI = evenArc(:,2); % y-axis imaginary coordinates


dt = 1/N; % Discrete sampling in time equal to reciprocal of number of dots making up the input image

% Custom rainbow colour for the circles representing the cn-amplitudes 
rgb = customcolormap([0 0.17 0.33 0.5 0.67 0.83 1],...
    {'#ff0000', '#ffa500', '#ffff00', '#00ff00',...
    '#00ffff', '#0000ff', '#7f00ff'},2*freqMax+1);
rgb = flip(rgb);

cn = zeros(2*freqMax+1,1); % Initialize Fourier-component complex amplitudes cn
% Generate Fourier-component complex amplitudes cn
for n = -freqMax:freqMax % Loop to calculate all cn's for all sampled frequencies 
    for jj = 1:N % For each cn, sum instead of integrating - the discrete part! 
        % Euler factor exp(-2pi i n t), where n is an integer from -freqMax
        % to +freqMax, and t = jj*dt goes from 0 to 1 in steps of dt = 1/N
        eulerFac = exp(-1i*n*2*pi*jj*dt);
        cn(n+1+freqMax,1) = cn(n+1+freqMax,1) + (1/(2*freqMax+1))*eulerFac*(xR(jj,1) + 1i*xI(jj,1)); % c(1) corresponds to c_-freqMax, etc...
    end
end

abscn = abs(cn); % Magnitudes of complex cn's
[abscnSort,I] = sort(abscn); % Order abscn and get the indices I
cnSort = flip(cn(I)); % Flip vector to go from large to small
Isort = flip(I); % List of indices from largest cn to smallest
nSort = Isort - freqMax - 1; % Convert index numbers (1 to 2*freqMax+1) to frequency multiple, n (-freqMax to +freqMax) 

% Now generate figure over one cycle
xy = zeros(N+1,1); % Initialize complex real-space values that will be generated by the inverse FT 
kk = 0;
for tt = 0:dt:1 % Go through one cycle in time. Each Fourier component will rotate n times in this time
    kk = kk+1; % Index from 1 to N
    for n = -freqMax:freqMax
        xy(kk,1) = xy(kk,1) + cn(n+1+freqMax,1)*exp(1i*2*pi*n*tt); % Told ya! 
    end
end

kk = 0;
% Do plotting part separately for the sake of cleanness
% Part 1: Draw out the image over one cycle, including the circles and
% vectors

A = 2*freqMax*1.0005;

for tt = 0:dt:1
    hold off
    kk = kk+1;
    plot(real(xy(1:kk)),imag(xy(1:kk)),'color','w','linewidth',4)
    hold on

    for jj = 1:2*freqMax+1 % Plot out all jj circles and vectors 
        % colIndex generated to rush through the reds and oranges and dwell
        % on the greens, blues, and violets, as these are the colours of
        % the tinsywinsy circles and can hardly be seen otherwise 
        colIndex = floor(A*(1 - exp(-7*jj/freqMax)));
        if (jj==1)
            centrex = 0;
            centrey = 0;
        else
            centrex = centrex + real(cnSort(jj-1)*exp(1i*2*pi*nSort(jj-1)*tt));
            centrey = centrey + imag(cnSort(jj-1)*exp(1i*2*pi*nSort(jj-1)*tt));
        end
        % Plot the Fourier component (amplitude AND phase) 
        plot([centrex centrex+real(cnSort(jj)*exp(1i*2*pi*nSort(jj)*tt))],...
            [centrey centrey+imag(cnSort(jj)*exp(1i*2*pi*nSort(jj)*tt))],...
            'color','w','linewidth',1)

        radius = abs(cnSort(jj));
        % Plot the circle on which the Fourier component rotates 
        if (radius>=0.001)
            plot(radius*cos(th)+centrex,radius*sin(th)+centrey,'color',rgb(colIndex,:),'linewidth',0.8)
        else
            plot(centrex, centrey, 'o', 'MarkerEdgeColor', 'none', 'MarkerFaceColor',rgb(colIndex,:),'MarkerSize',0.8);
        end
    end
    axis equal
    axis off 
    % In GRABIT, I tend to define the x- and y-limits as -1 to +1 (square) 
    xlim([-1 1]);
    ylim([-1 1]);
    % Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame);
end

% % Part 2: Fade out circles
for fade = 1:-1/30:0
    hold off
    plot(real(xy(1:N+1)),imag(xy(1:N+1)),'color','w','linewidth',4)
    hold on
    for jj = 1:2*freqMax+1 % Plot out all jj circles and vectors
        colIndex = floor(A*(1 - exp(-7*jj/freqMax)));
        if (jj==1)
            centrex = 0;
            centrey = 0;
        else
            centrex = centrex + real(cnSort(jj-1)*exp(1i*2*pi*nSort(jj-1)*tt));
            centrey = centrey + imag(cnSort(jj-1)*exp(1i*2*pi*nSort(jj-1)*tt));
        end
        vector = plot([centrex centrex+real(cnSort(jj)*exp(1i*2*pi*nSort(jj)*tt))],...
            [centrey centrey+imag(cnSort(jj)*exp(1i*2*pi*nSort(jj)*tt))],...
            'color','w','linewidth',0.7);
        vector.Color(4) = fade;

        radius = abs(cnSort(jj));
        if (radius>=0.001)
            circ = plot(radius*cos(th)+centrex,radius*sin(th)+centrey,'color',rgb(colIndex,:),'linewidth',0.5);
            circ.Color(4) = fade;
        else
            spot = plot(centrex, centrey, 'o', 'MarkerEdgeColor', 'none', 'MarkerFaceColor',rgb(colIndex,:),'MarkerSize',0.5);
            spot.Color(4) = fade;
        end
    end
    axis equal
    axis off
    xlim([-1 1]);
    ylim([-1 1]);
    % Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame);
end

% Part 3: Subtract from high to low frequencies the Fourier components 
% until you have only a circle. Plot this for each stage of "simplification" 

% Speed up first part of getting rid of high-frequency components,
% as one can hardly see any change. Then slow down for the second half 
reduce1 = N+1:-2:N+1-floor(N/2); % First half
reduce2 = N-floor(N/2):-1:1; % Second half 
reduceTot = [reduce1 reduce2];
for reduce = reduceTot
    xy = zeros(N+1,1);
    kk = 0;
    for tt = 0:dt:1
        kk = kk+1; % Index from 1 to N+1
        for ll = 1:reduce
            xy(kk,1) = xy(kk,1) + cnSort(ll,1)*exp(1i*2*pi*nSort(ll)*tt);
        end
    end
    hold off
    plot(real(xy(1:N+1)),imag(xy(1:N+1)),'color','w','linewidth',4)
    axis equal
    axis off
    xlim([-1 1]);
    ylim([-1 1]);
    % Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame);
end

% Output the movie as an mp4 file and close it
close(vid); 
% Ta-da! 