Subsections
3 C言語でフーリエ変換
3.1 プログラム
実際のプログラムを以下に示します．このプログラムは，- 式(19)を解くために，配列のデータを作成す る部分
- 連立方程式を計算する部分
- 連立方程式から，フーリエ級数の関数を計算する部分
- 関数をグラフにする部分
#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に示す結果が得られます．
ホームページ: Yamamoto's laboratory
著者: 山本昌志 Yamamoto Masashi
2005-11-25