From 3af0ed0421dd8189cf2dc1d5c8dc9970e40b6e39 Mon Sep 17 00:00:00 2001
From: Yarmo Mackenbach <yarmo@protonmail.com>
Date: Wed, 23 Jan 2019 09:54:29 +0100
Subject: [PATCH] Initial commit

---
 examples.m   |  53 +++++++++++++++++++
 fourierViz.m | 144 +++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 197 insertions(+)
 create mode 100644 examples.m
 create mode 100644 fourierViz.m

diff --git a/examples.m b/examples.m
new file mode 100644
index 0000000..63d22c1
--- /dev/null
+++ b/examples.m
@@ -0,0 +1,53 @@
+%% Sine
+
+frequency = 1;
+amplitude = 1;
+phase = 0; %pi/2;
+
+fourierViz(frequency, amplitude, phase, 'baseMultiply', 0.25, 'export', false);
+
+
+%% Square simple
+
+frequency = 1:2:3;
+amplitude = 2./(pi*comps);
+phase = [0 0.25]*2*pi; % Phases expressed in cycles
+
+fourierViz(frequency, amplitude, phase, 'baseMultiply', 1, 'export', false);
+
+
+%% Square
+
+frequency = 1:2:41;
+amplitude = 2./(pi*frequency);
+phase = zeros(numel(frequency), 1);
+
+fourierViz(frequency, amplitude, phase, 'baseMultiply', 4, 'export', false);
+
+
+%% Square complex
+
+frequency = 1:2:101;
+amplitude = 2./(pi*frequency);
+phase = zeros(numel(frequency), 1);
+
+fourierViz(frequency, amplitude, phase, 'baseMultiply', 20, 'export', false);
+
+
+%% Sawtooth
+
+frequency = 1:15;
+amplitude = 2./(pi*frequency);
+phase = zeros(numel(frequency),1);
+
+fourierViz(frequency, amplitude, phase, 'baseMultiply', 3, 'export', false);
+
+
+%% Beat
+
+frequency = [1 1.2];
+amplitude = [1 1];
+phase = [0 0];
+
+fourierViz(frequency, amplitude, phase, 'baseMultiply', 1.5, 'export', false);
+
diff --git a/fourierViz.m b/fourierViz.m
new file mode 100644
index 0000000..2e321f1
--- /dev/null
+++ b/fourierViz.m
@@ -0,0 +1,144 @@
+function fourierViz(varargin)
+%FOURIERVIZ Fourier visualization in Matlab
+% 
+%   By default, the fastest component has a frequency of 1 Hz
+
+    % Handle input
+    p = inputParser;
+    p.addOptional('frequency', 1);
+    p.addOptional('amplitude', 1);
+    p.addOptional('phase', 0);
+    p.addParamValue('baseMultiply', 2);
+    p.addParamValue('export', false);
+    p.addParamValue('filename', 'fourierVizLastRun.gif', @ischar);
+    p.parse(varargin{:});
+    p = p.Results;
+    
+    % Component frequency scale factor
+    freqScaleFactor = max(p.frequency) / p.baseMultiply;
+    
+    % Normalize magnitudes
+    p.amplitude = p.amplitude / sum(p.amplitude);
+    
+    % Figure
+    doPlot = true;
+    fh = figure('Position', [200 100 1600 800]);
+    opts = struct;
+    opts.margin = [0 0 0 0];
+    opts.spaceBetweenPlots = 0;
+    opts.nX = 1;
+    opts.nY = 1;
+    opts.iX = 1;
+    opts.iY = 1; 
+    mk.subplot(opts);
+    hold on;
+    set(gca,'xtick',[]);
+    set(gca,'xticklabel',[]);
+    set(gca,'ytick',[]);
+    set(gca,'yticklabel',[]);
+    xlim([-1.2 3.2]);
+    ylim([-1.2 1.2]);
+    
+    % Time
+    dt = 0.016;
+    tick = 0;
+    
+    % Base frequency & period
+    baseFrequency = sym(p.frequency(1)/freqScaleFactor);
+    for ii = 2:numel(p.frequency)
+        baseFrequency = gcd(baseFrequency, sym(p.frequency(ii)/freqScaleFactor));
+    end
+    baseFrequency = double(baseFrequency);
+    basePeriod = 1/baseFrequency;
+    
+    % Circle
+    xCircle = 0:2*pi/100:2*pi;
+    
+    % Oscilloscope
+    osc = nan(6e2,1);
+    oscilloscopePrefilled = false;
+    
+    % Colors
+    clrGrey = struct('color', [0.8784 0.8784 0.8784]);
+    clrGreen = struct('color', [0.1804 0.4902 0.1961]);
+    clrBlue = struct('color', [0.0824 0.3961 0.7529]);
+    clrRed = struct('color', [0.7765 0.1569 0.1569]);
+    
+    while doPlot
+        tick = tick + 1;
+        
+        % Has the figure been closed
+        if ~ishandle(fh)
+            doPlot = false;
+            continue;
+        end
+        
+        % Update time
+        t = dt * tick;
+        
+        % Update oscilloscopePrefilled
+        if ~oscilloscopePrefilled && t>=basePeriod
+            if ~any(isnan(osc))
+                oscilloscopePrefilled = true;
+            end
+            tick = 1;
+            t = dt * tick;
+        end
+        
+        % Init plot coordinates
+        x = 0;
+        y = 0;
+        
+        % Plot
+        if oscilloscopePrefilled
+            cla;
+            plot(cos(xCircle), sin(xCircle), '-k', clrGrey);
+        end
+        
+        % Compute components
+        for ii = 1:numel(p.frequency)
+            if mod(ii,2) == 0
+                clr = clrGreen;
+            else
+                clr = clrBlue;
+            end
+            alpha = ((t * 2*pi * p.frequency(ii)) / freqScaleFactor) + p.phase(ii);
+            x2 = x + (p.amplitude(ii)*cos(alpha));
+            y2 = y + (p.amplitude(ii)*sin(alpha));
+            if oscilloscopePrefilled
+                plot([x x2], [y y2], '-k', clr, 'LineWidth', 2);
+                plot(x+cos(xCircle)*p.amplitude(ii), y+sin(xCircle)*p.amplitude(ii), '-k', clr, 'LineWidth', 2);
+            end
+            x = x2;
+            y = y2;
+        end
+        osc = [y; osc(1:end-1)];
+        if oscilloscopePrefilled
+            plot(x, y, 'or', clrRed);
+            plot([x, 1.5], [y y], '-r', clrRed, 'LineWidth', 1);
+            plot(linspace(1.5, 3.2, numel(osc)), osc, '-r', clrRed, 'LineWidth', 2);
+        end
+        
+        % Export to gif
+        if oscilloscopePrefilled && p.export && t < basePeriod
+            if ishandle(fh)
+                frame = getframe(fh);
+                im = frame2im(frame);
+            end
+            [imind, cm] = rgb2ind(im, 256);
+
+            if tick == 1
+                imwrite(imind, cm, p.filename, 'gif', 'Loopcount', inf, 'DelayTime', 0.016); 
+            else
+                imwrite(imind, cm, p.filename, 'gif', 'WriteMode', 'append', 'DelayTime', 0.016); 
+            end
+        end
+        
+        % Pause time
+        if oscilloscopePrefilled
+            pause(dt);
+        end
+    end
+
+end
+