フーリエ級数の可視化プログラムの実装[Processing]

もはや2週間前のネタを真似ただけで新規性に乏しいのですが……(´・ω・`)
一応実装にあたって元ネタを書いておかないと申し訳ないのでツイートを貼らせていただきます……

これ( https://t.co/uMm0inKXeV )にインスパイアされて、円が１０個のバージョンを作ってみたらキモくなった pic.twitter.com/lUkBNNldy9

— どね (@donnay1224) 2016, 2月 5

k_1(x)=のところに好きな関数（数列）を入れて遊べるフーリエ級数視覚化マシーンを作りましたhttps://t.co/GmQo5NoZbz pic.twitter.com/vHrQ32FdWw

— 鯵坂もっちょ (@motcho_tw) 2016, 2月 6

という経緯です．

今回の実装はProcessingによるものなので上記ツイートのような数式で描くことができませんが，プログラミング言語というだけあって自由度は高いものとなっています．

円の数とか任意に決められますし．関数も自由に決められますし．

追記：鰺坂もっちょ先生がトゥギャってくれました→世はまさに大フーリエ時代 - Togetterまとめ

・ソースコード

例として矩形波を近似しています．関数funcの中身を[0,1)区間で変えてやればその形になりますが，無限大に発散させたりするとダメです．
あと数字が小さすぎる(1.0とか)だと小さすぎて見えませんので気をつけてください．

	// フーリエ級数可視化クラス
	class FourierPlotter {
	int n_data; // データの個数
	int n_approx; // 近似次数
	float[][] radius; // 半径．radius[0][k]にsin成分のk次の成分が入る
	float[] data; // データ
	float centerX; // 中心X座標
	float centerY; // 中心Y座標
	float phase; // 位相のズレ
	public float x; // 先端部分のX座標
	public float y; // 先端部分のY座標
	// コンストラクタ
	FourierPlotter(float[] data, int n_approx,
	float centerX, float centerY, float phase) {
	// 引数をメンバ変数に代入した後，フーリへ級数展開
	this.data = data;
	this.n_data = data.length;
	this.n_approx = n_approx;
	this.phase = phase;
	this.centerX = this.x = centerX;
	this.centerY = this.y = centerY;
	this.calcRadius();
	}
	// フーリエ級数展開を計算する関数
	void calcRadius() {
	// フーリエ級数の係数は半径に当たる
	this.radius = new float[2][this.n_approx];
	float dt = 2.0*PI / this.n_data; // 微小長さdt
	// バイアス項(0次の項)は無視している
	for ( int i = 1; i < this.n_approx; i++ ) {
	for (int j = 0; j < this.n_data; j++ ) {
	float t = 2.0*PI*j/this.n_data;
	// 普通に矩形積分する
	this.radius[0][i] += sin(i*t)*this.data[j]*dt/PI; // sin成分の計算
	this.radius[1][i] += cos(i*t)*this.data[j]*dt/PI; // cos成分の計算
	}
	}
	}
	// 描画を行う関数
	void visualize(float theta) {
	// x,y座標をとりあえず中心に持ってくる
	this.x = this.centerX;
	this.y = this.centerY;
	float tempX, tempY;
	for (int i = 0; i < this.n_approx; i++ ) {
	for (int j= 0; j<2; j++) {
	// 円を描く
	ellipse(x, y, 2.0*this.radius[j][i], 2.0*this.radius[j][i]);
	tempX = x;
	tempY = y;
	// 先端部分の位置座標を計算(j=1の時，余弦なので位相がPI/2ずれる)
	x += this.radius[j][i]*cos(i*theta+PI*j/2 + this.phase);
	y += this.radius[j][i]*sin(i*theta+PI*j/2 + this.phase);
	line(x, y, tempX, tempY); // 各円の中心をを直線で結ぶ
	}
	}
	}
	}
	float scale = 1.0; // 拡大率
	float dx=0.0; // x方向の移動
	float dy=-0.0; // y方向の移動
	int n_data = 100; // データの個数
	FourierPlotter fp; // 描画クラス
	float[] points; // 過去に通った座標
	int count; // ループカウンタ
	int n_approx = 20; // 近似次数
	void setup() {
	size(640, 480);
	// 関数funcを離散的なデータにしてFourierPlotterに渡す
	float[] data = new float[n_data];
	for (int i = 0; i < n_data; i++ )
	data[i] = func( (float)(i)/n_data );
	fp = new FourierPlotter(data, n_approx, -100, 0, 0);
	// 諸々を初期化
	count = 0;
	points = new float[1000];
	}
	void draw() {
	background(0);
	translate(width/2-dx, height/2-dy); // 原点を中心に移動
	scale(scale, scale);
	fill(0, 0, 0, 0);
	strokeWeight(3);
	stroke( 0, 255, 0);
	float theta = count * 2*PI / n_data;
	fp.visualize(theta); // 描画
	points[count%points.length] = fp.y; // y座標を保存
	// 過去のデータの描画
	stroke(255,0,0);
	line(fp.x,fp.y,100.0,fp.y);
	plot(points,100,0,count%points.length);
	count +=1; // ループカウンタを増やす
	// フレームを保存する場合はコメントインする
	//saveFrame("./frames/#####.png");
	}
	// [0,1)区間で定義される区分的に滑らかな関数
	// 近似する関数の対象なのでここを変えればいい
	float func(float x) {
	// 矩形波
	return x < 0.5 ? -50:50;
	// サイン波
	//return 100*sin(2*x*PI);
	}
	// 配列arrをプロットする関数
	void plot(float[] arr,float x,float y,int startIndex){
	int N = arr.length;
	float y2 = arr[startIndex%N];
	for(int i = 1 ; i < N ; i++ ){
	line(x+N-i+1,y2,x+N-i,y+arr[(startIndex+i)%N]);
	y2 = arr[(startIndex+i)%N] + y;
	}
	}
	void keyPressed() {
	// キー操作による拡大・移動
	if (key == CODED) {
	if (keyCode == UP) {
	dy -= 10;
	} else if (keyCode == DOWN) {
	dy += 10;
	} else if (keyCode == LEFT) {
	dx -= 10;
	} else if (keyCode == RIGHT) {
	dx += 10;
	}
	} else if (key == '+') {
	scale *= 1.05;
	} else if (key=='-') {
	scale /= 1.05;
	}
	}

view raw VisualizingFourierSeries.pde hosted with ❤ by GitHub