• 締切済み

Cプログラムによる画像の高速フーリエ変換FFT

Cプログラムにる画像処理について教えてください。逆高速フーリエ変換(IFFT)によって周波数領域から実空間領域に変換する際、周波数帯域を(低周波数領域に)制限してIFFTを実行できるという文献があるのですが可能でしょうか?通常の逆離散フーリエ変換(IDFT)では可能でした。 以下のコードの修正で可能でしょうか? 例) 512*512の空間周波数領域データを212*212(仮に)の低周波数領域に帯域制限してIFFTで画像に戻したいのです。 以下の以下の一次元FFTコードを用いてX方向とY方向にIFFTしたいと考えています。 void FFT(int ir, int nx, float *xr, float *xi, float *si, float *co, unsigned short *brv) // 1次元フーリエ変換 // int ir; 順変換(1)と逆変換(-1) // int nx; 1次元FFTのデータ数 // float *xr; 実部のデータ xr[nx] // float *xi; 虚部のデータ xi[nx] // float *si; FFT用のサインデータ si[nx/2] // float *co; FFT用のコサインデータ co[nx/2] // unsigned short *brv; FFT用の入れ替えデータ brv[nx] { int i, j, n1, n2=nx, j3, j4, k, l, ll, d=1, g; float a, b, c, s; for(l = 1; l <= nx/2; l *= 2, d += d) { g = 0; ll = n2; n2 /= 2; for(k = 1; k <= n2; k++) { n1 = k-ll; c = co[g]; s = -ir*si[g]; g += d; for(j = ll; j <= nx; j += ll) { j3 = j+n1-1; j4 = j3+n2; a = xr[j3]-xr[j4]; b = xi[j3]-xi[j4]; xr[j3] += xr[j4]; xi[j3] += xi[j4]; xr[j4] = c*a+s*b; xi[j4] = c*b-s*a; } } } 通常のDFTでの帯域制限はこのように行い可能でした。 void fourier1d(int ir, float *fr, float *fi, int nx) { int i, j, n = 1; float *gr, *gi; double u, x; gr = (float *)malloc((unsigned long)nx*sizeof(float)); gi = (float *)malloc((unsigned long)nx*sizeof(float)); for(i = 0 ; i < nx ; i++) { u = i-nx/2; gr[i] = gi[i] = 0; for(j = 150 ; j < 362 ; j++) { x = j-nx/2; gr[i] += (float)( fr[j]*cos(2*PI*u*x/nx)+ir*fi[j]*sin(2*PI*u*x/nx)); gi[i] += (float)(-ir*fr[j]*sin(2*PI*u*x/nx)+fi[j]*cos(2*PI*u*x/nx)); } } if(ir == -1) n = 212; // 逆変換はデータ数で割る for(i = 0 ; i < nx ; i++) { fr[i] = gr[i]/n; fi[i] = gi[i]/n; } free(gr); free(gi); }

みんなの回答

  • Tacosan
  • ベストアンサー率23% (3656/15482)
回答No.2

文章は逆フーリエ変換だけどコードは (逆じゃない) フーリエ変換だ. どっちが正しいんだろう. いずれにしても, 「FFT ないし IFFT の途中で帯域制限」はできないことはないだろうけどする必要があるかどうかが疑問. #1 でも言われてるけどおそらくかなり面倒なので, よほど狭い帯域に制限するのでなければ普通にやった方がお手軽かつ高速な可能性もあります. # でもって, そこまで狭帯域に制限するなら DFT の方が速かったりする.

回答No.1

あんまりソース見ていないですけど、 (1)帯域制限をしてから(2)IFFTするということですか? それなら簡単にできます。 FFTとDFTは演算量が違うだけで結果は同じですから。 圧縮では普通に使われる技術ですね(DCTだったりもしますが)。 それとも、 IFFTを実行しつつ帯域制限もして(制限した帯域の)無駄な演算も行わないということでしょうか? やったこと無いですが、めんどうそうですね、 制限する帯域によっては効率よくできそうですが、普通のIFFT演算ではなくなりそうです。

kiku570827
質問者

お礼

IFFTを実行しつつ帯域制限もして(制限した帯域の)無駄な演算も行わないということです。帯域制限位相限定相関法(BLPOC)というものに使用します。IFFTを用いてやりたいのですがコード作成が難しくてできません。

関連するQ&A

  • 高速フーリエ変換(FFT)のプログラム

    高速フーリエ変換(FFT)のプログラム 高速フーリエ変換のプログラムが必要になり、visual c++で作成しています。 http://www.kiso.tsukuba.ac.jp/~watanabe/FFT.htm こちらのサイトの一番下のプログラムを参考にして作成しています。 ここで質問なのですが信号データが8点のプログラムが上記のサイトにはあり、私は作成したプログラムが正しいかどうかチェックするため4点の信号データのFFTを求めたいのですが、以下のようにx0=2,x1=2,x3=8,x4=-4の値をフーリエ変換するプログラムにしたところ、この4点フーリエ変換の正規の値と、このプログラムの出力値が違ってしまいます。私の知識不足も否めませんがどこが間違っているのかわからずこまっています。 よろしければ教えてください。 void four1(float datar[],float datai[],unsigned long nn, int isign){ unsigned long n,mmax,m,j,istep,i; double wtemp,wr,wpr,wpi,wi,theta; float tempr,tempi; n=nn; j=0; for (i=0;i<n ;i++) { if (j > i) { SWAP(datar[j], datar[i]); SWAP(datai[j], datai[i]); } m=n >> 1; while (m >= 1 && j > m) { j -= m; m >>= 1; } j += m; } mmax=1; while (n > mmax) { istep=mmax << 1; theta=isign*(6.28318530717959/istep); wtemp=sin(0.5*theta); wpr = -2.0*wtemp*wtemp; wpi=sin(theta); wr=1.0; wi=0.0; for (m=0;m<mmax;m++) { for (i=m;i<n;i+=istep) { j=i+mmax; tempr=wr*datar[j]-wi*datai[j]; tempi=wr*datai[j]+wi*datar[j]; datar[j]=datar[i]-tempr; datai[j]=datai[i]-tempi; datar[i] += tempr; datai[i] += tempi; } wr=(wtemp=wr)*wpr-wi*wpi+wr; wi=wi*wpr+wtemp*wpi+wi; } mmax=istep; } } void realft() { float datar[4]={2.0,2.0,8.0,-4.0}; float datai[4]={0.0,0.0,0.0,0.0}; int i; four1(datar,datai,4,-1); for(i=-1;i<3;i++){ printf("i = %d:\tReal: %f, \tImag: %f\n", i, datar[i]/4,datai[i]/4); } }

  • C言語でのFFTについて

    http://tsuyu.cocolog-nifty.com/blog/2007/03/publi.html に掲載されているVBAのFFTプログラムをC言語に書き換えて 実行しているのですがうまくいきません。 どこが、間違っているか教えてください。 ======以下FFTのサブルーチンソースコード===== void FFT(float Xr[], float Xi[]) { i=0,j=0,k=0,l=0,m=0,n=0,p=0,q=0; n=DN; m = log10(n)/log10(2); Table(c,s); l=n,h=1; for(g=1;g<=m;g++){ l/=2,k=0; for(q=1;q<=h;q++){ p=0; for(i=k;i<=l+k-1;i++){ j=i+l; a=Xr[i]-Xr[j], b=Xi[i]-Xi[j]; Xr[i] = Xr[i] + Xr[j], Xi[i] = Xi[i] + Xi[j]; if(p==0){ Xr[j]=a,Xi[j]=b; }else{ Xr[j] = a * c[p] + b * s[p], Xi[j] = b * c[p] - a * s[p]; } p+=h; } k+=l+l; } h+=h; } j=n/2; for(i=1;i<=n-1;i++){ k=n; if(j<i){ //ビットリバース swap(&Xr[i],&Xr[j]); swap(&Xi[i],&Xi[j]); } k=k/2; while(j>=k){ j=j-k; k /=2; } j = j + k; } }

  • Cプログラムです

    #include<stdio.h> #include<math.h> int main(void) { double D[512]; short data[16000]; double datafft[512]; double horizon[256]; double dB[256]; double power[256]; double Xr[512],Xi[512]; double w, pai=3.1415926; double er, ei; double r,x,r1,x1,r2,x2 ; int i,k,n,N=500; for(k = 0; k < N; k++) { Xr[k] = 0.0; Xi[k] = 0.0; for(n = 0; n < N; n++) { w = 2 * pai * k * n / N; er = cos(w); ei = sin(-w); r = r1 * r2 - x1 * x2; x = r1 * x2 + r2 * x1; Xr[k] = Xr[k] + er / N; Xi[k] = Xi[k] + ei / N; } } for(i=0;i<512;i++) { datafft[i] = (double)data[i+1000]; power[i] = Xr[i] * Xr[i] + Xi[i] * Xi[i]; } for(i=0;i < 256;i++) { horizon[i] = 16000.0 * ((double)i /512); dB[i] = 10 * log10(power[i]); printf("%f %f \n", horizon[i], dB[i]); } } 正規乱数を発生させたいんですけど、変な値が出てしまいます。 正しくはどうすれば良いか教えてもらえませんか?

  • FFT(高速フーリエ変換)の考え方

    ///////////質問(1)//////////////////// 離散フーリエ変換が以下の式で表せることは理解できます。      N-1       nk          -2πi / N A =  Σ  a   W ,   W = e    n   k=0   k                 n:基本周波数のn倍 N:データ数 ここで、実際のプログラム上では、 (1)akに サンプリングしたデータとサンプリング周期を掛けたもの を入れ、 (2)それをWと掛けあわせて行くと思います。 では、Wには何を代入すればいいのでしょうか? Wを実部・虚部にわける、もしくは、 絶対値・位相にわけて計算する方法で あっているのでしょうか? Anは周波数n成分の面積を表していると思うのですが、面積に虚数概念がで てくるのも変な話なような気がして、質問させていただきました。 ///////////質問(2)////////////////////    N/2-1    2nk   n N/2-1     2nk A = Σ a  W + W  Σ  a  W  n  k=0  2k        k=0  2k+1 Cooley-Tukey 型 FFTは質問(1)の式を上記のように2項に分解することで、 計算を減らすことができるとのことです。 しかし、私には理解できません。 左側の項でN/2回の掛け算を行い、右側の項でもN/2回の掛け算を行うのでは、 結局式を分解しただけで、何も変わっていないように思えます。 どのように考えればいいのでしょうか? アドバイスよろしくお願いします。

  • FFTのソースについて

    512個のデータに対してFFTをかけ、変換後の絶対値を得たいのですが矢印の部分でプログラムがとまり、特にエラーメッセージも出ません。 原因がわかりません。申し訳ありませんが、よろしくお願いします。 void FFT(int x[], int y[]) {    int d, dn, pow, m, j1, j2, exp2, j , k, ndv2, t, flags= 1;    int n_data = 0, dnumb = 512, nf, jk;    double w, arg, c, s, t1, t2, t0, dt, ana, anb, answer;    double [] xr = new double[512];    double [] xi = new double[512]; dn=dnumb; w=6.283185303/dnumb; pow=exp2(dnumb); *2の何乗かを計算するメソッドです for(t = 0; t < 512; t++){ if(t < 500) xr[t] = x[t]; else xr[t] = xr[499]; xi[t] = 0; } for(int i=1; i<=pow ; i++) { m=dn; dn=dn/2; arg=0; for(int j=1; j<=dn; j++) { c=Math.cos(arg); s=-flags*Math.sin(arg); arg=arg+w; k=m; while(k<=dnumb) { j1=k-m+j; j2=j1+dn; →→→→→→→  t1=xr[j1]-xr[j2]; t2=xi[j1]-xi[j2]; xr[j1]=xr[j1]+xr[j2]; xi[j1]=xi[j1]+xi[j2]; xr[j2]=c*t1+s*t2; xi[j2]=c*t2-s*t1; k=k+m; } } w=2*w; } j=1; ndv2=dnumb/2; for(int i=1; i<=(dnumb-1); i++) { if(i<j) { t1=xr[j]; t2=xi[j]; xr[j]=xr[i]; xi[j]=xi[i]; xr[i]=t1; xi[i]=t2; } k=ndv2; while(k<j) { j=j-k; k=k/2; } j=j+k; } for(t = 0; t < 512; t++){ ana = Math.pow(xr[t], 2); anb = Math.pow(xi[t], 2); answer = Math.pow(ana + anb, 0.5); } if (flags==1) { for(int j=1; j<=dnumb; j++) { xr[j]=xr[j]/dnumb; xi[j]=xi[j]/dnumb; } } }

    • ベストアンサー
    • Java
  • FFTを使って信号を周波数変換する方法を教えてください。

    音信号を周波数変換するプログラムを作成しています。 FFTを使って実現しようとしているのですが、うまくいきません。 現時点で作った方法では、 (1)FFTする。結果は配列x_re[NFFT]、x_im[NFFT]に格納。(x_re:実数部 x_im:虚数部 NFFT:ポイント数、配列の内容としては周波数の低い順に結果データが並んでいる) (2) x_re,x_imの内容をずらす。(例えば以下のように配列内容を1つずらせばIFFTをした時周波数が高くなるはず) for(i=0;i<NFFT-1;i++){ x_re[i+1]=x_re[i]; x_im[i+1]=x_im[i]; } (3)配列x_re[NFFT]、x_im[NFFT]に対しIFFTする。 FFT・IFFTが正しく動作するのは確認しています。 (動作実績も結構あります。(2)を省略し(1)(3)だけとすれば出力結果は入力結果と同じ(出力音声を聴いた感じで)になるので、正しく動作していると思います。) よって(2)が間違っていると思います。 (2)をどのようにすれば周波数変換できるのか教えてください。また今の所FFTを使おうとしているのですが、別に入力信号を周波数変換できればOK(周波数を上げたり下げたりしたい)なので、その方法があれば教えてください(><)

  • FFTをc言語でプログラミング

    助けてください↓c言語でFFTをプログラミングしているのですがバグが発見できずに困り果て居ます。 プログラミングは以下の通りです。画像はカラーではなく、グレースケールのMANDRILLに対して行っています。 for(n=0;n<M;n++){    for(i=0;i<N;i++){ for(j=0;j<N;j++){ Imbatx[n][i][j]=0; } } } //行FFT //まずはソート for(i=0;i<N;i++){ int x=1; for(j=0;j<N;j++){ if(j<N/2){        //左半分 if(j%2!=0) //jが奇数 batx[0][i][j]=imagein[i][j+N/2-1]; if(j%2==0) //jが偶数 batx[0][i][j]=imagein[i][j]; } if(N/2<=j && j<N){ //右半分 if(j%2!=0) //jが奇数 batx[0][i][j]=imagein[i][x+N/2-1];      if(j%2==0) //jが偶数 batx[0][i][j]=imagein[i][x]; x++; } } } for(i=0;i<N;i++){ for(n=1;n<M+1;n++){ k=0; for(j=0;j<N;j++){ mul=(int)pow(2.0,(double)n); if(j%mul < mul/2){ //偶数列 batx[n][i][j] = batx[n-1][i][j]+cos(2*PI*k/mul)*batx[n-1][i][j+mul/2]- sin(2*PI*k/mul)*Imbatx[n-1][i][j+mul/2]; Imbatx[n][i][j]=Imbatx[n-1][i][j]+sin(2*PI*k/mul)*batx[n-1][i][j+mul/2]+ cos(2*PI*k/mul)*Imbatx[n-1][i][j+mul/2]; k++; } else{ //奇数列 batx[n][i][j] = cos(2*PI*k/mul)*batx[n-1][i][j] + batx[n-1][i][j-mul/2]- sin(2*PI*k/mul)*Imbatx[n-1][i][j]; Imbatx[n][i][j]=Imbatx[n-1][i][j-mul/2]+sin(2*PI*k/mul)*batx[n-1][i][j]+cos(2*PI*k/mul)*Imbatx[n-1][i][j]; k++; } Re_countx[i][j] = batx[M][i][j]; Im_countx[i][j] = Imbatx[M][i][j]; } } } for(i=0;i<N;i++){ for(j=0;j<N;j++){ FFT[i][j]=pow(pow(Re_countx[i][j],2)+pow(Im_countx[i][j],2),0.5); FFT[i][j]=log10(FFT[i][j]); if(FFT[i][j]>max) max=FFT[i][j]; imageout[i][j]=(unsigned char)255*FFT[i][j]/max; } } 結果的に出力される画像は以下の画像になります。 ごらんのように縞模様が出てしまいます。 分かる方どうかお願いいたします。

  • 一次元のフーリエ変換(FFT)

    フーリェ変換について勉強しているものです。このページ(​http://www.kurims.kyoto-u.ac.jp/~ooura/fftman/ftmndl.html)で配布している、fft.tgz (71 KB)というファイルのdfst()という関数を使おうと思っているんですが、その関数の返り値(返り配列と言った方が正しい?)の詳しい意味がよくわからないのです。 以下はこの関数のマニュアル(readme.txt)の一部です。 dfst()の大まかな処理内容を示しています。 --------------------------------------------------------- 以下の変換は : S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n このように分割される : S[2*k] = sum_j=1^n/2-1 (a[j]-a[n-j])*sin(pi*j*k/(n/2)), ◎ S[2*k+1] = sum'_j=1^n/2 (a[j]+a[n-j])*sin(pi*j*(k+1/2)/(n/2)) 注意:sum_■^▲はΣを表しており、■から▲までの部分和となっています。 nは要素数です。 ------------------------------------------------------ この◎の配列はどのような意味を持つのでしょうか?

  • 離散フーリエ変換について

    G(n/Nτ)= Σ[k=0->N-1] {τ*f(kτ)*e^(-i2πkn/N)} 上記の式は離散フーリエ変換の式らしいですが、 これを関数化しようとしてつまずいています。 どのように解釈すれば、関数化できるのでしょうか? 特に、複素数iがよくわかりません。 e^iが点を左に90度回転させるくらいはわかります。 e^(-i2πkn/N)を関数powで表現できなくて困っていますが、多分見当違いだとは思います。 フーリエ級数展開はわかります。 最終的にはFFTを行いたいのですが、 その理解の前に離散フーリエ変換ができないといけないと思っています。 よろしくお願いします。 void GraphClass::ScatteredFourierConvert() { /* N-1 G(n/Nτ)= Σ {τ*f(kτ)*e^(-i2πkn/N)} k=0 k=何番目の値か? τ=値を読んだ間隔 N=値を読んだ数 */ int Tau=10; int N=RangeSize.cx*2/Tau; double Sum=0; double e=2.7182818; for(int k=0;k<N-1;k++) { Sum+=Tau*Data[k*Tau]*pow(e,-i*2*PI*k*n/N); } }

  • 離散フーリエ変換について

    double DiscreteFourierTransformClass::DiscreteFourierTransformProcess(int n) { double e=2.7182818; double Pai=3.1415926535; int Tau=10; int N=DataLen/Tau; double Sum=0; for(int k=0;k<N;k++) { Sum+=Tau*Data[k*Tau]*pow(e,i*2*Pai*k*n/N); } return Sum; } Tau*Data[k*Tau]=帯 pow(e,i*2*Pai*k*n/N)=ほしい周波数の乗算 ということで、離散フーリエ展開と同じなのでしょうか? パソコン上で計算するには、離散フーリエ展開で計算して、 FFTというのは、pow(e,i*2*Pai*k*n/N)この回転子が3,4個程度のcos,sinで決まるから その値で計算すれば高速になるよということでしょうか? 離散フーリエ変換 G(n/N)=Σ[k=0::N-1]{τ*f(kτ)e^-i2πkn/N} k:何番目の値か τ:値を読んだ間隔 N:値を読んだ数 n:基本周波数の何倍か G(n/N)->ほしいnの成分を得る