FFTWを使う

はじめに

FFTWの使い方を説明します。
DFT (discrete fourier transform) や DCT (discrete cosine transform) が計算できます。

download & make

FFTWをdownload / makeします。

console

mkdir fftw
cd fftw/
wget http://www.fftw.org/fftw-3.3.8.tar.gz
tar xzvf fftw-3.3.8.tar.gz 
cd fftw-3.3.8/
./configure
make
sudo make install

code example

FFTWを使ってDCT IVを計算する例を示します。

次の関数でDCTのplanを作成します。
　fftw_plan_r2r_1d

引数fftw_r2r_kindでDCT typeが指定できます。DCT IVは次を指定します。
　FFTW_REDFT11 (DCT IV)

DCT IV の式は次の通りです。

X_k = 2 \sum^{N-1}_{n=0}{x_n \cos \frac{(2n+1)(2k+1)}{4N}} \\

DCTの実行はfftw_executeで行います。

DCT IVを2度実行して復元できることを確認します。
ただ定義式そのままでは復元できないため計算後に2Nで割ります。

makefile

CFLAGS=-I. -Wall -Werror -O2
INCS=
OBJS=test.o
LIBS=-lfftw3 -lm
TARGET=test

%.o: %.c $(INCS)
	$(CC) $(CFLAGS) -c -o $@ $<

$(TARGET): $(OBJS)
	$(CC) $(CFLAGS) -o $@ $^ $(LIBS)

clean:
	rm -rf $(TARGET) *.o

test.c

# include <stdio.h>
# include <inttypes.h>
# include <fftw3.h>

inline uint64_t RDTSCP()
{
    uint32_t hi, lo, aux=0;
    __asm__ volatile("rdtscp" : "=a" (lo), "=d" (hi), "=c" (aux));
    return ((uint64_t)hi << 32) | lo;
}

uint64_t _rdtscp_start, _rdtscp_end;
inline void tic()
{
  _rdtscp_start = RDTSCP();
}

inline void toc(const char* str)
{
  _rdtscp_end = RDTSCP();
  printf("profile:%s %" PRIu64 " clock\n", str, _rdtscp_end-_rdtscp_start);
}

int main(int argc, char* argv[])
{
  int i, w=32, len=1024;
  fftw_plan p;
  static double in[1024], out[1024];
  
  for (i=0; i<1024; i++) {
    in[i] = i;
  }
  printf("input\n");
  for (i=0; i<len; i++) {
    printf("%+05.3e ", in[i]);
    if ((i+1)%32 == 0) {
      printf("\n");
    }
  }
  
  p = fftw_plan_r2r_1d(len/*int n*/, in, out,
    FFTW_REDFT11 /* fftw_r2r_kind */, FFTW_ESTIMATE /* flags */);
  tic();
  fftw_execute(p);
  toc("");
  
  printf("output forward\n");
  for (i=0; i<len; i++) {
    printf("%+05.3e ", out[i]);
    if ((i+1)%w == 0) {
      printf("\n");
    }
  }
  
  p = fftw_plan_r2r_1d(len/*int n*/, out, in,
    FFTW_REDFT11 /* fftw_r2r_kind */, FFTW_ESTIMATE /* flags */);
  tic();
  fftw_execute(p);
  toc("");
  
  printf("output inverse\n");
  for (i=0; i<len; i++) {
    printf("%+05.3e ", in[i]/2048.0);
    if ((i+1)%w == 0) {
      printf("\n");
    }
  }
  
  return 0;
}

実行結果

console

input
+0.000e+00 +1.000e+00 +2.000e+00 +3.000e+00 +4.000e+00 +5.000e+00 +6.000e+00 +7.000e+00 +8.000e+00 +9.000e+00 +1.000e+01 +1.100e+01 +1.200e+01 +1.300e+01 +1.400e+01 +1.500e+01 +1.600e+01 +1.700e+01 +1.800e+01 +1.900e+01 +2.000e+01 +2.100e+01 +2.200e+01 +2.300e+01 +2.400e+01 +2.500e+01 +2.600e+01 +2.700e+01 +2.800e+01 +2.900e+01 +3.000e+01 +3.100e+01 
...
+9.920e+02 +9.930e+02 +9.940e+02 +9.950e+02 +9.960e+02 +9.970e+02 +9.980e+02 +9.990e+02 +1.000e+03 +1.001e+03 +1.002e+03 +1.003e+03 +1.004e+03 +1.005e+03 +1.006e+03 +1.007e+03 +1.008e+03 +1.009e+03 +1.010e+03 +1.011e+03 +1.012e+03 +1.013e+03 +1.014e+03 +1.015e+03 +1.016e+03 +1.017e+03 +1.018e+03 +1.019e+03 +1.020e+03 +1.021e+03 +1.022e+03 +1.023e+03 
profile: 130306 clock
output forward
+4.845e+05 -5.393e+05 +2.329e+05 -2.080e+05 +1.378e+05 -1.283e+05 +9.762e+04 -9.274e+04 +7.556e+04 -7.259e+04 +6.162e+04 -5.963e+04 +5.202e+04 -5.059e+04 +4.501e+04 -4.393e+04 +3.966e+04 -3.883e+04 +3.545e+04 -3.478e+04 +3.205e+04 -3.150e+04 +2.924e+04 -2.878e+04 +2.689e+04 -2.650e+04 +2.488e+04 -2.455e+04 +2.316e+04 -2.287e+04 +2.166e+04 -2.140e+04 
...
+1.025e+03 -1.025e+03 +1.025e+03 -1.025e+03 +1.024e+03 -1.024e+03 +1.024e+03 -1.024e+03 +1.024e+03 -1.024e+03 +1.024e+03 -1.024e+03 +1.024e+03 -1.024e+03 +1.024e+03 -1.024e+03 +1.024e+03 -1.024e+03 +1.024e+03 -1.024e+03 +1.024e+03 -1.024e+03 +1.024e+03 -1.024e+03 +1.024e+03 -1.024e+03 +1.024e+03 -1.024e+03 +1.024e+03 -1.024e+03 +1.024e+03 -1.024e+03 
profile: 144551 clock
output inverse
+1.081e-13 +1.000e+00 +2.000e+00 +3.000e+00 +4.000e+00 +5.000e+00 +6.000e+00 +7.000e+00 +8.000e+00 +9.000e+00 +1.000e+01 +1.100e+01 +1.200e+01 +1.300e+01 +1.400e+01 +1.500e+01 +1.600e+01 +1.700e+01 +1.800e+01 +1.900e+01 +2.000e+01 +2.100e+01 +2.200e+01 +2.300e+01 +2.400e+01 +2.500e+01 +2.600e+01 +2.700e+01 +2.800e+01 +2.900e+01 +3.000e+01 +3.100e+01 
...
+9.920e+02 +9.930e+02 +9.940e+02 +9.950e+02 +9.960e+02 +9.970e+02 +9.980e+02 +9.990e+02 +1.000e+03 +1.001e+03 +1.002e+03 +1.003e+03 +1.004e+03 +1.005e+03 +1.006e+03 +1.007e+03 +1.008e+03 +1.009e+03 +1.010e+03 +1.011e+03 +1.012e+03 +1.013e+03 +1.014e+03 +1.015e+03 +1.016e+03 +1.017e+03 +1.018e+03 +1.019e+03 +1.020e+03 +1.021e+03 +1.022e+03 +1.023e+03