トップ | puarts.com
ロゴ
「フーリエ変換」に関連する記事一覧
0  

FFT についてとてもわかりやすい資料があったのでメモ代わりに貼っておきます。

http://lab.sdm.keio.ac.jp/maenolab/makino/documents/20090827_Fourier.pdf

一度理解しても、数年後にまた使う時に細かいところを忘れているので、こういう資料があるとすぐに思い出せて助かります。

Reliable Software Frequency Analyzer
http://www.relisoft.com/freeware/index.htm

音声にFFTを手軽に適用してみることができるソフトを見つけましたのでメモしておきます。しかも、ソースコードつきです。

ファイル入力もマイク入力もどちらも出来ていい感じです。

画像に試しにFFTを適用してみたいという時に手軽に使えるフリーソフトないかなぁと思っていて、MeeSoft Image Analyzerというソフトを見つけました。

使った感じだと大分バグが多いのが嫌ですが、とりあえず自分の作ったFFTのプラグラムとちゃんと同じ絵が出るか試すのに使うには手軽でいいです。

MeeSoft Image Analyzer
http://meesoft.logicnet.dk/

C++でFFTを行う上で複素数を扱うデータのデータ構造をどうするのが良いのか調べるために2種類のデータ構造でFFTの処理時間を検証してみました。 まずひとつ目はFFTWで使われている構造体と同じように実数と虚数が交互に入っている配列型のデータ構造。こちらが一般的なんでしょうか。こちらはデータ構造Aとします。

もうひとつは実数と虚数のデータ配列を分けて2つの配列を使うデータ構造です。こちらをデータ構造Bとします。

配列が2つに分離するとメモリ確保のコストと実数と虚数に同時にメモリアクセスする際のストライドが増えてしまいますが、実数と虚数を別々にオペレーションしなければならないときはこのように分かれていると高速に処理できるので、こちらが有効な場合も考えられます。

今回はFFTにおける処理時間をこの2つのデータ構造により、どの程度変化が出るのか検証しました。
複素数の配列の要素数は2^23個で試行回数は30回で平均をとりました。以下が検証結果です。

環境
OS:  Windows 7 64bit
CPU:  core-i7 920
Memory:  6GB

データ構造A:  3.255 sec
データ構造B:  5.305 sec

相互に実数と虚数が格納されている配列型のデータ構造Aの方が2秒ほど速い結果になりました。実数と虚数を別々に処理することがない限りはメモリアクセス時のストライドはAの方が小さいので当然と言えば当然なのですが、FFTを実装する上では複素数のデータはデータ構造Aのようにして作った方が明確になったので、こちらで作った方が良さそうです。

使用したコードも記しておきます。

実数と虚数を交互に格納したデータ構造を用いたコード


int FFT1d(
        double *io_complex,        
        int size, int exponent,        
        bool ifft)
{
        int i1, i2,l,l1,l2;
        double c1,c2,tx,ty,t1,t2,u1,u2,z;
        int complex_size = size * 2;
                
        //** argument check -------------------- 
        if(size < 2){
                WSPERROR("size is less than 2\n");        
                return -1;
        }
        if(exponent<=0)
        {
                int i, j;
                for(exponent=0, i=size; (i=i/2)!=0; ++exponent){}
                for(i=0, j=1; i<exponent; ++i, j*=2){}
                if(size != j){
                        WSPERROR("size is not power of 2, %d != %d\n", size, j);        
                        return -1;
                }
        }
        //** ------------------------------------ 
        //printf("size: %d, exponent: %d\n", size, exponent);
        
        //** reverse bit
        i2 = size >> 1;
        {
                double *re_ptr, *im_ptr;
                double *re_ptr_j, *im_ptr_j, *end_ptr;
                int i, j, k;
                re_ptr = io_complex;
                end_ptr = re_ptr + complex_size - 2;
                for (j=0; re_ptr!=end_ptr; re_ptr+=2) {
                        re_ptr_j = io_complex + j*2;
                        if (re_ptr < re_ptr_j){
                                im_ptr = re_ptr + 1;
                                im_ptr_j = re_ptr_j + 1;
                                //** swap complex[i] and complex[j]
                                tx = *re_ptr;
                                ty = *im_ptr;
                                *re_ptr = *re_ptr_j;
                                *im_ptr = *im_ptr_j;
                                *re_ptr_j = tx;
                                *im_ptr_j = ty;
                        }
                                
                        for(k = i2; k<=j; k>>=1){ j-=k; }
                        j+=k;
                }
        }
        
        //** Compute FFT
        {
                double *real_ptr, *imag_ptr, *end_ptr;
                double *real_ptr_2, *imag_ptr_2,        
                double *end_ptr_2, *real_ptr_i1, *imag_ptr_i1;
                int i, j, k, l2_step, l1_step;
                c1 = -1.0;        
                c2 = 0.0;
                l2 = 1;
                for (l=0; l<exponent; ++l) {
                        l1 = l2;
                        l2 <<= 1;
                        u1 = 1.0;        
                        u2 = 0.0;
                        l2_step = l2 * 2;
                        l1_step = l1 * 2;
                        real_ptr = io_complex;
                        end_ptr = real_ptr + l1_step;
                        end_ptr_2 = io_complex + complex_size;
                        for(; real_ptr!=end_ptr; real_ptr+=2) {
                                real_ptr_2 = real_ptr;
                                for(; real_ptr_2<end_ptr_2; real_ptr_2+=l2_step){
                                        imag_ptr_2 = real_ptr_2 + 1;
                                        real_ptr_i1 = real_ptr_2 + l1_step;
                                        imag_ptr_i1 = real_ptr_i1 + 1;
                                        t1 = u1 * (*real_ptr_i1) - u2 * (*imag_ptr_i1);
                                        t2 = u1 * (*imag_ptr_i1) + u2 * (*real_ptr_i1);
                                        *real_ptr_i1 = *real_ptr_2 - t1;        
                                        *imag_ptr_i1 = *imag_ptr_2 - t2;
                                        *real_ptr_2 += t1;
                                        *imag_ptr_2 += t2;
                                }
                                z =         u1 * c1 - u2 * c2;
                                u2 = u1 * c2 + u2 * c1;
                                u1 = z;
                        }
        
                        c2 = sqrt((1.0 - c1) / 2.0);
                        if (ifft==false){ c2 = -c2; }
                        c1 = sqrt((1.0 + c1) / 2.0);
                }
        }
        //** Scale array size for forward fft
        if(ifft==false) {
                double *ptr, *end_ptr;
                ptr = io_complex;
                end_ptr = ptr + complex_size;
                for (; ptr!=end_ptr; ptr+=2) {
                        *ptr /= size;
                        *(ptr+1) /= size;
                }
        }
        return 0;
}

実数と虚数を別の配列で分けたコード


int FFT1d(double *io_real,
    double *io_imag,
    int size, int m, bool ifft)
{
        int i1, i2,l,l1,l2;
        double c1,c2,tx,ty,t1,t2,u1,u2,z;
                
        //** argument check --------------------------------- 
        if(size < 2){
                WSPERROR("size is less than 2\n");
                return -1;
        }
        if(m<=0)
        {
                int i, j;
                for(m=0, i=size; (i=i/2)!=0; ++m){}
                for(i=0, j=1; i<m; ++i, j*=2){}
                if(size != j){
                        WSPERROR("size is not power of 2, %d != %d\n", size, j);        
                        return -1;
                }
        }
        //** -------------------------------- 
        //printf("size: %d, m: %d\n", size, m);
        
        //** reverse bit
        i2 = size >> 1;
        {
                double *real_ptr, *imag_ptr, *end_ptr;
                double *real_ptr_j, *imag_ptr_j;
                int i, j, k;
                real_ptr = io_real;
                imag_ptr = io_imag;
                end_ptr = io_real + size - 1;
                for (j=0; real_ptr!=end_ptr; ++real_ptr, ++imag_ptr) {
                        real_ptr_j = io_real + j;
                        if (real_ptr < real_ptr_j){
                                imag_ptr_j = io_imag + j;
                                //** swap complex[i] and complex[j]
                                tx = *real_ptr;
                                ty = *imag_ptr;
                                *real_ptr = *real_ptr_j;
                                *imag_ptr = *imag_ptr_j;
                                *real_ptr_j = tx;
                                *imag_ptr_j = ty;
                        }
                                
                        for(k = i2; k<=j; k>>=1){ j-=k; }
                        j+=k;
                }
        }
        
        //** Compute FFT
        {
                double *real_ptr, *imag_ptr, *end_ptr;
                double *real_ptr_2, *imag_ptr_2;
                double *end_ptr_2, *real_ptr_i1, *imag_ptr_i1;
                int i, j, k;
                c1 = -1.0;        
                c2 = 0.0;
                l2 = 1;
                for (l=0; l<m; ++l) {
                        l1 = l2;
                        l2 <<= 1;
                        u1 = 1.0;        
                        u2 = 0.0;
                        real_ptr = io_real;
                        imag_ptr = io_imag;
                        end_ptr = io_real + l1;
                        end_ptr_2 = io_real + size;
                        for(; real_ptr!=end_ptr; ++real_ptr, ++imag_ptr) {
                                real_ptr_2 = real_ptr;
                                imag_ptr_2 = imag_ptr;
                                for(; real_ptr_2<end_ptr_2; real_ptr_2+=l2, imag_ptr_2+=l2){
                                        real_ptr_i1 = real_ptr_2 + l1;
                                        imag_ptr_i1 = imag_ptr_2 + l1;
                                        t1 = u1 * (*real_ptr_i1) - u2 * (*imag_ptr_i1);
                                        t2 = u1 * (*imag_ptr_i1) + u2 * (*real_ptr_i1);
                                        *real_ptr_i1 = *real_ptr_2 - t1;        
                                        *imag_ptr_i1 = *imag_ptr_2 - t2;
                                        *real_ptr_2 += t1;
                                        *imag_ptr_2 += t2;
                                }
                                z =         u1 * c1 - u2 * c2;
                                u2 = u1 * c2 + u2 * c1;
                                u1 = z;
                        }
        
                        c2 = sqrt((1.0 - c1) / 2.0);
                        if (ifft==false){ c2 = -c2; }
                        c1 = sqrt((1.0 + c1) / 2.0);
                }
        }
        //** Scale size for forward FFT
        if(ifft==false) {
                double *real_ptr, *imag_ptr, *end_ptr;
                real_ptr = io_real;
                imag_ptr = io_imag;
                end_ptr = real_ptr + size;
                for (; real_ptr!=end_ptr; ++real_ptr, ++imag_ptr) {
                        *real_ptr /= size;
                        *imag_ptr /= size;
                }
        }
        return 0;
}

OpenCV2.2のcv::gpu::dftを使って4x4の複素行列をフーリエ変換、逆フーリエ変換をする例です。なかなかネット上にサンプルコードがなくて少し使い方に苦労したので、同じようにOpenCVのdftをGPU機能つきで使用したい方の参考になればと思います。

This is a sample code to apply fourier transform and inverse fourier transform with cv::gpu::dft function in OpenCV 2.2.

Source code

#include<opencv2/gpu/gpu.hpp>

void cvGpuTest_DFT_IDFT(){
 int w=4, h=4;
 int len=w*h*2;
 cv::Mat dst;
 cv::gpu::GpuMat gpuMat;

 //** make source data -----
 float *srcData;
 srcData=new float[len];
 for(int i=0; i<len/2; i++){
  srcData[i*2]=i;
  srcData[i*2+1]=0;
 }
 cv::Mat srcMat(w, h, CV_32FC2, srcData);

 //** Print Source Data ----
 printf("SRC------------\n");
 for(int y=0; y<h; y++){
  for(int x=0; x<w; x++){
   printf("%d %d: %f %f\n", x, y, srcMat.at<float>(y, x*2), srcMat.at<float>(y, x*2+1));
  }
 }

 //** DFT ------------------
 gpuMat = srcMat;

 cv::gpu::dft(gpuMat, gpuMat, cv::Size(w, h));

 dst=gpuMat;

 //** Print DFT Result ------
 printf("DFT------------\n");
 for(int y=0; y<h; y++){
  for(int x=0; x<w; x++){
   printf("%d %d: %f %f\n", x, y, dst.at<float>(y,x*2), dst.at<float>(y, x*2+1));
  }
 }

 //** IDFT -----------------
 gpuMat=dst;

 cv::gpu::dft(gpuMat, gpuMat, cv::Size(w, h),

              cv::DFT_INVERSE);

 dst=gpuMat;

 //** Print IDFT Result -----
 printf("IDFT------------\n");
 for(int y=0; y<h; y++){
  for(int x=0; x<w; x++){
   printf("%d %d: %f %f\n", x, y, dst.at<float>(y,x*2)/(w*h), dst.at<float>(y, x*2+1));
  }
 }
}

int main()
{
 cvGpuTest_DFT_IDFT();
 return 0;
}


Result

SRC------------ 

0 0: 0.000000 0.000000

1 0: 1.000000 0.000000

2 0: 2.000000 0.000000

3 0: 3.000000 0.000000

0 1: 4.000000 0.000000

1 1: 5.000000 0.000000

2 1: 6.000000 0.000000

3 1: 7.000000 0.000000

0 2: 8.000000 0.000000

1 2: 9.000000 0.000000

2 2: 10.000000 0.000000

3 2: 11.000000 0.000000

0 3: 12.000000 0.000000

1 3: 13.000000 0.000000

2 3: 14.000000 0.000000

3 3: 15.000000 0.000000

DFT------------ 

0 0: 120.000000 0.000000

1 0: -8.000000 8.000000

2 0: -8.000000 0.000000

3 0: -8.000000 -8.000000

0 1: -32.000000 32.000000

1 1: 0.000000 0.000000

2 1: 0.000000 0.000000

3 1: 0.000000 0.000000

0 2: -32.000000 0.000000

1 2: 0.000000 0.000000

2 2: 0.000000 0.000000

3 2: 0.000000 0.000000

0 3: -32.000000 -32.000000

1 3: 0.000000 0.000000

2 3: 0.000000 0.000000

3 3: 0.000000 0.000000

IDFT------------ 

0 0: 0.000000 0.000000

1 0: 1.000000 -0.000000

2 0: 2.000000 0.000000

3 0: 3.000000 0.000000

0 1: 4.000000 0.000000

1 1: 5.000000 -0.000000

2 1: 6.000000 0.000000

3 1: 7.000000 0.000000

0 2: 8.000000 0.000000

1 2: 9.000000 -0.000000

2 2: 10.000000 0.000000

3 2: 11.000000 0.000000

0 3: 12.000000 0.000000

1 3: 13.000000 -0.000000

2 3: 14.000000 0.000000

3 3: 15.000000 0.000000

C言語でフーリエ変換した画像を周波数成分の振幅値で除算して位相部分のみ抜き出して画像に出力するプログラムを作りました。FFTにはFFTW、画像の入出力にはOpenCVを使っています。
下の画像が入力と出力です。

I made a C program for extration of phase channel from gray scale image. In this progaram, I used FFTW for FFT, and OpenCV for image I/O.

Here is an input image.
inputImg.jpg

Here is a result image.
4a375fbe.jpeg

ソースコードも置いておきます。

Here is source code
phaseExtraction.c

フーリエ変換の原理はだいたいわかるんですが、2次元フーリエ変換の場合、元画像がどう変化すると周波数空間はどのように変化するかということは直感的につかみ難いです。

ということで、フーリエ変換で周波数空間と元画像にどのような関係があるのかを自作プログラムで動画をフーリエ変換して実験してみました。下のような6角形が回転した画像を変換すると周波数空間の画像もそれに合わせて回転しました。

周波数空間と元画像の関係のgifアニメ

非常におもしろい性質です。

ご覧の通り、元画像を周波数空間の画像に重ねると、周波数空間の画像が元画像の法線のようにも見えます。

何かこういった性質を利用したアートができないものかと考え中です。とりあえず、もう少しいろいろと実験は必要です。

C言語で DFT のプログラムを書いてみました。

0  

0.0315 sec
にほんブログ村 ゲームブログ ファイアーエムブレムへ にほんブログ村 デザインブログ コンピュータグラフィックスへ

Copyright(C)2006-2018 wsp All Rights Reserved