Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Fourier series animation

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).

animation1 animation2

This was inspired by the following similar animations:

% 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
% 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
@talk2laton

This comment has been minimized.

Copy link

@talk2laton talk2laton commented Jan 31, 2015

I love that. I will create mine too

@JMD8881

This comment has been minimized.

Copy link

@JMD8881 JMD8881 commented Feb 2, 2015

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

@amroamroamro

This comment has been minimized.

Copy link
Owner Author

@amroamroamro amroamroamro commented Feb 2, 2015

@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...

@JMD8881

This comment has been minimized.

Copy link

@JMD8881 JMD8881 commented Feb 3, 2015

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

@Adnankhalil

This comment has been minimized.

Copy link

@Adnankhalil Adnankhalil commented Sep 11, 2015

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?

@lmendo

This comment has been minimized.

Copy link

@lmendo lmendo commented Oct 28, 2016

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 :-)

@benjaminranderis

This comment has been minimized.

Copy link

@benjaminranderis benjaminranderis commented Sep 14, 2018

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?

@ncclementi

This comment has been minimized.

Copy link

@ncclementi ncclementi commented Mar 22, 2019

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment