工学系のための複素数の話  第６巻
「複素数とフーリエ変換」プログラム集

６．１  フーリエ変換

【リスト6-1】関数系の直交を確認するためのプログラム
#include "stdio.h"
#include "Complex.h"
Complex g(int N, int n, int k){return expJ(2.0*PI*n*k/N);}
int main(void){
	int k1, k2,N=10;
	while(scanf("%d,%d",&k1,&k2)!=EOF){
		Complex T=Complex();Complex Z;
		for(int i=0;i<N;i++)T = T + g(N, i,k1)*g(N, i,k2);
		printCom(T);printf("\n");
	}
	getchar();
	return 0;
}
  
【リスト6-2】
#include "stdio.h"
#include "Complex.h"
#define NUMDT 64
Complex X[NUMDT+1], Y1[NUMDT+1], Y2[NUMDT+1];
void setData(Complex X[],int num){
	double DTH=2*PI/num;
	for(int i=0;i<(num+1);i++) X[i]=Complex(cos(DTH*i),0)+0.5;
}
int CSVFlag = true;
void printDT(Complex DT[], int num){
	for(int i=0;i<(num+1);i++){
		if(CSVFlag)	printf("\n %4d , %lf, %lf",i, DT[i].r, DT[i].i);
		else{ printf("\n %4d :",i); printCom(DT[i]);}
	}
}
void DFT(Complex X[], Complex Y1[], int num){
	double DN = -PI*2/num;
	for(int k=0;k<num;k++){
		Y1[k] = Complex();
		for(int N=0; N<num; N++)
			Y1[k] = Y1[k] + X[N]*expJ(DN*N*k);
	}
	Y1[num]=Y1[0];
}
void IDFT(Complex Y1[], Complex Y2[], int num){
	double DN = PI*2/num;
	for(int k=0;k<num;k++){
		Y2[k] = Complex();
		for(int N=0; N<num; N++)
			Y2[k] = Y2[k] + Y1[N]*expJ(DN*N*k);
	}
	for(int k=0;k<num;k++)Y2[k]=Y2[k]/num;
	
	Y2[num]=Y2[0];
}
int main(void){
	setData(X, NUMDT);
	printf("\n\n ■ Data");
	printDT(X, NUMDT);
	DFT(X, Y1, NUMDT);
	printf("\n\n ■ DFT Result");
	printDT(Y1, NUMDT);
	IDFT(Y1,Y2, NUMDT);
	printf("\n\n ■ IDFT Result");
	printDT(Y2, NUMDT);
	
	getchar();
	return 0;
}

【DFTとIDFTの処理共通化】
 void DFT_IDFT(Complex Y1[], Complex Y2[], int num){
 	double DN = -PI*2/num;int n=abs(num);
 	for(int k=0;k<n;k++){
 		Y2[k] = Complex();
 		for(int N=0; N<n; N++)
 			Y2[k] = Y2[k] + Y1[N]*expJ(DN*N*k);
 	}
 	if(num<0) for(int k=0;k<n;k++)Y2[k]=Y2[k]/num;
 	Y2[n]=Y2[0];
 } 

６．２  高速フーリエ変換

【リスト6-3】
#include "stdio.h"
#include "Complex.h"
#define swap(type, x,y) do{type t; t=x; x=y; y=t;}while(0)
#define NUMDT 64
Complex X[NUMDT+2]; //, Y1[NUMDT+1], Y2[NUMDT+1];
void setData(Complex X[],int num){
	double DTH=2*PI/num;
	for(int i=0;i<(num+1);i++) X[i]=Complex(cos(DTH*i),0)+0.5;
}
int CSVFlag = true;
void printDT(Complex DT[], int num){
	for(int i=0;i<(num+1);i++){
		if(CSVFlag)	printf("\n %4d , %lf, %lf",i, DT[i].r, DT[i].i);
		else{ printf("\n %4d :",i); printCom(DT[i]);}
	}
}
void FFT1(Complex x[], int num){
	Complex tmp;double argV;
	int n,j,jnh,jp,jx,k,ne,NH,NH2;
	n=abs(num);	argV=-2.0*PI/num;
	for(j=0;j<n;j++){
		for(jp=j, jx=0, k=1; k<n; k<<=1,jp>>=1 )jx=(jx<<1)|(jp & 1);
		if(j<jx)swap(Complex,x[j],x[jx]);
	}
	for(NH=2, ne=n/2; ne>=1;ne>>=1, NH<<=1)
		for(NH2=NH/2,jp=0;jp<n; jp+=NH)
			for(j=jp, jx=0;j<jp+NH2;j++,jx+=ne){
				jnh=j+NH2;
				tmp=x[jnh]*expJ(argV*jx);
				x[jnh]=x[j]-tmp;
				x[j] =x[j]+tmp;
			}
	if(num<0) for(j=0;j<n;j++) x[j]= X[j]/n;
	x[n]=x[0];
}
int main(void){
	setData(X, NUMDT);
	printf("\n\n ■ Data");printDT(X, NUMDT);
	FFT1(X, NUMDT);
	printf("\n\n ■ DFT Result");printDT(X, NUMDT);
	FFT1(X,-NUMDT);
	printf("\n\n ■ IDFT Result");printDT(X, NUMDT);
	getchar();
	return 0;
}
【リスト6-4】周波数間引きFFTの関数
・
・（前略）【リスト6-3】を参照
・
void FFT2(Complex x[], int num){
	Complex tmp;double argV;
	int n,j,jnh,jp,jx,k,ne,NH,NH2;
	n=abs(num); argV=-2.0*PI/num;
	for(NH=n/2, ne=1; ne<n;ne<<=1, NH>>=1)
		for(NH2=NH*2,jp=0;jp<n; jp+=NH2)
			for(j=jp, jx=0;j<jp+NH;j++,jx+=ne){
				jnh=j+NH; tmp=x[j];
				x[j] =tmp+x[jnh];
				x[jnh]=(tmp-x[jnh])*expJ(argV*jx);
			}
	for(j=0;j<n;j++){
		for(jp=j, jx=0, k=1; k<n; k<<=1,jp>>=1 )jx=(jx<<1)|(jp & 1);
		if(j<jx)swap(Complex,x[j],x[jx]);
	}
	if(num<0) for(j=0;j<n;j++) x[j]= X[j]/n;
	x[n]=x[0];
}
int main(void){
	setData(X, NUMDT);
	printf("\n\n ■ Data");printDT(X, NUMDT);
	FFT2(X, NUMDT);
	printf("\n\n ■ DFT Result");printDT(X, NUMDT);
	FFT2(X,-NUMDT);
	printf("\n\n ■ IDFT Result");printDT(X, NUMDT);
	getchar();
	return 0;
}

【リスト6-5】  回転因子とビット逆順の値をあらかじめ計算
Complex wTab[NUMDT/2];	//回転因子計算用テーブル
int bTab[NUMDT]; 	//ビット逆順交換用テーブル
static int n_old=0;
int chRAD2(int num){
	int pow2=1,j;	//要素数が2のべき乗かどうかを判定
	for(j=0;j<32;j++) if(num==pow2) break; else pow2<<=1;
	if(j<32)return 0;
	printf("\n データ数が2のべき乗になっていません");
	return 1;
}
int FFT(Complex x[], int n_FFT){
	Complex xtmp; int j,jnh,k,jx,num=abs(n_FFT), m, NH, NH2;
	if(num !=n_old){   //1回目の呼び出しのときのみテーブル作成
		if(chRAD2(num))return 1;
		double argV=-2.0*PI/num; //回転因子計算用テーブル作成
		for(j=0;j<num/2;j++) wTab[j]=expJ(argV*j);
		bTab[0]=0; //ビット逆順交換用テーブル作成
		for(m=1, NH=num/2; m<num;m<<=1,NH>>=1)
			for(k=0;k<m;k++)bTab[k+m]=bTab[k]+NH;
		n_old=num;
	}
	//FFT計算開始
	for(m=1,NH=num/2; m<num; m<<=1,NH>>=1)
		for(NH2=NH<<1, k=0; k<num; k += NH2)
			for(j=k, jx=0; j<k+NH; j++,jx += m){
				jnh=j + NH; xtmp = x[j]; x[j] = xtmp+x[jnh];
				x[jnh]=(xtmp-x[jnh])*wTab[jx];
			}
	//ビット逆順入れ換え
	for(j=0; j<num;j++) if(j<bTab[j]) swap(Complex, x[j],x[bTab[j]]);
	if(n_FFT<0) for(j=0;j<num; j++)	x[j]=x[j]/num;//逆FFT
	return 0;
}

６．３  ２次元フーリエ変換

【リスト6-6】２次元DFT
#include "stdio.h"
#include "Complex.h"
#define swap(type, x,y) do{type t; t=x; x=y; y=t;}while(0)
#define GHDFT 64
#define GWDFT 64
Complex X[GWDFT+1][GHDFT+1],tmDFT[GWDFT+1][GHDFT+1]; 
Complex Y1[GWDFT+1][GHDFT+1],Y2[GWDFT+1][GHDFT+1]; 
void setData(Complex X[][GHDFT+1],int num){
	double DTH=2*PI/num;
	for(int i=0;i<(num+1);i++) 
		for(int j=0;j<(num+1);j++)
			if((i>(num/4) && i<(3*num/4))&&
				(j>(num/4) && j<(3*num/4)))
				X[i][j]=Complex(1,0);
			else X[i][j]=Complex(0,0);
}
int CSVFlag = true;
void printDT(Complex DT[][GHDFT+1], int num, int image){
	for(int i=0;i<(num+1);i++){
		printf("\n %4d", i);
		for(int j=0;j<(num+1);j++){
			double D=DT[i][j].r; if(image) D=DT[i][j].i;
			printf(", %lf",D);
		}
	}
}
void fprintDT(FILE *fn, Complex DT[][GHDFT+1], int num, int image){
	for(int i=0;i<(num+1);i++){
		fprintf(fn,"\n %4d", i);
		for(int j=0;j<(num+1);j++){
			double D=DT[i][j].r; if(image) D=DT[i][j].i;
			fprintf(fn,", %lf",D);
		}
	}
}
void DFT(Complex beforDFT[][GHDFT+1], Complex afterDFT[][GHDFT+1]){
	double ARG=-2.0*PI/GHDFT;
	for(int k=0;k<=GWDFT;k++)for(int N=0;N<=GWDFT;N++){//中央への移動
		double DD=((k+N)%2)!=0 ? -1.0: 1.0;
		afterDFT[k][N]=beforDFT[k][N]*DD;
	}
	for(int i=0;i<GWDFT;i++){//Y方向変換
		for(int k=0;k<GHDFT;k++){
			tmDFT[i][k]=Complex(0,0);
			for(int N=0;N<GHDFT;N++){
				tmDFT[i][k]=tmDFT[i][k]+afterDFT[i][N]*expJ(ARG*N*k);
			}
		}
		tmDFT[i][GHDFT]=tmDFT[i][0];
	}
	tmDFT[GWDFT][GHDFT]=tmDFT[0][0];
	ARG=-2.0*PI/GWDFT;
	for(int i=0;i<GHDFT;i++){//X方向変換
		for(int k=0;k<GWDFT;k++){
			afterDFT[k][i]=Complex(0,0);
			for(int N=0;N<GWDFT;N++){
				afterDFT[k][i]=afterDFT[k][i]+tmDFT[N][i]*expJ(ARG*N*k);
			}
		}
		tmDFT[i][GHDFT]=tmDFT[i][0];
	}
}
//逆DFT
void IDFT(Complex afterDFT[][GHDFT+1],  Complex reverceDFT[][GHDFT+1]){
	double ARG=2.0*PI/GHDFT; //Y方向変換
	for(int i=0;i<GWDFT;i++){
		for(int k=0;k<GHDFT;k++){
			reverceDFT[i][k]=Complex(0,0);
			for(int N=0;N<GHDFT;N++){
				reverceDFT[i][k]=reverceDFT[i][k]+afterDFT[i][N]*expJ(ARG*N*k);
			}
		}
		reverceDFT[i][GHDFT]=reverceDFT[i][0];
	}
	reverceDFT[GWDFT][GHDFT]=reverceDFT[0][0];
	for(int i=0;i<=GWDFT;i++){
		for(int k=0;k<=GHDFT;k++){
			reverceDFT[i][k]=reverceDFT[i][k]/GHDFT;
		}
	}
	ARG=2.0*PI/GWDFT; //X方向変換
	for(int i=0;i<GHDFT;i++){ 
		for(int k=0;k<GWDFT;k++){
			tmDFT[k][i]=Complex(0,0);
			for(int N=0;N<GWDFT;N++){
				tmDFT[k][i]=tmDFT[k][i]+reverceDFT[N][i]* expJ(ARG*N*k);
			}
		}
		tmDFT[GWDFT][i]=tmDFT[0][i];
	}
	for(int i=0;i<=GWDFT;i++){
		for(int k=0;k<=GHDFT;k++){
			tmDFT[i][k]=tmDFT[i][k]/GHDFT;
		}
	}
	for(int k=0;k<=GWDFT;k++)for(int N=0;N<=GWDFT;N++){//周辺への移動
		double DD=((k+N)%2)!=0 ? -1.0: 1.0;
		reverceDFT[k][N]=tmDFT[k][N]*DD;
	}
}
int main(void){
	FILE *fn=fopen("d:\\testdt.csv","w");
	setData(X,GWDFT );
	fprintf(fn,"\n\n ■ Data");	fprintDT(fn, X, GWDFT, false);
	DFT(X,Y1);
	fprintf(fn,"\n\n ■ DFT Result Re");fprintDT(fn,Y1, GWDFT,false);
	fprintf(fn,"\n\n ■ DFT Result Im");fprintDT(fn,Y1, GWDFT,true);
	IDFT(Y1, Y2);
	fprintf(fn,"\n\n ■ IDFT Result");fprintDT(fn,Y2, GWDFT,false);
	fclose(fn);	getchar();
	return 0;
}

【リスト6-7】  2次元FFT
void fft_2d(Complex beforDFT[][GHDFT+1],Complex afterDFT[][GHDFT+1],
					int n_fft){
	int n1, n2,num=labs(n_fft);  Complex u[GWDFT+1];
	for(n1=0; n1<num; n1++){
		for(n2=0; n2<num; n2++){
			double DD=((n1+n2) & 1) ? -1.0: 1.0;
			afterDFT[n1][n2] = beforDFT[n1][n2]*DD;
		}
		if(FFT(afterDFT[n1], n_fft))return;
	} 
	for(n2=0; n2<num; n2++){
		for(n1=0; n1<num; n1++) u[n1]=afterDFT[n1][n2];
		FFT(u, n_fft);
		for(n1=0; n1<num; n1++) afterDFT[n1][n2]=u[n1];
	}
} 
  	
【リスト6-8】  高速化2次元FFT
void fft_2d(Complex beforDFT[][GHDFT+1],Complex afterDFT[][GHDFT+1],
						 int n_fft){
	int n1, n2,num=labs(n_fft),num1=num+1; Complex u[GWDFT+1];
	memmove(afterDFT, beforDFT,sizeof(Complex)*num1*num1);
	for(n1=0; n1<num; n1++) if(FFT(afterDFT[n1], n_fft))return;
	for(n2=0; n2<num; n2++){
		for(n1=0; n1<num; n1++) u[n1]=afterDFT[n1][n2];
		FFT(u, n_fft);
		for(n1=0; n1<num; n1++) afterDFT[n1][n2]=u[n1];
	}
}

【リスト6-9】  スペクトル位置の移動(GWDFT=64, GHDFT=64)
void setDspSpec(Complex TM[][GHDFT+1], Complex DT[][GHDFT+1]){
	int w2=GWDFT/2, h2=GHDFT/2;
	for(int x=0; x<GWDFT; x++)for(int y=0;y<GHDFT;y++){
		int x2=x+w2, y2=y+h2;
		if(x2>GWDFT) x2= x2-GWDFT; if(y2>GHDFT) y2= y2-GHDFT;
		TM[x2][y2]=DT[x][y];
	}
} 

【リスト6-10】  2次元高速逆フーリエ変換(GWDFT=64, GHDFT=64)
void rev_fft_2d(Complex afterDFT[][GHDFT+1],
						  Complex reverceDFT[][GHDFT+1], int n_fft){
	int num=labs(n_fft);
	fft_2d(afterDFT,reverceDFT, -num);

	for(int k=0;k<=num;k++)for(int N=0, N2=num-1;N<num/2;N++,N2--)
			swap(Complex, reverceDFT[k][N],reverceDFT[k][N2]);
	for(int k=0;k<=num;k++)for(int N=0, N2=num-1;N<GWDFT/2;N++,N2--)
				swap(Complex, reverceDFT[N][k],reverceDFT[N2][k]);
	}