Fourier series animation

README.md

MATLAB code to demonstrate Fourier series representation of periodic signals (as a sum of sinusoidal functions).

The animation shows an approximation of a square wave signal using the first 4-terms of its Fourier series. (Change the parameters near the top of the code to manipulate the animations and explore other variations).

This was inspired by the following similar animations:

• Fourier Series Animation using Circles
• Fourier series square wave circles animation
• Fourier series visualisation with d3.js
• Smooth motion of rotating circles

animation1.m

% animation params
periods = 2;                % how many full-period rotations to do
num_steps = 100 * periods;  % number of animation steps per period

% circles
freq = [1 3 5 7];           % frequency of each sine function (rotation speed)
r = 4./(freq*pi);           % radii of circles (amplitude of sines)
K = numel(freq);            % number of circles
R = max(r);                 % maximum radius
bounds = ceil(R);           % offset the circles to left of origin

% circles points (co-centric)
t = linspace(0, 2*pi, 50).';
x = bsxfun(@times, r, cos(t)) - bounds;  % shifted with offset
y = bsxfun(@times, r, sin(t));

% PLOT: open wide figure
pos = get(groot, 'DefaultFigurePosition');
hFig = figure('Position',pos .* [1 1 1.5 0.9]);
movegui(hFig, 'center')

% PLOT: circles
hCircles = line(x, y, 'LineWidth',2);       % plot circles
line([0 0], [-bounds bounds], 'Color','k')  % plot y=0 axis
grid on, box on
axis equal
axis([-2*bounds ceil(2*pi*periods) -bounds bounds])
xlabel('Angle \theta')
title('Sine Functions')
ax = gca;
ax.XTick = ax.XTick(ax.XTick >= 0);         % start x-ticks at zero

% PLOT: prepare animated graphics objects
hRadii = line(nan(2,K), nan(2,K), 'LineWidth',1);
hLines = line(nan(2,K), nan(2,K), 'LineWidth',1, ...
    'LineStyle','--', 'Marker','.', 'MarkerSize',15);

% PLOT: sine functions curves
xx = linspace(0, 2*pi*periods, num_steps);  % X-periods in N-steps
hCurves = gobjects(K,1);
for k=1:K
    hCurves(k) = animatedline('Color',hRadii(k).Color, 'LineWidth',2);
end
legend(hCurves, ...
    strcat('4sin(', num2str(freq(:)), '\theta)/', num2str(freq(:)), '\pi'), ...
    'Orientation','Horizontal', 'Location','SouthOutside')

% GIF
%{
% gifsicle --colors 256 --careful --delay=5 --optimize --crop 90,45+675x260 -o out1_.gif out1.gif
[im,map] = rgb2ind(frame2im(getframe(hFig)), 256);
imwrite(im, map, 'out1.gif', 'gif', 'DelayTime',0.1, 'LoopCount',Inf);
%}

% animation
theta = zeros(1,K);            % current rotation angles for each circle (in radians)
step = (xx(2)-xx(1)) .* freq;  % step sizes for each circle (arclen = angle / radius)
for i=1:num_steps
    % break if figure is closed
    if ~isgraphics(hFig), break; end
    
    % for each circle
    for k=1:K
        % rotating point on k-th circle
        x = r(k) * cos(theta(k));
        y = r(k) * sin(theta(k));

        % update graphics
        set(hLines(k), 'XData',[x-bounds; xx(i)], 'YData',[y; y])
        set(hRadii(k), 'XData',[0; x]-bounds, 'YData',[0; y])
        addpoints(hCurves(k), xx(i), y)
    end

    % refresh plot
    pause(0.02)    % drawnow
    
    % GIF
    %{
    [im,map] = rgb2ind(frame2im(getframe(hFig)), 256);
    imwrite(im, map, 'out1.gif', 'gif', 'DelayTime',0.1, 'WriteMode','append');
    %}

    % increment angles
    theta = theta + step;
end

animation2.m

% animation params
periods = 2;                % how many full-period rotations to do
num_steps = 150 * periods;  % number of animation steps per period

% circles
freq = [1 3 5 7];           % frequency of each sine function (rotation speed)
r = 4./(freq*pi);           % radii of circles (amplitude of sines)
K = numel(freq);            % number of circles
R = max(r);                 % maximum radius
bounds = ceil(sum(r));      % offset the circles to left of origin

% circles points (co-centric)
t = linspace(0, 2*pi, 50).';  % NPTS=50
x = bsxfun(@times, r, cos(t));
y = bsxfun(@times, r, sin(t));

% add center (so that the radius line is drawn)
x = x([1 1:end],:); x(1,:) = 0;
y = y([1 1:end],:); y(1,:) = 0;

% circle points in homogeneous form (one slice per circle)
% (size = 4 x NUMPTS x K)
pts = cat(2, reshape(x,[],1,K), reshape(y,[],1,K));
pts(:,3,:) = 0;               % z=0
pts(:,4,:) = 1;               % w=1
pts = permute(pts, [2 1 3]);  % transpose

% PLOT: open wide figure
pos = get(groot, 'DefaultFigurePosition');
hFig = figure('Position',pos .* [1 1 1.5 0.9]);
movegui(hFig, 'center')

% PLOT: circles with radius lines (initially untransformed and all centered
% on origin), then translate circles to the left
hCircles = line(x, y, 'LineWidth',2);
hTr = hgtransform('Matrix',makehgtform('translate',[-bounds 0 0]));
set(hCircles, 'Parent',hTr)
line([0 0], [-bounds bounds], 'Color','k')  % plot y=0 axis
grid on, box on
axis equal
axis([-2*bounds ceil(2*pi*periods) -bounds bounds])
xlabel('Angle \theta')
ax = gca;
ax.XTick = ax.XTick(ax.XTick >= 0);         % start x-ticks at zero

% PLOT: trace path from smallest circle
hPath = animatedline('Color',hCircles(end).Color, 'LineStyle','-', ...
    'MaximumNumPoints',ceil(num_steps/periods)+1, 'Parent',hTr);

% PLOT: connecting line between trace and curve
hLine = line(nan, nan, 'Color',[0.5 0.5 0.5], ...
     'LineStyle','-', 'Marker','.', 'MarkerSize',10);

% PLOT: resulting Fourier series (sum of harmonics)
xx = linspace(0, 2*pi*periods, num_steps);  % X-periods in N-steps
yy = nan(size(xx));
hCurve = line(xx, yy, 'Color','m');
labels = cellstr(strcat('4sin(', num2str(freq(:)), '\theta)/', num2str(freq(:)), '\pi'));
legend(hCurve, strjoin(labels,' + '), ...
    'Orientation','Horizontal', 'Location','SouthOutside')

% GIF
%{
% gifsicle --colors 256 --careful --delay=5 --optimize --crop 90,30+675x300 -o out2_.gif out2.gif
[im,map] = rgb2ind(frame2im(getframe(hFig)), 256);
imwrite(im, map, 'out2.gif', 'gif', 'DelayTime',0.1, 'LoopCount',Inf);
%}

% animation
theta = zeros(1,K);                     % current rotation angles for each circle (in radians)
step = (2*pi*periods).*freq/num_steps;  % step sizes for each circle
for i=1:num_steps+1
    % break if figure is closed
    if ~isgraphics(hFig), break; end

    % compute and store transformation matrices (translations and rotations)
    tt = [theta(1) diff(theta)];
    Tx = cell(1,K);
    Rz = cell(1,K);
    for k=1:K
        Tx{k} = makehgtform('translate',[r(k) 0 0]);
        Rz{k} = makehgtform('zrotate',tt(k));
    end

    % for each circle
    for k=1:K
        % chain of transformations
        T = Rz{1};
        for j=1:k-1
            T = T * Tx{j} * Rz{j+1};
        end

        % apply combined transform on k-th circle
        pt = T * pts(:,:,k);
        pt = pt(1:2,:);  % extract x/y coords from homogeneous form

        % update k-th circle
        set(hCircles(k), 'XData',pt(1,:), 'YData',pt(2,:))
    end

    % update trace line
    pt = pt(:,2);  % take first point from last circle (skipping the center)
    addpoints(hPath, pt(1), pt(2))

    % update curve line
    yy = [pt(2) yy(1:end-1)];  % prepend point, and rollover
    set(hCurve, 'YData',yy)    % or fliplr(yy)

    % update connecting horizontal line
    set(hLine, 'XData',[pt(1)-bounds; 0], 'YData',[pt(2); pt(2)])

    % refresh plot
    pause(0.02)    % drawnow

    % GIF
    %{
    [im,map] = rgb2ind(frame2im(getframe(hFig)), 256);
    imwrite(im, map, 'out2.gif', 'gif', 'DelayTime',0.1, 'WriteMode','append');
    %}

    % increment angles
    theta = theta + step;
end

out1_.gif

out2_.gif

any reason I'd be getting the following error:

Attempt to reference field of non-structure array.

Error in animation1 (line 34)
ax.XTick = ax.XTick(ax.XTick >= 0); % start x-ticks at zero

@JMD8881:
Right, this syntax is from the new HG2 graphics systems introduced in the latest MATLAB R2014b.
I'm also using new functions like animatedline

You could replace that part with the old syntax of get(gca,..) and set(gca,...). But you'll also have to replace the use of animatedline with regular line graphics object and update them manually...

ok, thanks! I'll just try it again when we update to the newest version.

hey ! can you tell me Is fourier series is a transform or not .. as i know that fourier series analyze time domain signal in frequency domain so it is a transform am i right?

You really should make this more visible than just a gist. I come to it over and over, and always find it hard to find :-)

This is a great animation. I can't seem to get the same GIF 'quality' that you get. How did you make such a nice GIF?

Does anyone know if there is something like this but in python? I don't have matlab
Thanks

This is a inspiring animation. However, what would happen when two wave functions intersect in a Fourier series representation of periodic signals?