フーリエ級数
よくある円でフーリエ級数を表現するアレ
デモ
クリックで再描画
原理
フーリエ級数とは、周期関数を三角関数の無限和で表したもの。
\[ f(t) = \frac {a_0} 2 + \sum_{k=1}^{\infty} {a_k \cos{2\pi fkt} + b_k \sin{2\pi fkt}} \]矩形波と三角波はそれぞれ以下のようにフーリエ級数展開できる(奇関数なので定数項とcosの係数は0)。
矩形波
\begin{eqnarray} f(t) &=& \frac 4 {\pi} \sum_{n=1}^{\infty} \frac {\sin{\{(2n-1)2\pi ft}\}} {(2n-1)} \\ &=& \frac 4 \pi \left( \frac {\sin {2\pi ft}} 1 + \frac {\sin {\left(3 \times 2\pi ft \right)}} 3 + \frac {\sin {\left( 5 \times 2\pi ft \right)}} 5 + \cdots \right) \end{eqnarray}三角波
\begin{eqnarray} f(t) &=& \frac 8 {\pi^2} \sum_{n=1}^{\infty} \sin { \left\{ \frac{(2n-1)\pi} 2 \right\} } \frac{\sin{\{(2n-1)2\pi ft\}}}{(2n-1)^2} \\ &=& \frac 8 {\pi^2} \left( \frac {\sin {2\pi ft}} 1 - \frac {\sin {\left(3 \times 2\pi ft \right)}} 9 + \frac {\sin {\left( 5 \times 2\pi ft \right)}} {25} + \cdots \right) \end{eqnarray}さて、中心\( C\left(c_x, c_y\right)\)、半径\(r\)の円弧上を、点\(P\left(p_x(t), p_y(t)\right)\)が反時計回りに周期\(\frac 1 f\)で等速円運動しているとすると、
時刻tにおけるPの座標は
\[ (c_x + r\cos(2\pi ft), c_y + r\sin(2\pi ft)) \cdots (1) \]時計回りの場合は
\[ (c_x + r\cos(2\pi ft), c_y - r\sin(2\pi ft)) \]である。ここで円\(k\): 中心\(\left(c_{kx}, c_{ky}\right)\), 半径\(r_k (k = 0, 1, 2, 3, ...) \)を考え、 円\(k\)の中心を、円\(k-1\)上を周期\(f_{k-1}\)で等速円運動する点\(P_{k-1}\)と同じ座標になるように定義する。
また、\(r_0\)と\(r_k\)の比を\(g(k)\), \(f_0\)と\(f_k\)の比を\(h(k)\)で表す。
さらに回転方向を表すために、kに対して1または-1を返す関数\(s(k)\)を用意する。
\begin{eqnarray}\\ c_{kx}(t) &=& \begin{cases} 0 & (k=0) \\ p_{(k-1)x}(t) & (k = 1,2,3, ...) \\ \end{cases} \\ c_{ky}(t) &=& \begin{cases} 0 & (k=0) \\ p_{(k-1)y}(t) & (k = 1,2,3, ...) \\ \end{cases} \\ r_k &=& \begin{cases} r_0 & (k=0) \\ \frac {r_0} {g(k)} & (k = 1,2,3, ...) \end{cases} \\ f_k &=& \begin{cases} 1 & (k=0) \\ f_0 h(k) & (k = 1,2,3, ...) \end{cases} \\ \end{eqnarray}このとき(1)を使うと、kについての漸化式が得られる。
\begin{eqnarray} c_{kx}(t) &=& \begin{cases} 0 & (k=0) \\ c_{(k-1)x}(t) + \frac {r_0} {g(k)} \cos { \frac{2\pi f_0 t} {h(k)} } & (k = 1,2,3, ...) \\ \end{cases} \\ c_{ky}(t) &=& \begin{cases} 0 & (k=0) \\ c_{(k-1)y}(t) + s(k) \frac {r_0} {g(k)} \sin { \frac{2\pi f_0 t} {h(k)} } & (k = 1,2,3, ...) \\ \end{cases} \\ \end{eqnarray}円nの円周上の点の座標は
\begin{eqnarray} p_{nx}(t) &=& c_{(n+1)x}(t) \\ &=& r_0\sum_{k=1}^{n+1} {\frac 1 {g(k)} \cos { 2\pi f_0 h(k) t }} \\ p_{ny}(t) &=& c_{(n+1)y}(t) \\ &=& r_0\sum_{k=1}^{n+1} {s(k)\frac 1 {g(k)} \sin { 2\pi f_0 h(k) t} } \end{eqnarray}\(n\to\infty\)の極限をとると(本当は収束することの証明が必要)
\begin{eqnarray} p_{\infty x}(t) &=& r_0\sum_{k=1}^{\infty} {\frac 1 {g(k)} \cos { 2\pi f_0 h(k) t }} \\ p_{\infty y}(t) &=& r_0\sum_{k=1}^{\infty} {s(k)\frac 1 {g(k)} \sin { 2\pi f_0 h(k) t} } \end{eqnarray}この式をよく見ると、\(p_{\infty x}(t)\)は偶関数のフーリエ級数展開、\(p_{\infty y}(t)\)は奇関数のフーリエ級数展開になっており
\begin{eqnarray} r_0 &=& \frac 4 \pi \\ s(k) &=& 1 \\ g(k) &=& 2k - 1 \\ h(k) &=& 2k -1 \end{eqnarray}とすれば\(p_{ny}(t)\)は矩形波、
\begin{eqnarray} r_0 &=& \frac 8 {\pi^2} \\ s(k) &=& \sin {\frac{(2n-1)\pi} {2}} \\ g(k) &=& 2k - 1 \\ h(k) &=& (2k - 1)^2 \end{eqnarray}とすれば\(p_{ny}(t)\)は三角波になる。
実装
まず一般化しておく
abstract class FourierCircle {
final float R0;
float cx, cy, r, n, px, py;
int time;
public FourierCircle(float cx, float cy, float n ,int time, float r) {
this.cx = cx;
this.cy = cy;
this.n = n;
this.R0 = r;
this.r = r * r0() / radiusRatio(n);
this.time = time;
px = this.cx + this.r * cos(2 * PI * freqRatio(n) * time / PERIOD);
py = this.cy + sign(n) * this.r * sin(2 * PI * freqRatio(n) * time / PERIOD);
}
/* abstract methods */
public abstract float r0();
public abstract float sign(float n);
public abstract float freqRatio(float n);
public abstract float radiusRatio(float n);
public abstract FourierCircle getChild();
}
FourierCircleを継承して矩形波、三角波を作成
class SquareWave extends FourierCircle {
/* square wave */
public SquareWave(float cx, float cy, float n ,int time, float r) {
super(cx, cy, n , time, r);
}
public float r0() {
return 4 / PI;
}
public float sign(float n) {
return 1;
}
public float freqRatio(float n) {
return 2 * n - 1;
}
public float radiusRatio(float n) {
return 2 * n - 1;
}
public SquareWave getChild() {
return new SquareWave(px, py, n + 1 ,this.time, R0);
}
}
class TriangleWave extends FourierCircle {
/* triangle wave */
public TriangleWave(float cx, float cy, float n ,int time, float r) {
super(cx, cy, n , time, r);
}
public float r0() {
return 8 / PI / PI;
}
public float sign(float n) {
return n % 2 == 1 ? 1: -1;
}
public float freqRatio(float n) {
return 2 * n - 1;
}
public float radiusRatio(float n) {
return (2 * n - 1) * (2 * n - 1);
}
public TriangleWave getChild() {
return new TriangleWave(px, py, n + 1 ,this.time, R0);
}
}