フーリエ級数を可視化するプログラムです．例として矩形波の近似を行っています．

VisualizingFourierSeries.pde

// フーリエ級数可視化クラス
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;
  }
}