FFT で三角関数(sin)を１回しか呼ばない方法

FFT では sin/cos 関数を何度もよびださないように、予めテーブルを作ったりしますが、テーブル生成時の三角関数を一回で済ませます。

sin テーブル生成の例

C++

void make_table1(double *table, unsigned int size)
{
  for (unsigned int i = 0; i < size; i++)
    table[i] = sin(M_PI*2 * i/size);
}

を代替します。座標$\left(x,y\right)$ を $\left(1,0\right)$ から $\left(cos45°,sin45°\right)$ まで回転するときの座標が $\left(cos,sin\right)$ なので、それを使います。計算量は減りますが、精度も落ちます。

C++

void make_table2(double *table, unsigned int size)
{
  assert((size >= 4) && (size == (size & -size)));

  unsigned int h_sz = size >> 1;
  unsigned int q_sz = h_sz >> 1;
  unsigned int o_sz = q_sz >> 1;
  unsigned int qhsz = h_sz + q_sz;

  double ds = sin(M_PI / h_sz); // sin 関数は1回だけ
  double dc = sqrt(1 - ds * ds);
  double px = 1, py = 0;

  table[0    + 0   ] =  0;
  table[0    + q_sz] = +1;
  table[h_sz + 0   ] =  0;
  table[h_sz + q_sz] = -1;
  for (unsigned int i = 1; i < o_sz; i++)
    {
      double nx = px * dc - py * ds;
      double ny = px * ds + py * dc;

      if (false)
        { // 少しだけ演算精度がよくなる
          double r = sqrt(nx * nx + ny * ny);
          nx /= r;
          ny /= r;
        }

      table[0    + i] = +ny;
      table[q_sz - i] = +nx;
      table[q_sz + i] = +nx;
      table[h_sz - i] = +ny;
      table[h_sz + i] = -ny;
      table[qhsz - i] = -nx;
      table[qhsz + i] = -nx;
      table[size - i] = -ny;

      px = nx;
      py = ny;
    }

  double ir2 = sqrt(2.0) / 2;

  table[0    + o_sz] = +ir2;
  table[q_sz + o_sz] = +ir2;
  table[h_sz + o_sz] = -ir2;
  table[qhsz + o_sz] = -ir2;
}

テーブル生成テスト プログラム

sample.cpp

# include <cassert>
# include <cstdio>
# include <cmath>

using namespace std;

// 全て sin 関数に頼る
void make_table1(double *table, unsigned int size)
{
  for (unsigned int i = 0; i < size; i++)
    table[i] = sin(M_PI*2 * i/size);
}

void make_table2(double *table, unsigned int size)
{
  assert((size >= 4) && (size == (size & -size)));

  unsigned int h_sz = size >> 1;
  unsigned int q_sz = h_sz >> 1;
  unsigned int o_sz = q_sz >> 1;
  unsigned int qhsz = h_sz + q_sz;

  double ds = sin(M_PI / h_sz); // sin 関数は1回だけ
  double dc = sqrt(1 - ds * ds);
  double px = 1, py = 0;

  table[0    + 0   ] =  0;
  table[0    + q_sz] = +1;
  table[h_sz + 0   ] =  0;
  table[h_sz + q_sz] = -1;
  for (unsigned int i = 1; i < o_sz; i++)
    {
      double nx = px * dc - py * ds;
      double ny = px * ds + py * dc;

      if (false)
        { // 少しだけ演算精度がよくなる
          double r = sqrt(nx * nx + ny * ny);
          nx /= r;
          ny /= r;
        }

      table[0    + i] = +ny;
      table[q_sz - i] = +nx;
      table[q_sz + i] = +nx;
      table[h_sz - i] = +ny;
      table[h_sz + i] = -ny;
      table[qhsz - i] = -nx;
      table[qhsz + i] = -nx;
      table[size - i] = -ny;

      px = nx;
      py = ny;
    }

  double ir2 = sqrt(2.0) / 2;

  table[0    + o_sz] = +ir2;
  table[q_sz + o_sz] = +ir2;
  table[h_sz + o_sz] = -ir2;
  table[qhsz + o_sz] = -ir2;
}

void compare(const double *l, const double *r, unsigned int n)
{
  double dmin = +1;
  double dmax = +0;

  printf("--------------------------------\n");
  for (unsigned int i = 0; i < n; i++)
    {
      double d = abs(l[i] - r[i]);

      printf("%4d/%4d: %+.8e, %+.8e, %+.8e\n", i, n, l[i], r[i], d);

      if (d > dmax) dmax = d;
      if (d < dmin) dmin = d;
    }
  printf("dmin=%e, dmax=%e\n", dmin, dmax);
}

int main(int argc, char **argv)
{
  unsigned int size = (1 << 10);
  double table1[size];
  double table2[size];
  make_table1(table1, size);
  make_table2(table2, size);
  compare(table1, table2, size);
  return 0;
}

誤差は
if(false) の場合： dmin=0.000000e+00, dmax=4.218847e-15
if(true) の場合： dmin=0.000000e+00, dmax=7.771561e-16

---

配列が $2^n$ ならば、最大の $n$ を決めて sin 値もテーブル化すると、sin 関数の呼び出しはゼロになりますが...