• 締切済み

一次元のフーリエ変換(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は要素数です。 ------------------------------------------------------ この◎の配列はどのような意味を持つのでしょうか？

みんなの回答

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

えと.... 質問の意味がよくわからないんですが, 単に「入力を離散サイン変換した結果」じゃダメ?

0-o
質問者

補足

返答ありがとうございます。 私はまだ勉強中ですので詳しいことは言えませんが、おっしゃったことはズバリその通りなんですが、例えば音声のデータを対象としたとき、変換結果の配列からどの周波数の音が一番大きいかな～、なんて見るわけですが、それなら 「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))」 でもいいわけです。 私が分からないのは、この二つを比較して何が違うのか、または何の目的でこの形にするのか、です。

  1. カテゴリ
  2. [技術者向] コンピューター
  3. プログラミング・開発
  4. C・C++・C#

関連するQ&A

  • フーリェ変換（特に一次元のFFT）

    フーリェ変換について勉強しているものです。このページ（http://www.kurims.kyoto-u.ac.jp/~ooura/fftman/ftmndl.html）で配布している、fft.tgz (71 KB)というファイルの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); } }

  • 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)について

    二次元のフーリエサイン変換を考えています。 正方形の中心にだけ高い数値を持つ状況をつくり、その波数分布を示したい(最終目標は違いますが・・・）のですがうまくいきません。 下記のように書いてみましたが、コンパイルは通っても、変換した値が一様に0となってしまいます。 よろしければ、見てやってください。 #include<stdio.h> #include<math.h> #define max 100 #define pi 3.141592 main() { //変数i,jが位置、ki,kjが波数に当たります int i,j,ki,kj; //Nはデータ点の数 int N=30; double a[max][max]; double A[max][max]; //aの分布をデルタ関数的なものとします。 for(i=0;i<N;i++){ for(j=0;j<N;j++){ a[i][j]=0; }} a[N/2+1][N/2+1]=100; for(ki=0;ki<N;ki++){ for(kj=0;kj<N;kj++){A[ki][kj]=0.0; for(i=0;i<N;i++){ for(j=0;j<N;j++){ A[ki][kj] += (2/N)*(2/N)*a[i][j]*sin((2*pi*ki*i)/N)*sin((2*pi*kj*j)/N); }}}} return(0); }

  • ＦＦＴを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; } } 結果的に出力される画像は以下の画像になります。 ごらんのように縞模様が出てしまいます。 分かる方どうかお願いいたします。

  • 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回の掛け算を行うのでは、 結局式を分解しただけで、何も変わっていないように思えます。 どのように考えればいいのでしょうか？ アドバイスよろしくお願いします。

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

    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の成分を得る

  • 逆離散フーリエ変換

    逆フーリエ変換でもとの波形を作りたいです。 http://www.geocities.jp/supermisosan/fourier.html を参考にしました。 ///参考URLの逆フーリエ変換プログラム #include <stdio.h> #include <math.h> #define max 10000 //Nの最大限度数 #define pi 3.1415926535 //円周率 int main() { int i,k,n,N; double Ref,Imf; double ReF[max+1],ImF[max+1]; FILE *fourier; FILE *inversef; //逆フーリエ変換したいデータをあらかじめfourier.dataに保存しておく fourier=fopen("fourier.data","r"); inversef=fopen("inversef.data","w"); //データの読み込み。 for(N=0;N<max;N++) { if(fscanf(fourier,"%d %lf %lf",&i,&ReF[N],&ImF[N]) == EOF) { N--; break; } } //フーリエ逆変換 for(k=0;k<N;k++) { Ref=Imf=0.0; for(n=0;n<N;n++) { Ref+=(ReF[n]*cos(2*pi*k*n/N)-ImF[n]*sin(2*pi*k*n/N))/N; Imf+=(ReF[n]*sin(2*pi*k*n/N)+ReF[n]*sin(2*pi*k*n/N))/N; } fprintf(inversef,"%d %f %f\n",k,Ref,Imf); } fclose(fourier); fclose(inversef); return 0; } ////// Refが実数部分の値で、Imfが虚数部分の値になるのですが、 どうしたら元の波形データが得られるのかがわかりません。 教えて頂けたら幸いです。

  • フーリエ変換のC言語プログラムについて

    正弦波（およびガウス性雑音)をフーリエ変換(離散)→逆フーリエ変換するというプログラムを組みました。正弦波をフーリエ変換すると実部は２回ピークがくるはずなのですが、すべて「0.000000」または「-0.000000」と表示されてしまいます。虚部は正常なようで実装の仕方もさほど違わないので、何が問題なのかわからずにいます。念のためコードはすべて載せますが、該当箇所は関数Fourierの fp = fopen("reX.txt", "w"); //書き込む あたりです。問題点を教えていただけないでしょうか。お願いします。 //gauss.txt, sin.txt:発生させたガウス性雑音＆正弦波 //reX, imX:フーリエ変換の実部と虚部 //re-1, im-1:逆フーリエ変換の実部と虚部 #include <stdio.h> #include <stdlib.h> #include <math.h> #include <time.h> #define PI 3.14159265358979323846 #define N 256 //n:長さ, w:角周波数, p:位相(phase), a:振幅 void SinCurve(int n, double w, double p, double a) { FILE *fp; double x; int t; fp = fopen("sin.txt", "w"); //書き込むので"w" if(fp == NULL) { printf("file open error\n"); } else { for(t = 0; t < n; t++) { x = a * sin( w*(double)t + p ); fprintf(fp, "%f\n", x); } } fclose(fp); } //n:長さ, s:分散, m:平均 void Gauss(int n, double s, double m) { FILE *fp; double x, x1, x2, y1; int t; srand((unsigned) time(NULL)); fp = fopen("gauss.txt", "w"); //書き込むので"w" if(fp == NULL) { printf("file open error\n"); } else { for(t = 0; t < n; t++) { x1 = ( (double)rand() + 1.0 ) / ( (double)RAND_MAX + 2.0); x2 = ( (double)rand() + 1.0 ) / ( (double)RAND_MAX + 2.0); y1 = pow(-2.0*log(x1), 0.5) * cos(2.0*PI*x2); y1 = s * y1 + m; fprintf(fp, "%f\n", y1); } } fclose(fp); } //ファイル名sのデータをフーリエ変換し、実部と虚部をreX, imXに保存 void Fourier(int num, char *s) { FILE *fp; int k, n; double largeX, x[N+100], t; fp = fopen(s, "r"); //読み込み if(fp == NULL) { printf("file open error\n"); } else { // printf("%s\n", s); for(k = 0; k < num; k++) { fscanf(fp, "%lf", &x[k]); printf("x[%d]=%f\n", k, x[k]); } } fp = fopen("reX.txt", "w"); //書き込む if(fp == NULL) { printf("file open error\n"); } else { for(k = 0; k < num; k++) { largeX = 0.0; t = 2.0*PI*(double)k / (double)N; for(n = 0; n < num; n++) { largeX += x[n] * cos((double)n*t); // printf("%f\n", largeX); } fprintf(fp, "%f\n", largeX); printf("reX[%d]=%f\n", k, largeX); } } fp = fopen("imX.txt", "w"); //書き込む if(fp == NULL) { printf("file open error\n"); } else { for(k = 0; k < num; k++) { largeX = 0.0; t = 2.0*PI*k / (double)N; for(n = 0; n < num; n++) { largeX -= x[n] * sin(n*t); } fprintf(fp, "%f\n", largeX); } } fclose(fp); } void InverseFourier(int num) { FILE *fp; int k, n; double a[N+100], b[N+100], x, t; //a:reX, b:imX fp = fopen("reX.txt", "r"); //読み込み if(fp == NULL) { printf("file open error\n"); } else { for(k = 0; k < num; k++) { fscanf(fp, "%lf", &a[k]); // printf("a[%d]=%f\n", k, a[k]); } } fp = fopen("imX.txt", "r"); //読み込み if(fp == NULL) { printf("file open error\n"); } else { for(k = 0; k < num; k++) { fscanf(fp, "%lf", &b[k]); // printf("b[%d]=%f\n", k, b[k]); } } fp = fopen("re-1.txt", "w"); //読み込み if(fp == NULL) { printf("file open error\n"); } else { for(n = 0; n < num; n++) { x = 0.0; t = 2.0*PI*(double)n / (double)N; for(k = 0; k < num; k++) { x +=a[k] *cos(k*t) - b[k] *sin(k*t); } x /= (double)N; fprintf(fp, "%f\n", x); // printf("x[%d]=%f\n", n, x); } } /* fp = fopen("im-1.txt", "w"); //読み込み if(fp == NULL) { printf("file open error\n"); } else { for(n = 0; n < num; n++) { x = 0.0; for(k = 0; k < num; k++) { t = 2.0*PI*(double)k / (double)N; x = x + a[k] *sin(n*t) + b[k] *cos(n*t); } x /= (double)N; fprintf(fp, "%f\n", x); } } */ fclose(fp); } int main(void) { SinCurve(N, PI/8.0, 0.0, 1.0); // Gauss(N, 1.0, 0.0); Fourier(N, "sin.txt"); // Fourier(N, "gauss.txt"); InverseFourier(N); return 0; }

  • EXCELを使ってFFT(高速フ－リエ変換)をする方法2

    添付ファイルの波形のスペクトル及びパワースペクトルを求め、それらの違いについて考察せよ。ただし、データ数1024個でFFTを行うこと。 1.波形１：sin波　１周期（振幅1） A:1(A1)-1024(A1024)の番号が並んでいる B：B1=SIN(2*PI()*(A1-1)/1024+PI()/4)～B1024=SIN(2*PI()*(A1024-1)/1024+PI()/4) この場合EXCELでどのような操作を行うとFFT(高速フ－リエ変換)をした波形になるのでしょうか？結果は複素数になるのでグラフ化のとき何かしなければならないと思うのですが詳しい方教えてください．お願いします. 波形のスペクトルとパワースペククトルの違いもよくわかりません。

注目のQ&A

カテゴリ

専門家に質問してみよう

職業から探して質問する

あなたにピッタリな商品が見つかる！ OKWAVEセレクト

質問する