fourierViz/fourierViz.m
2019-01-23 09:54:29 +01:00

145 lines
4.1 KiB
Matlab

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