Fourier Series Approximations of Common Waveforms

Seth DeLand
/
Fourier Series Approximations of Common Waveforms
on 8 Nov 2023
19
46
0
0
1948
function drawframe(nn)
% Approximate a triangle, square, and sawtooth wave using Fourier Series. 
% Number of Coefficients
o = 8;
l = lines;
% Triangle Wave
w(1).o = 1:2:o; 
w(1).c = (8/pi^2).*((-1).^((w(1).o-1)/2))./(w(1).o.^2);
w(1).a = @(t) 2*w(1).o*pi*t;
w(1).x0 = sum(w(1).c(1:end-1));
w(1).rbg = l(1,:);
% Square Wave
w(2).o = 1:o;
w(2).c = (4/pi)*1./(2*w(2).o-1); 
w(2).a = @(t) (2*w(2).o - 1)*2*pi*t;
w(2).x0 = sum(w(2).c(1:end-1));
w(2).rbg = l(2,:);
% Sawtooth Wave
w(3).o = 1:o;
w(3).c = (2/pi).*(1./w(3).o);
w(3).a = @(t) 2*w(3).o*pi*t;
w(3).x0 = sum(w(3).c(1:end-1));
w(3).rbg = l(3,:);
dx = 0.2;
%% Graphics Initialization
% Create the axes and title
clf;
yl = [-2 2];
xl = [-3 12];
n = 3;
tl = tiledlayout(3,1,'TileSpacing','tight');
ti = ["Triangle","Square","Sawtooth"];
for i = 1:n
    ax(i) = nexttile;
    set(ax(i),"YLim",yl,"XLim",xl,...
        "XAxisLocation","origin","YAxisLocation","origin",...
        "XTickLabel","","YTickLabel","","TickLength",[0 0],...
        "XGrid","on","YGrid","on",...
        "DataAspectRatio",[1 1 1],"PlotBoxAspectRatioMode","manual");
    title(ti(i))
    
    % Initialize an array of hgtransform objects.  We will use these transforms
    % to position the circles that approximate the Fourier series.  
    h(i).h = repelem(hgtransform,length(w(i).o),1);
    
    % Create the circles and transform them to the appropriate locations.
    ct(i).t = zeros(4,4,length(w(i).o));
    for n = 1:length(w(i).o)
        c = circle(0,0,w(i).c(n));
        if n==1 
            % The first circle always stays at the origin, so it gets an
            % identity translation.
            h(i).h(n) = hgtransform("Parent",ax(i));
            ct(i).t(:,:,n) = makehgtform("translate",[0 0 0]);
        else
            % Every other circle will need to be translated by the radius of
            % the circle before it.  We parent each circle to its predecessor.
            h(i).h(n) = hgtransform("Parent",h(i).h(n-1));
            ct(i).t(:,:,n) = makehgtform("translate",[w(i).c(n-1) 0 0]);
        end
        % Parent the circle to the hgtransform, and apply the translation.
        set(c,'Parent',h(i).h(n));
        set(h(i).h(n),"Matrix",ct(i).t(:,:,n));
    end
    
    % Initialize the line for the approximated signal.
    a(i).l = line(w(i).x0,0,'LineStyle','-','LineWidth',2,'Color',w(i).rbg);
    % Initialize the horizontal line that connects the last circle with the
    % approximated signal.
    z(i).l = line([w(i).x0 w(i).x0],[0 0],'LineStyle','--','Color',[0.5 0.5 0.5],'Marker','.','MarkerSize',10);
end
%% Animation
% Animation is done with an outer loop that loops over time and an inner
% loop that calculates the successive terms in the Fourier series.
for t = 0:nn-1 % Loop through time
    for i = 1:3 % Loop through each wave
        th = w(i).a(t/24);
        for n = 1:length(th) % Loop through term in Fourier series      
            if n>1
                % Our function calculates absolute angles, but because our
                % transforms work relative to each other we need to subtract
                % the previous angle for terms of order>1.
                cr = makehgtform("zrotate",th(n)-th(n-1));
            else
                cr = makehgtform("zrotate",th(n));
            end
            % Apply the rotation and then the translation. Transform order is
            % right-to-left of matrix multiplication order.
            set(h(i).h(n),"Matrix",ct(i).t(:,:,n)*cr);
        end
        % Calculate the x- and y-coordinate of the point on the last circle.
        % There is probably a way to calculate this from the transform objects,
        % but I find it easier to calculate it from the terms in the series.
        % This is a matrix multiply: coeffs is [1 x order] and theta' is 
        % [order x 1].
        x = w(i).c*cos(th)';
        y = w(i).c*sin(th)';
    
        % Append the newest point to the approximated line
        a(i).l.XData = [a(i).l.XData a(i).l.XData(end)+dx];
        a(i).l.YData = [y a(i).l.YData];
    
        % Update the coordinates of the horizontal line
        z(i).l.XData(1) = x;
        z(i).l.YData = [y y];
    end
    
    % Update the figure
end
end
%% Helper Functions
function c = circle(xc,yc,r)
    % Create a circle with center {xc,yc} with radius r. The circle also
    % has a radial line drawn from its center to {xc,0}.
    nPts = 1E2;
    th = linspace(0,2*pi,nPts);
    g = [0.5 0.5 0.5];
    x = xc+r*sin(th);
    y = yc+r*cos(th);
    c(1) = line(x,y,'Color',g);
    c(2) = line([xc xc+r],[yc yc],'Color',g,'Marker','.','MarkerSize',10,'LineStyle',':');
end
Animation

Remix Tree

Fourier Series Approximations of Common Waveforms
Seth DeLand