高速フーリエ変換(FFT)のプログラム
- 高速フーリエ変換(FFT)のプログラムをvisual c++で作成中
- 信号データのFFTを求めたいが、プログラムの出力値が正規の値と違ってしまう
- どこが間違っているのかわからず困っている
- ベストアンサー
高速フーリエ変換(FFT)のプログラム
高速フーリエ変換(FFT)のプログラム 高速フーリエ変換のプログラムが必要になり、visual c++で作成しています。 http://www.kiso.tsukuba.ac.jp/~watanabe/FFT.htm こちらのサイトの一番下のプログラムを参考にして作成しています。 ここで質問なのですが信号データが８点のプログラムが上記のサイトにはあり、私は作成したプログラムが正しいかどうかチェックするため４点の信号データの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); } }
- MRHANABI
- お礼率0% (0/1)
- Visual Basic
- 回答数1
- ありがとう数2
- みんなの回答 （1）
- 専門家の回答
質問者が選んだベストアンサー
while (m >= 1 && j > m) { じゃなくて while (m >= 1 && j >= m) { だろうね。 それから for(i=-1;i<3;i++){ これは for(i=0;i<4;i++){ だろうね。
関連するQ&A
- プログラム内への他プログラムの導入方法
前回は不適切な質問を出してしまい、誠に申し訳ありませんでした。 関係者各位様には、この場を借りて謝らせてください。 改めて質問させてください。 あるプログラムＡの中に別のプログラムＢを組み込みたい(結果を代入できれば十分です)のですが、プログラムＢが自分で作ったプログラムでないため、いまいち勝手がつかめません。 プログラムＢ外で、この結果を使いたい場合、どのようにプログラムを組めばいいか、教えて欲しいです。 よろしければ、今一度皆様の知識をお貸しください。 プログラムＢ→ #include <stdio.h> #include <math.h> #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr void four1(float data[], unsigned long nn, int isign); /* test driver for four() */ int main(void) { /* sample data */ /* Note:sample[0] is not used.*/ /* It is due to the implementation of the function four1().*/ float original[]={0, /* not used*/ 0,0, /* 0 + 0i */ 1,0, /* 1 + 0i */ 2,0, /* 2 + 0i */ 3,0}; /* 3 + 0i */ /* fourier-transformed original[] */ float fft[]={0,0,0,1,0,2,0,3,0}; /* inverse-fourier-transformed original[]*/ float invfft[]={0,0,0,1,0,2,0,3,0}; int i; int nn = 4; /* nmber of sample data */ printf("Original data:\n"); for (i=1; i<2*nn; i+=2) { printf("%d th datum: %8.5f + %8.5fi \n", (i+1)/2, original[i], original[i+1]); } four1(fft, nn, 1); /* Fourier transform */ printf("Fourier Transformed data:\n"); for (i=1; i<2*nn; i+=2) { printf("%d th datum: %8.5f + %8.5fi \n", (i+1)/2, fft[i], fft[i+1]); } four1(invfft, nn, 1); /* Fourier transform */ four1(invfft, nn, -1); /* Inverse Fourier transform */ printf("Fourier Transformed data:\n"); for (i=1; i<2*nn; i+=2) { printf("%d th datum: %8.5f + %8.5fi \n", (i+1)/2, invfft[i]/nn, invfft[i+1]/nn); } four1(fft, nn, 1); /* Fourier transform */ printf("Fourier Transformed data:\n"); for (i=1; i<2*nn; i+=2) { printf("%d th datum: %8.5f + %8.5fi \n", (i+1)/2, fft[i], fft[i+1]); } } void four1(float data[], 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 << 1; j=1; for (i=1;i<n;i+=2) { if (j > i) { SWAP(data[j],data[i]); SWAP(data[j+1],data[i+1]); } m=nn; while (m >= 2 && j > m) { j -= m; m >>= 1; } j += m; } mmax=2; while (n > mmax) { istep=mmax << 1; theta=isign*(6.28318530717959/mmax); wtemp=sin(0.5*theta); wpr = -2.0*wtemp*wtemp; wpi=sin(theta); wr=1.0; wi=0.0; for (m=1;m<mmax;m+=2) { for (i=m;i<=n;i+=istep) { j=i+mmax; tempr=wr*data[j]-wi*data[j+1]; tempi=wr*data[j+1]+wi*data[j]; data[j]=data[i]-tempr; data[j+1]=data[i+1]-tempi; data[i] += tempr; data[i+1] += tempi; } wr=(wtemp=wr)*wpr-wi*wpi+wr; wi=wi*wpr+wtemp*wpi+wi; } mmax=istep; } } #undef SWAP
- ベストアンサー
- C・C++・C#
- Cプログラムによる画像の高速フーリエ変換FFT
Cプログラムにる画像処理について教えてください。逆高速フーリエ変換（IFFT)によって周波数領域から実空間領域に変換する際、周波数帯域を（低周波数領域に）制限してIFFTを実行できるという文献があるのですが可能でしょうか？通常の逆離散フーリエ変換（IDFT)では可能でした。 以下のコードの修正で可能でしょうか？ 例） ５１２＊５１２の空間周波数領域データを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); }
- 締切済み
- C・C++・C#
- フーリエ変換やＦＦＴのプログラム
Ｃ言語で書かれているＦＦＴやフーリエ変換のプログラムのあるお勧めのサイトがあればおしえていただけないでしょうか？ よろしくお願いいたします。
- ベストアンサー
- C・C++・C#
- FFT（高速フーリエ変換）の考え方
///////////質問(1)//////////////////// 離散フーリエ変換が以下の式で表せることは理解できます。 N-1 nk -2πi / N Ａ ＝ Σ ａ Ｗ ， Ｗ ＝ ｅ 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 Ａ ＝ Σ ａ Ｗ ＋ Ｗ Σ ａ Ｗ n k=0 2k k=0 2k+1 Cooley-Tukey 型 FFTは質問(1)の式を上記のように２項に分解することで、 計算を減らすことができるとのことです。 しかし、私には理解できません。 左側の項でN/2回の掛け算を行い、右側の項でもN/2回の掛け算を行うのでは、 結局式を分解しただけで、何も変わっていないように思えます。 どのように考えればいいのでしょうか？ アドバイスよろしくお願いします。
- 締切済み
- 数学・算数
- 多次元高速フーリエ変換について
高速フーリエ変換fftによって、計算量のオーダーが n^2 からnlogn まで落とせるんですよね？ それで、3次元のフーリエ変換って、 1次元のフーリエ変換を３回やれば n^2*nlogn=n^3lognのオーダーでできると思うのですが、 これ以上速いオーダーではできませんか？
- ベストアンサー
- 数学・算数
- 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・C++・C#
- ＦＦＴをc言語でプログラミング
助けてください↓ｃ言語でＦＦＴをプログラミングしているのですがバグが発見できずに困り果て居ます。 プログラミングは以下の通りです。画像はカラーではなく、グレースケールの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; } } 結果的に出力される画像は以下の画像になります。 ごらんのように縞模様が出てしまいます。 分かる方どうかお願いいたします。
- 締切済み
- C・C++・C#
- 離散フーリエ変換について
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); } }
- ベストアンサー
- C・C++・C#
補足
早速の回答ありがとうございます。 値が合いました。ありがとうございます。 一つ質問なのですが、x0～x3までの信号データをフーリエ変換して求まる周波数スペクトル値はX-1～X2で与えられると教科書等に書いてあります。 iの値が0～3をとるとするとX-1の値がi=3のときに与えられますよね。 この場合X-1＝1/4*{x0+jx1-x2-jx3}=X3 なのでしょうか。X-1がi=3の場合に対応するだけなのでしょうか。 理解が浅く意味がわからない質問になってしまっているかもしれませんがよろしくお願いします。