Fourier-analyse: Lama's tekenen met cirkels

Wat bedoel ik precies met tekenen? Het idee is dat we een pad kunnen nemen dat iets representeert, zoals een vis, het $\pi$-symbool, of vrijwel alles wat we kunnen tekenen door een potlood op papier te zetten en een afbeelding te schetsen zonder het potlood op te tillen tot de voltooiing. Met dit pad kunnen we een reeks vectoren die in cirkels draaien met verschillende frequenties aan elkaar koppelen (tip-aan-staart). De punt van de laatste vector tekent dan ons oorspronkelijke ontwerp, of in ieder geval iets dat er zeer sterk op lijkt.

Je kunt bijvoorbeeld een rechte lijn tekenen met behulp van een vorm die juist het tegenovergestelde is van een rechte lijn. Door de cirkels aan te passen, kun je de getekende vorm veranderen. Door bijvoorbeeld de straal van de buitenste cirkel te verkleinen en de draaisnelheid te verhogen, kunnen we een vierkant tekenen.

Wat gebeurt er dan als we een derde cirkel toevoegen? Hiermee kunnen de curves geavanceerder worden. Met een derde cirkel kunnen we een curve tekenen die lijkt op een vis. Voor een extremer voorbeeld van wat mogelijk is met enkel cirkels die aan elkaar zijn gekoppeld, is er een lama die wordt getekend met $1024$ cirkels, elk met verschillende frequenties, stralen en starthoeken. In het algemeen kan (bijna) elke gesloten curve worden getekend met cirkels, ook wel epicykels genoemd.

Over het algemeen geldt: hoe meer cirkels we toevoegen, hoe complexere tekeningen we kunnen produceren. Ook zal de tekening bij meer cirkels beter lijken op het oorspronkelijke "input"-ontwerp. Om te begrijpen wat die animatie creëert, moeten we kijken naar de onderliggende wiskunde en intuïtie achter de Fourier-transformatie (en -reeks). Zoals eerder vermeld, zijn er drie eigenschappen van een cirkel die we kunnen aanpassen om verschillende curves te tekenen:

• De snelheid waarmee hij draait (frequentie)
• Hoe groot de cirkel is (straal)
• De hoek waarmee hij begint (fase)

Al deze drie eigenschappen kunnen worden gerepresenteerd door één enkel complex getal.

Fourier-transformatie

Een Fourier-transformatie is een transformatie van het tijdsdomein naar het frequentiedomein. Wat betekent dit? Als we een functie hebben, $\sin(2\pi\times 3t)$, dan zou de Fourier-transformatie van die functie simpelweg een enkele piek (spike) zijn bij de frequentie $\pm 3\text{hz}$.

In essentie doen we het volgende: we zoeken welke frequenties aanwezig zijn in datgene wat we willen representeren, en vervolgens combineren we die sinusgolven tot een benadering van het oorspronkelijke object.

Complexe getallen

Een complex getal bestaat uit twee delen: het "reële" deel en het "imaginaire" deel. Dit voegt in feite een extra as toe aan de getallen die we al kennen om ze te kunnen representeren. Een complex getal wordt geschreven als $3+4i$, waarbij $3$ het reële deel is en $4i$ het imaginaire deel.

Een andere nuttige eigenschap is dat complexe getallen gemakkelijk kunnen worden beschouwd als tweedimensionale vectoren, waarbij het reële deel $x$ is en het imaginaire deel $y$ in een standaard punt $(x, y)$ in de tweedimensionale Euclidische ruimte. Deze representatie is zeer handig vanwege de connectie tussen trigonometrie en complexe exponentiële functies. De formule die de link tussen deze twee legt, staat bekend als de identiteit van Euler.

Lineaire algebra

Lineaire algebra is een belangrijk concept om te begrijpen hoe we afbeeldingen tekenen met vectoren die in cirkels draaien. Namelijk: hoe we vectoren kunnen representeren in termen van andere vectoren. In de standaard 2D $(x, y)$ Euclidische ruimte kunnen we alle vectoren representeren in termen van de eenheidsvectoren $\textbf{x}$ en $\textbf{y}$. Deze worden gewoonlijk $\hat{\textbf{i}}$ en $\hat{\textbf{j}}$ genoemd, met $\hat{\textbf{i}} = (1, 0)$ (de $x$-as) en $\hat{\textbf{j}} = (0, 1)$ (de $y$-as).

Om dit te doen gebruiken we twee operaties: vermenigvuldiging met een scalar en optelling. Bijvoorbeeld, de vector $\vec{a} = (5,-8)$ kan worden uitgedrukt als:

$$ \begin{aligned}\vec{a} &= 5\cdot \hat{\textbf{i}} - 8\cdot \hat{\textbf{j}} \\ &= 5\cdot (1, 0) + -8\cdot (0, 1)\\ &=(5, 0) + (0, -8)\\ &= (5, -8)\end{aligned} $$

Deze vermenigvuldiging met een scalar en optelling van vectoren om een andere vector te produceren, staat bekend als de lineaire combinatie van een set vectoren.

$\hat{\textbf{i}}$ en $\hat{\textbf{j}}$ zijn echter niet de enige vectoren die gebruikt kunnen worden. In feite kunnen elke twee vectoren worden geschaald en opgeteld om elke andere vector te vormen, zolang ze geen veelvoud van elkaar zijn. Vectoren die niet tot elkaar geschaald kunnen worden, zijn lineair onafhankelijk en kunnen dus basisvectoren zijn van de gehele 2D-ruimte.

Het inproduct (of dot product) is de methode waarmee we bepalen met welk getal vectoren moeten worden geschaald om ze in termen van de gekozen vectoren te representeren. Als we een vector $\vec{\textbf{w}}$ willen uitdrukken in termen van vectoren $\vec{u}$ en $\vec{v}$, moeten we bepalen welke getallen we ermee moeten vermenigvuldigen om de lineaire combinatie $\vec{w} = c1\cdot \vec{u}+c2\cdot \vec{v}$ te vormen. Een andere manier om hierover na te denken is dat het inproduct de vraag beantwoordt: "hoeveel van $\vec{u}$ zit er in $\vec{w}$?".

Het inproduct van twee willekeurige vectoren is gedefinieerd als:

$$ \begin{aligned} \langle\vec{v1},\vec{v2}\rangle &= {v1}x\times {v2}x + {v1}y\times {v2}y\\ &= \|\vec{v1}\|\times \|\vec{v2}\| \times \cos(\theta) \end{aligned} $$

Waarbij $\theta$ de hoek tussen $v1$ en $v2$ is en $\|\vec{v}\|$ de grootte (magnitude) van de vector is. De magnitude, ook wel de norm genoemd, is gedefinieerd als $\|\vec{v}\| = \sqrt{{v}x^2 + {v}y^2}$.

Het inproduct heeft enkele belangrijke eigenschappen:

1. Het inproduct van een vector met zichzelf is het kwadraat van de magnitude: $\vec{v1}\cdot \vec{v1} = \|\vec{v_1}\|^2$.
2. Het inproduct van vectoren die loodrecht op elkaar staan (orthogonaal), is $0$, omdat $\cos(\frac{\pi}{2}) = 0$.

Inproducten van continue functies

De definitie van het inproduct kan worden uitgebreid naar continue functies. Het belangrijkste verschil zit in het type sommatie. Voor een 2D-inproduct is de sommatie discreet over een eindig aantal termen:

$$ \langle\vec{a},\vec{b}\rangle = \sum{n =1}^{N}{an\times b_n} $$

Voor complexe vectoren van lengte $n$, $\mathbb{C}^n$, moeten we het complex geconjugeerde van de tweede vector nemen:

$$ \langle\vec{a},\vec{b}\rangle = \sum{n = 1}^{N}{an\times \overline{b_n}} $$

Voor sommaties over continue intervallen, waar het interval oneindig klein kan zijn, gebruiken we een integraal:

$$ \langle f(t),g(t)\rangle = \frac{1}{T}\int_0^T{f(t)\overline{g(t)}dt} $$

De enige restrictie is dat de functies $f(t)$ en $g(t)$ kwadratisch integreerbaar moeten zijn over het interval $[0,T]$. Dit betekent dat de integraal van het kwadraat van de absolute waarde eindig is:

$$ f(t) \in L^2(0,T)\quad\text{als}\quad\int_0^T{|f(t)|^2dt}\quad\text{eindig is.} $$

Noot: $L^2$ is de verzameling kwadratisch integreerbare functies en $L^2(0, T)$ is de verzameling kwadratisch integreerbare functies over het interval $[0, T]$.

Complexe exponentiële functies als sinus en cosinus

Complexe exponentiële functies bieden een handige manier om met sinussen en cosinussen om te gaan via de identiteit van Euler:

$$ e^{ix} = \cos(x) + i \times \sin(x) $$

Waarbij $x$ beschrijft hoe ver men rond een cirkel reist. Bijvoorbeeld, $e^{\pi i}$ gaat $\pi$ radialen rond de cirkel en eindigt bij $-1 + 0i$, wat leidt tot de beroemde identiteit: $e^{\pi i} = -1$.

Waarom zorgt het verheffen van $e$ tot een imaginair getal voor een rotatie rond een eenheidscirkel? Het antwoord ligt in de natuurkunde. Als de positie $p$ op tijdstip $t$ is $p(t)=e^{it}$, dan is de snelheid de eerste afgeleide: $v(t)=\frac{dp}{dt}=i\times e^{it}$.

Vermenigvuldigen met $i$ is gelijk aan een rotatie van $90^{\circ}$ tegen de klok in. Omdat de snelheidsfunctie altijd loodrecht staat op de richting van de huidige positie, is de curve die wordt getekend een cirkel.

Daarnaast geldt dat $\frac{1}{2\pi}\int_0^{2\pi}{\cos(t)\sin(t)\,dt}=0$, omdat $\cos(t)$ en $\sin(t)$ orthogonale functies zijn, wat ze een basis maakt voor de Fourier-transformatie.

In wezen voeren we een basisverandering uit van de tweedimensionale complexe getallenruimte naar een oneindig dimensionale ruimte met basisfuncties $(...,e^{i\frac{-4\pi}{T}t}, e^{i\frac{-2\pi}{T}t}, 1, e^{i\frac{2\pi}{T}t}, e^{i\frac{4\pi}{T}t},...)$. Omdat dit een lineaire operatie is die van de ene functieruimte naar de andere mapt, noemen we dit een lineaire operator.

Een blokgolf benaderen

De blokgolf is een soort "hello world" in de wereld van signaalverwerking. Een blokgolf wordt gedefinieerd als een $[-\pi, \pi]$ periodieke golf waarbij de helft van de periode "aan" is (1) en de andere helft "uit" (0):

$$ \text{square}(t)= \text{sgn}(\sin{(t)}) $$

Waarbij sgn de tekenfunctie is: $$ \text{sgn}(x)=\begin{cases}1 &\text{als}\; x>0 \\ 0 &\text{als}\; x=0 \\ -1 &\text{als}\; x<0 \end{cases} $$

We kunnen integreren (het inproduct van de blokgolf-functie en de $k$-de basisterm) om de $k$-de Fourier-coëfficiënt van de blokgolf te vinden:

$$ ck=\int{-\pi}^{\pi}{\text{sgn}(\sin (t))e^{-i k t}dt} $$

Omdat de tekenfunctie discontinu is, moeten we de integraal splitsen per geval. Na berekening komen we uit op de volgende stuksgewijze functie voor $c_k$:

$$ c_k = \begin{cases} 0 &\text{als}\; k=0 \\ \cfrac{i\left(e^{-ik\pi}+e^{ik\pi}-2\right)}{2k\pi} &\text{als}\;k\not = 0 \end{cases} $$

Om de blokgolf te benaderen, sommeren we elke $c_k$ coëfficiënt vermenigvuldigd met $e^{i k t}$. Na vereenvoudiging via trigonometrische identiteiten krijgen we een formule zonder imaginaire termen:

$$ \frac{4}{\pi}\sin(t) + \frac{4}{3\pi}\sin(3t) + \dots $$

Naarmate we meer sinusgolven toevoegen, nadert de benadering de oorspronkelijke blokgolf. Dit noemen we de Fourier-reeks.

Lama's tekenen

Dezelfde principes kunnen worden gebruikt om vrijwel alles te tekenen, zolang het gerepresenteerd kan worden als een gesloten curve die redelijk vloeiend is. Om een functie te verkrijgen die een lama representeert, kunnen we een SVG-bestand gebruiken. We bemonsteren het SVG-pad op gelijkmatig verdeelde intervallen om een lijst met samplepunten te krijgen, waarna we een Discrete Fourier Transformatie (DFT) uitvoeren.

Sampling van SVG-punten

Hieronder staat de implementatie voor het bemonsteren van een SVG naar punten in TypeScript:

const numSamples = 400;
let pts = [];
async function entry() {
  const svgUrl = './images/llama.svg';
  const response = await fetch(svgUrl);
  const text = await response.text();
  pts = await sampleSvgPoints(text, numSamples);
}
entry();

async function sampleSvgPoints(xml: string, numPts: number): Promise<[number, number][]> {
  const parser = new DOMParser();
  const doc = parser.parseFromString(xml, 'application/xml');
  let samplePts = [];
  const paths = Array.from(doc.querySelectorAll('path'));
  const totalLength = paths.reduce((acc, v) => acc + v.getTotalLength(), 0);
  const step = (totalLength / numPts);
  paths.forEach(path => {
    const sample = [];
    const n = M.floor(path.getTotalLength());
    for (let i = 0; i < n; i += step)
      sample.push(toPt(path.getPointAtLength(i)));
    samplePts = samplePts.concat(sample);
  });
  while (samplePts.length > numPts)
    samplePts.pop();
  return samplePts;
}

function toPt(pt: SVGPoint): [number, number] {
  return [pt.x, pt.y];
}

Discrete Fourier Transformatie (DFT)

Vervolgens passen we de DFT toe om de Fourier-coëfficiënten te vinden:

$$ ck=\frac{1}{N}\sum{j=0}^{N}{p(j)e^{-i\frac{2\pi}{N}kj}} $$

Waarbij $N$ het totaal aantal samplepunten is. De DFT zet $N$ gelijkmatig verdeelde complexe punten in het tijdsdomein om naar hun representatie in het frequentiedomein.

De implementatie in TypeScript ($\mathcal{O}(N^2)$):

import * as M from 'mathjs';
export type DFTData = M.Complex & {
  freq: number;
  radius: number;
  phase: number;
};
export async function dft(cps: M.Complex[]): Promise<DFTData[]> {
  const X = [];
  const N = cps.length;
  for (let i = 0; i < N; i++) {
    let sum = M.complex('0');
    for (let n = 0; n < N; n++) {
      const xn = cps[n];
      const theta = -(2 * M.pi * i * n) / N;
      const c = M.complex(`${M.cos(theta)}+${M.sin(theta)}i`);
      sum = M.add(sum, M.multiply(xn, c)) as M.Complex;
    }
    const { re, im } = M.multiply(1 / N, sum) as M.Complex;
    X.push({
      re,
      im,
      freq: i,
      radius: M.sqrt((re * re) + (im * im)),
      phase: M.atan2(im, re)
    });
  }
  return X;
}

Visualisatie

Nu we de coëfficiënten hebben, kunnen we de curve parametreren om deze te plotten:

function parameterize(cn: M.Complex[]) {
  return async (t: number): Promise<M.Complex> =>
    await sigma(k =>
      M.multiply(M.complex(cn[k - 1].re, cn[k - 1].im),
      M.exp(M.multiply(M.i, (k - 1) * t))), [1, numSamples]) as M.Complex;
}

Naarmate het aantal sinusvormen toeneemt, wordt de representatie van de lama nauwkeuriger. Het meeste werk wordt gedaan door de lage frequenties; de hogere frequenties vullen de kleinere details in, waardoor de benadering scherper wordt.

Om de animatie met de draaiende cirkels te maken, gebruiken we de p5.js bibliotheek:

function draw() {
  for (let i = 0; i < complex.length; i++) {
    const cn = complex[i];
    const prevx = x;
    const prevy = y;
    const { freq, radius, phase } = cn;
    x += radius * M.cos(freq * time + phase);
    y += radius * M.sin(freq * time + phase);
    stroke(230, 85);
    noFill();
    ellipse(prevx, prevy, radius * 2);
    stroke(230, 142);
    line(prevx, prevy, x, y);
  }
  points.push(createVector(x, y));
  stroke(230);
  noFill();
  beginShape();
  for (let i =0 ; i < points.length; i++)
    vertex(points[i].x, points[i].y);
  endShape();
  if (time > M.pi * 2)
    points.pop();
  time += ((2 * M.pi) / pts.length);
}