c言語でDFTのプログラムを作成し、データ数N個に対して処理時間を測定した結果

このQ&Aのポイント
  • c言語でDFTのプログラムを作成し、データ数N個に対して処理時間を測定した結果、処理が早すぎる問題が発生している。
  • 元信号には5sin(20πt)の値を入れており、フーリエ変換を行うと対称なデータが得られる。
  • 実行した結果、データ数N/2を中心に対称なデータが得られるものの、処理が早すぎて正常に終了しているか不安。
回答を見る
  • ベストアンサー

c言語でDFTのプログラムを作成したのですが

c言語でDFTのプログラムを作成しました。 以下にソースを載せます。 #include<stdio.h> #include<math.h> #include<stdlib.h> #include<time.h> #define PI 3.141592653589793 #define N 64 //データ数 DFT(double result[]){ int i,k; double A[N],B[N],T=0; //A[N]:実数部,B[N]:虚数部 double a,b; for(k=0;k<N;++k){ a=b=0; for(i=0;i<N;++i){ a=a+result[i]*cos(2*PI*i*k/N);      b=b+(-1.0)*result[i]*sin(2*PI*i*k/N); } A[k]=a/N; B[k]=b/N; } for(i=0;i<N;++i){ printf("[%f秒]:Re:%f,Im:%f\n",T,A[i],B[i]); //変換後の値を表示 T=T+(0.1/N); } } main(){ int i; double T=0; double Sampdata; double result[N]; for(i=0;i<N;++i){ Sampdata=5*sin(20*PI*T);      //0~0.1秒間をN個にサンプリング result[i]=Sampdata; //サンプリングデータを代入 T=T+(0.1/N); } clock_t start,end; //処理時間計測開始 start=clock(); DFT(result); end=clock(); printf("%.2f秒かかりました\n",(double)(end-start)/CLOCKS_PER_SEC); //処理時間表示 } 元信号には5sin(20πt)の値を入れています。この信号は周期は0.1secです。 これでフーリエ変換を行うとデータ数N/2を中心に対称なデータが出てくるのですが、処理が終わるのが早い気がするんです。 例えば2^15個のデータで実行しても2分もかからずに処理が終わってしまいます。一応、対称性が出てるとはいえ、終わるのが早すぎる気がするのですが、おかしい所があれば教えていただけると嬉しいです。 よろしくお願いします。

質問者が選んだベストアンサー

  • ベストアンサー
  • kmee
  • ベストアンサー率55% (1857/3366)
回答No.1

計算に使っているコンピュータの性能はどれくらいなのでしょう? このプログラムで一番時間がかかると思われるのは > for(i=0;i<N;++i){ > printf("[%f秒]:Re:%f,Im:%f\n",T,A[i],B[i]); //変換後の値を表示> > T=T+(0.1/N); > } この出力の部分です。 手許のcore i7 2.66GHz gcc -O0 (最適化無し) DFT(result); を2^15回ループ 結果出力部分をコメントアウト で22秒程でした。 CPUの速度、コンパイラの最適化等でもっと速くなります。 2分かからないくらいなら、別に変では無いと思います。 心配なら、この計算結果を、別の方法で計算したもの(実績のあるDFT計算ライブラリ、Excel等の別ソフト等)と比較してみては?

fcb9203
質問者

補足

回答ありがとうございます。 私の周りでは2^15でDFTを行った場合に30分も時間がかかったという意見があったのですがこんなにかかるものなのでしょうか? 私のと処理の仕方が違うかもしれませんが、それでも2^15のデータの処理に30分もかかるというのが自分と違い、不安になってこの質問をしました。 その周りの人にもプログラムを教えてもらえるような状況ではないのですがどうでしょうか? しつこいですがお願いします。

関連するQ&A

  • 離散フーリエ変換(DFT)のプログラム

    タイトルの通りのレポートを出されたのですが、その問題に似たサンプルソースすら何をいっているのかわからない状態です。ひとまず、サンプルソースが何をいっているのか理解したいので、いくつか教えてください。 ソースです。質問はその後に書かせてもらいました。 #include<stdio.h> #include<math.h> #define N 10 #define F 0.1 #define PI 3.14151692 #define SQR(x) ((x)*(x)) void func() { int n; FILE *fp; fp=fopen("temporal.data","w"); for(n=0;n<N;n++) fprintf(fp,"%lf\n",cos(2.0*PI*F*(double)n)); fclose(fp); } void get_data(double x[]) { int n; FILE *fp; fp=fopen("temporal.data","r"); for(n=0;n<N;n++) fscanf(fp,"%lf",&x[n]); fclose(fp); } void dft(double x[],double X_r[],double X_i[]) { int k,n; for(k=0;k<N;k++){ X_r[k]=X_i[k]=0.0; for(n=0;n<N;n++){ X_r[k]+=x[n]*cos(2.0*PI*(double)n*(double)k/(double)N); X_i[k]-=x[n]*sin(2.0*PI*(double)n*(double)k/(double)N); } } for(k=0;k<N;k++) printf("X[%d]=%lf+j%lf\n",k,X_r[k],X_i[k]); } void amplitude(double X_r[],double X_i[]) { int k; FILE *fp; double amp; fp=fopen("amp.data","w"); for(k=0;k<N;k++){ amp=sqrt(SQR(X_r[k])+SQR(X_i[k])); fprintf(fp,"%lf\n",amp); } fclose(fp); } main() { double x[N],X_r[N],X_i[N]; func(); get_data(x); dft(x,X_r,X_i); amplitude(X_r,X_i); } 文字数の制限があるみたいなので、質問を別にさせてもらいます。

  • C言語のプログラムで質問です。

    C言語のプログラムで質問です。 下のプログラム(最小二乗法の計算)を実行したところ -1.#IND00 というエラーが出てしまいます。 どこを直せばいいのでしょうか、教えてください。 #include <stdio.h> #include <math.h> /* gauss33.c */ #define N 3 main(){ double A[N][N],Aa[N][N]; double b[N],x[N], bb[N], e[N]; int n=N; int i, j, k; double akk, aik, s; double y[N]; double xx,yy; for(i=0;i<n;i++){ /*変数の初期化*/ x[i]=y[i]=0; for(j=0;j<n;j++) A[i][j]=0; } for(i=0;i<5;i++){ /*データ点は5点*/ printf("\n(x,y)="); scanf("%lf,%lf",&xx,&yy); A[0][0]+=xx*xx*xx*xx; /*Σx^4*/ A[0][1]+=xx*xx*xx; /*Σx^3*/ A[0][2]+=xx*xx; /*Σx^2*/ A[0][1]=A[1][0]; A[0][2]=A[1][1]=A[2][0]; A[1][2]+=xx; /*Σx*/ A[1][2]=A[2][1]; A[2][2]=n; y[0]+=xx*xx*yy; /*Σx^2y*/ y[1]+=xx*yy; /*Σxy*/ y[2]+=yy; /*Σy*/ } /* save original coefficients */ for(i=0; i<n; i++){ for(j=0; j<n; j++){ Aa[i][j]=A[i][j]; } bb[i]=b[i]; } /* forward operation */ for(k=0; k<n-1; k++){ akk=1/A[k][k]; for (i=k+1; i<n; i++){ aik=-A[i][k]*akk; for (j=k+1; j<n; j++){ A[i][j]+=aik*A[k][j]; } b[i]+=aik*b[k]; } for(j=k+1; j<n; j++){ A[k][j]*=akk; } b[k]*=akk; } /* backward operation */ x[n-1]=b[n-1]/A[n-1][n-1]; for(k=n-2; k>=0; k--){ s=0.0; for (j=k+1; j<n; j++){ s+=A[k][j]*x[j]; } x[k]=b[k]-s; } /* chek */ for(i=0; i<n; i++){ s=0.0; for(j=0; j<n; j++){ s+=Aa[i][j]*x[j];} e[i]=s-bb[i]; printf("\nx(%d)=%f error=%f\n",i, x[i], e[i]); } }

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

    正弦波(およびガウス性雑音)をフーリエ変換(離散)→逆フーリエ変換するというプログラムを組みました。正弦波をフーリエ変換すると実部は2回ピークがくるはずなのですが、すべて「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; }

  • C++言語のプログラムをfortranに変換

    const int N = 100; const double q = 10.0, dt = 0.00001, Dm = 30.0, t0 = 2.0, K = 1.0, pi = 3.141592, f = 3.0; double C[N], dC[N]; double dx = q/N; for (int i = 0; i < N; ++i) C[i] = 0; // 初期条件 for (double t = 0; t < t0; t += dt) {  C[0] = 0; // 境界条件1  C[N - 1] = K*sin(2*pi*f*t); // 境界条件2  for (int i = 1; i < N - 1; ++i) dC[i] = (Dm*(C[i + 1] - 2*C[i] + C[i - 1])/(dx*dx))*dt;  for (int i = 1; i < N - 1; ++i) C[i] += dC[i]; } for (int i = 0; i < N; ++i) cout << C[i] << endl; // t = t0 このプログラムをfortranに変換できる方いますか?

  • C言語のプログラムで質問です。

    C言語のプログラムで質問です。 これは、2元1次連立方程式の解を求めるプログラムです。 このプログラムを (1)3元1次連立方程式の解を求めるプログラムにする (2)係数行列、定数行列(6、7行目)をキーボードからの入力にする。 ようにしたいのですが、どうすればよいでしょうか。 前半の部分を変えれば良いようなのですが分かりません。教えてください。 #include <stdio.h> #include <math.h> /* gauss22.c */ #define N 2 main(){ double A[N][N]={1.,4.,3.,2.}, Aa[N][N]; /*簡単のため係数行列を予め指定*/ double b[N]={4.,5.}, x[N], bb[N], e[N]; /*簡単のため定数ベクトルを予め指定*/ int n=2; int i, j, k; double akk, aik, s; /* save original coefficients */ for (i=0; i<n; i++){ for (j=0; j<n; j++){ Aa[i][j]=A[i][j]; } bb[i]=b[i]; } /* forward operation */ for (k=0; k<n-1; k++){ akk=1/A[k][k]; for (i=k+1; i<n; i++){ aik=-A[i][k]*akk; for (j=k+1; j<n; j++){ A[i][j]+=aik*A[k][j]; } b[i]+=aik*b[k]; } for (j=k+1; j<n; j++){ A[k][j]*=akk; } b[k]*=akk; } /* backward operation */ x[n-1]=b[n-1]/A[n-1][n-1]; for (k=n-2; k>=0; k--){ s=0.0; for (j=k+1; j<n; j++){ s+=A[k][j]*x[j]; } x[k]=b[k]-s; } /* chek */ for (i=0; i<n; i++){ s=0.0; for (j=0; j<n; j++){ s+=Aa[i][j]*x[j];} e[i]=s-bb[i]; printf("x(%d)=%f error=%f?n",i, x[i], e[i]); } }

  • C言語についての質問です。

    C言語についての質問です。 このプログラムの前半にある A[N][N]={{1,0,0},{0,1,0},{0,0,1}} b[N]={4,5,3} という行列の成分をキーボードから入力するようにする にはどうすればいいでしょうか。 for や scanf や printf を使って、変えてくれないでしょうか。 #include <stdio.h> #include <math.h> /* gauss33.c */ #define N 3 main(){ double A[N][N]={{1,0,0},{0,1,0},{0,0,1}}; double b[N]={4,5,3}; double Aa[N][N]; double x[N], bb[N], e[N]; int n=N; int i, j, k; double akk, aik, s; /* input original coefficients */ /* save original coefficients */ for(i=0; i<n; i++){ for(j=0; j<n; j++){ Aa[i][j]=A[i][j]; } bb[i]=b[i]; } /* forward operation */ for(k=0; k<n-1; k++){ akk=1/A[k][k]; for (i=k+1; i<n; i++){ aik=-A[i][k]*akk; for (j=k+1; j<n; j++){ A[i][j]+=aik*A[k][j]; } b[i]+=aik*b[k]; } for(j=k+1; j<n; j++){ A[k][j]*=akk; } b[k]*=akk; } /* backward operation */ x[n-1]=b[n-1]/A[n-1][n-1]; for(k=n-2; k>=0; k--){ s=0.0; for (j=k+1; j<n; j++){ s+=A[k][j]*x[j]; } x[k]=b[k]-s; } /* chek */ for(i=0; i<n; i++){ s=0.0; for(j=0; j<n; j++){ s+=Aa[i][j]*x[j];} e[i]=s-bb[i]; printf("\nx(%d)=%f error=%f\n",i, x[i], e[i]); } }

  • C言語についてなのですが、

    C言語についてなのですが、 #include<stdio.h> #include<stdlib.h> #include<time.h> #include<search.h> int main(void) { int i,j,k,temp,n,count,time,list[65537]; clock_t startTime, endTime; printf("取得する乱数の個数を入力してください\n"); scanf("%d",&n); srand((unsigned)time(NULL)); printf("Before sort\n"); startTime = clock(); for(i = 0; i < n; i++) { list[i] = rand(); /* printf("%d\n", list[i]);*/ } count = 0; for (i = 1; i < n; i++) { for (j = i; j < n-i-1; j++) { count++; if(list[j] < list[j+1]) { temp = list[j]; list[j] = list[j+1]; list[j+1] = temp; } } } endTime = clock(); printf("\nAfter sort\n"); for(k = 0; k < n; k++) { /* printf("%d\n", list[k]);*/ } printf("\n比較回数:%d\n", count); printf("実行時間:%.4f秒\n", (double)(endTime - startTime) / CLOCKS_PER_SEC); return 0; } 上記のソースコードをcygwinで gcc -Wall -o k5-1-2 k5-1-2.c でコンパイルしようとすると k5-1-2.c:関数'main'内 k5-1-2.c:14:error:called object is not a function と表示されます。 いろいろなサイトを参考にして乱数取得用に srand((unsigned)time(NULL));を使うように書かれていたので使っているのですが、何かだめなのでしょうか?自分ではお手上げ状態で。

  • c言語 プログラム

    c言語の時間計測を調べていたら下記のようなプログラムを見つけました。 プログラムを読んでて疑問があったので質問します。 質問は2つあります for文のところに x = (double)k/loop が入っているのですがどのような意味なのでしょうか? また*1.e6の値はいくらなのでしょうか? サンプルプログラム #include<stdio.h> #include<time.h> #include<math.h> int main(){ clock_t stt,mid,end;//_測定時刻 double secs;//_経過秒数 double msec;//_マイクロ秒 int k,loop=10000000;//_測定増幅ループ double x;//_対象関数の引数 int clksec=CLOCKS_PER_SEC; char *ttl="elapse_for_exp(x)"; int n,m=7;//_観測ループ double sm,ss;//_合計、平方和 double mean;//_平均値 double stdv;//_標準偏差 sm=ss=0; for(n=0;n<m;n++){ stt= clock();//_測定開始 for(k=0;k<loop;k++){ x = (double)k/loop; } mid = clock();//_中間測定 for(k=0;k<loop;k++){ x=(double)k/loop; exp(x);//_時間測定対象の関数 } end=clock();//_測定終了 secs=(double)((end-mid)-(mid-stt))/clksec; msec=secs/loop*1.e6; printf("%s = %fmicrosec.\n",ttl,msec); sm += msec; ss += msec*msec; } printf("合計____%f\n",sm); printf("平方和__%f\n",ss); mean = sm/m; stdv=sqrt(ss/m-mean*mean); printf("平均値___=_%f\n",mean); printf("標準偏差_=_%f\n",stdv); return 0; }

  • 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]); } } 正規乱数を発生させたいんですけど、変な値が出てしまいます。 正しくはどうすれば良いか教えてもらえませんか?

  • C言語

    #include <stdio.h> #include <stdib.h> int main (void){ double a[5]={0.0,4.0,0.0,-5.0,1.0}; double x; int i,j,k,n; n=4; x=0.75; for(i=1;i<=n;i++) printf("%10.5f ,",a[i]); printf("\n"); for (i=1; i<=n+1; i++) printf("----------") printf("\n"); while(n>=1){ for(i=1; i<=n; i++) a[i]=a[i-1]*x+a[i]; for(i=1; i<=n; i++) prontf("%10.5f ,"a[i]); printf("\n"); n=n-1; } return 0; }

専門家に質問してみよう