3 C言語でフーリエ変換

3.1 プログラム

実際のプログラムを以下に示します。このプログラムは、

1. 式(19)を解くために、配列のデータを作成す る部分
2. 連立方程式を計算する部分
3. 連立方程式から、フーリエ級数の関数を計算する部分
4. 関数をグラフにする部分

に分かれています。実際のプログラムは、以下の通りで、次節に注意事項が載っ ているのでそれを見ながら理解してください。

	#include <stdio.h>
	#include <math.h>
	#define MAXN 200

	int gauss_jordan(int n, double a[][MAXN+10], double b[]);
	double func(double x, int n, double b[]);
	void mk_graph(int n, double x[], double x1, double x2,
                             double y[], double y1, double y2);

	/*===========================================================*/
	/*    main function                                          */
	/*===========================================================*/
	main(){
	
	   double a[MAXN+10][MAXN+10], b[MAXN+10];
	   double pi=4.0*atan(1);
	   int n, i, j, k;
	   int nf;
	   double x1, x2, dx, x[1010], y[1010];

	   n=100;
	
	   /* -------- 係数行列と同次項の設定 ---------- */
	   for(i=1 ; i<= n ; i++){
	      for(j=1 ; j<=n ; j++){
	         a[i][j]=cos((double)pi*(i-1)*(j-1)/n);
	      }
	      b[i]=(double)(i-1)/n;
	   }  	

	   /* -------- 連立方程式の計算 ---------------- */
	   if(gauss_jordan(n, a, b) == 0){
	      printf("singular matrix !!!\n");
	      exit(0);
	   };

	   /* -------- 元の関数の計算   ---------------- */
	   nf=1000;
	   x1=-2*pi;
	   x2=2*pi;

	   dx=(x2-x1)/nf;

	   for(i=0 ; i<=nf ; i++){
	      x[i]=x1+dx*i;
	      y[i]=func(x[i], n, b);
	   }

	   /* -------- グラフ作成 ----------------------- */
	   mk_graph(nf,x,x1,x2,y,-0.5,1.5);
   
	}

	/*===========================================================*/
	/*    Gauss-Jordan method                                    */
	/*===========================================================*/
	int gauss_jordan(int n, double a[][MAXN+10], double b[]){

	   int ipv, i, j;
	   double inv_pivot, temp;
	   double big;
	   int pivot_row, row[MAXN+10];

	   for(ipv=1 ; ipv <= n ; ipv++){

	      /* ---- 最大値探索 ---------------------------- */
	      big=0.0;
	      for(i=ipv ; i<=n ; i++){
	         if(fabs(a[i][ipv]) > big){
	            big = fabs(a[i][ipv]);
	            pivot_row = i;
	         }
	      }
	      if(big == 0.0){return(0);}	
	      row[ipv] = pivot_row;

	      /* ---- 行の入れ替え -------------------------- */
	      if(ipv != pivot_row){
	         for(i=1 ; i<=n ; i++){
	            temp = a[ipv][i];
	            a[ipv][i] = a[pivot_row][i];
	            a[pivot_row][i] = temp;
	         }
	         temp = b[ipv];
	         b[ipv] = b[pivot_row];
	         b[pivot_row] = temp;
	      }

	      /* ---- 対角成分=1(ピボット行の処理) ---------- */
	      inv_pivot = 1.0/a[ipv][ipv];
	      a[ipv][ipv]=1.0;
	      for(j=1 ; j <= n ; j++){
	         a[ipv][j] *= inv_pivot;
	      }
	      b[ipv] *= inv_pivot;

	      /* ---- ピボット列=0(ピボット行以外の処理) ---- */
	      for(i=1 ; i<=n ; i++){
	         if(i != ipv){
	            temp = a[i][ipv];
		    a[i][ipv]=0.0;
	            for(j=1 ; j<=n ; j++){
	               a[i][j] -= temp*a[ipv][j];
	            }
	            b[i] -= temp*b[ipv];
	         }
	      }

	   }

	   /* ---- 列の入れ替え(逆行列) -------------------------- */
	   for(j=n ; j>=1 ; j--){
	      if(j != row[j]){
	         for(i=1 ; i<=n ; i++){
	            temp = a[i][j];
	            a[i][j]=a[i][row[j]];
	            a[i][row[j]]=temp;
	         }
	      }
	   }

	   return(1);
	}

	/*===========================================================*/
	/*    function                                               */
	/*===========================================================*/
	double func(double x, int n, double b[]){

	   int i;
	   double y;
	
	   y=0;
	   for(i=1 ; i<=n ; i++){
	      y += b[i]*cos((i-1)*x);
	   }

	   return(y);
	}

	/*==========================================================*/
	/*   make a graph                                           */
	/*==========================================================*/
	void mk_graph(int n, double x[], double x1, double x2,
	                     double y[], double y1, double y2){

	   int i;
	   char *data_file;
	   FILE *out;
	   FILE *gp;

	   /* ------- make a data file ------------ */
	   data_file="gpdata.txt";
	   out = fopen(data_file, "w");
	   for(i=0; i<=n; i++){
	      fprintf(out, "%e\t%e\n", x[i], y[i]);
	   }
	   fclose(out);

	   /* == flowing lines make a graph by using gnuplot == */

	   gp = popen("gnuplot -persist","w");
  	
	   fprintf(gp, "reset\n");
	   fprintf(gp, "set terminal postscript eps color\n");
	   fprintf(gp, "set output \"graph.eps\"\n");
	   fprintf(gp, "set grid\n");

	   /* -------------  set x axis ----------------*/
	   fprintf(gp, "set xtics 1\n");
	   fprintf(gp, "set mxtics 10\n");
	   fprintf(gp, "set xlabel \"%s\"\n", "x");
	   fprintf(gp, "set nologscale x\n");
	   fprintf(gp, "set xrange[%e:%e]\n", x1, x2);

	   /* -------------  set y axis ---------------*/
	   fprintf(gp, "set ytics 1\n");
	   fprintf(gp, "set mytics 10\n");
	   fprintf(gp, "set ylabel \"%s\"\n", "y");
	   fprintf(gp, "set nologscale y\n");
	   fprintf(gp, "set yrange[%e:%e]\n", y1, y2);

	   /* -------------  plat graphs ---------------*/
	   fprintf(gp, "plot \"%s\" using 1:2 with line\n", data_file);

	   fprintf(gp, "set terminal x11\n");
	   fprintf(gp, "replot\n");
						
	   pclose(gp);
	}

3.2 注意点

このプログラムの内容について、以下に注意してください。

• 結果は実数を得たいのですが、整数同士の演算、特に除算は注意が必 要です。FORTRANで学習したように、整数同士の演算はたとえ除算であ ろうが整数になります。それを防ぐために、C言語では強制型変換(キャ スト)というものがあります。(型)式のように書くと、式の内容が示さ れた型に変換されます。例えば、以下のように使います。
  

  	
  	a[i][j]=cos((double)pi*(i-1)*(j-1)/n);
  

• グラフ作成の部分は関数にしました。その関数は、
  

  	
  	mk_graph(データ点数,x座標の配列,x軸左端,x軸右端,y座標の配列,y軸下端,y軸上端);
  

  の形になっています。これを呼び出せば、グラフを書きます。これをコ ピーペーストして使いましょう。

3.3 実行結果

このプログラムを実行すると、図2に示す結果が得られます。

図 2: 実行結果。フーリエ係数は100個

ホームページ: Yamamoto's laboratory
著者: 山本昌志