C言語による数値計算プログラミング演習

Size: px

Start display at page:

Download "C言語による数値計算プログラミング演習"

とらふみひがき
6 years ago
Views:

1 11. 離散フーリエ変換時間領域における連続関数 x(t) は, 周波数領域の連続関数 (f) を介して, フーリエ変換 (Fourier transform) : i2π ft ( f) = xte ( ) dt 逆フーリエ変換 (inverse Fourier transform): 2 x() t = ( f) e i π ft df と展開される x(t) は区間 [,T] 以外ではとすると, フーリエ変換における無限区間 [, ] での積分は有限区間 [,T] での積分となるこの区間を等分 ( 小区間幅 :T s =T/, 等分点 :t j =jt s (j=,1,,)) して複合台形則により数値積分を行い, さらに周波数 f k =k/t (k=,1,,-1) において評価すると, ( ) 1 i2 π kj/ 2 k = s j + s j= ( f ) T x( t ) e O T を得るこの式を T s で割って y(f k )=(f k )/T s とおいたものが, 離散フーリエ変換 (discrete Fourier transform, DFT) である x(t j )=x j, y(f k )=y k と記せば,DFT は, 時間領域における時間間隔 T s ( サンプリングタイムとよぶ ) の点のデータ Tを周波数領域における周波数間隔 f 1 =1/T( 基本周波数とよぶ ) の点のデータ x = [ x, x1, L, x 1] y= [ y, y1,l, y 1] T i2 π / w = e に変換するものであり, 回転因子を用いて, 以下のように書ける 1 kj yk = x jw ( k =,1, L, 1) j = 行列ベクトル表示を行うと, L w w L w 2 4 2( 1) y = Fx, F = 1 w w L w M M M M ( 1) 2( 1) ( 1)( 1) 1 w w L w となる逆離散フーリエ変換 (inverse discrete Fourier transform, IDFT) は,y から x を求めるプロセスであり, 次のようになる 1 x = y w ( j =,1, L, 1) 1 jk j k k = これは, 行列ベクトル表示により次のように書ける L 1 1 w w L w x = F F = w w w M M M M 1 w w L w 1 2 ( 1) ( 1) y, 1 L 周期性 : k k± j w = w ( j=,1, L ) ( 1) 2( 1) ( 1)( 1) 回転因子 w は 1 の乗根であるその k 乗あるいは-k 乗 (k=,1,,-1) もすべて 1 の乗根であり, これらはそれぞれ複素平面上の単位円の角度干 2πk/(k=,1,,-1) における点である回転因子の演算には, 次の性質およびが偶数のときには原点対称性 : k+ /2 w = w がある =8 のときの回転因子の例を図示する高速フーリエ変換 (fast Fourier transform, FFT) とは, 上記の変換を効率よく計算するアルゴリズムであり, 時間間引き型 FFT と周波数間引き型 FFT がある 1 k

2 11.1 時間間引き型 FFT 時間間引き型 FFT とは,DFT の右辺 ( 時間領域におけるデータ ) を以下のように x k が偶数と奇数のグループに並べ替えて, データ点数の DFT をデータ点数 /2 の DFT に分解していき, 演算量を減らすものである 1 /2 1 /2 1 kj (2 r ) k (2r+ 1) k () k (1) k j 2r 2r+ 1 k k j= r= r= y = x w = x w + x w = + w ( k =,1, L, 1) /2 1 () rk ' () () k' = x2 rw /2, k' + /2 = k' r = /2 1 (1) rk ' (1) (1) k' = x2r+ 1 w /2, k' + /2 = k' r = データ総数が =2 m の場合,p 段目ではデータ点数 p =2 p の DFT グループが n g =/ p =2 m-p 個ある p 段の最初において, 各グループ毎に前半には偶数列要素, 後半には奇数列要素を順に配置すると,p 段第 i グループ (i=,1,,n g -1) の前半 ( 偶数列要素 ) は p-1 段の第 2i グループ要素, 後半 ( 奇数列要素 ) は p-1 段の第 1 グループ要素に相当する p 段 i グループの第 k 要素を新たに p i,k (k=,1, p -1) と記すと,w k+p/2 p =-w k p より p p 1 j p 1 i, j = 2 i, j + w 2 1, (,, 1) p i+ j i = L ng p p 1 j p 1 ( j,, /2 1) i, j+ p/2 = 2 i, j w = L p 1, j p と分解されることになるこれを p=1 から始め,p=m まで繰り返せばよい =2 3 =8 の場合の例を示す最後の第 3 段では,8 点の DFT グループ ( 第グループ ) が 1 組存在し, この順番で周波数領域における値 y k となっている i, = 2 i, + w8 1, i,1= 2,1 i + w8 1, i, 2= 2 i, 2+ w8 1, i, 3= 2 i, 3+ w8 1, ( i = ) i, 4= 2 i, w8 1, i, 5= 2 i, 1 w8 1, i, 6= 2 i, 2 w8 1, i, 7= 2 i, 3 w8 1, 3 ここで右辺の ( 2,, 2,1, 2,2, 2,3 ) は偶数項 (x,x 2,x 4,x 6 ) をデータに持つ第 2 段の第 DFT グループ,( 2 1,, 2 1,1, 2 1,2, 2 1,3 ) は奇数項 (x 1,x 3,x 5,x 7 ) をデータに持つ第 2 段の第 1DFT グループであるこのように, 第 2 段では,4 点の DFT グループが 2 組存在し, さらに次のように分解される i, = 2 i, + w4 1, i,1= 2,1 i + w4 1,1 (,1) i = i, 2= 2 i, w4 1, i, 3= 2 i, 1 w4 1, 1 右辺の ( 1,, 1,1 ), (1 1,, 1 1,1 ), (1 2,, 1 2,1 ), (1 3,, 1 3,1 ) はそれぞれの偶数項 (x,x 4 ), 奇数項 (x 2,x 6 ), 偶数項 (x 1,x 5 ), 奇数項 (x 3,x 7 ) をデータに持つ第 1 段の第,1,2,3DFT グループであるつまり第 1 段では,2 点の DFT グループが 4 組存在し, さらに次のように分解される 1 i, = 2 i, + w2 1, (,1,2,3) 1 i = i,1= 2, i w2 1, 実際の計算は第 1 段から始めて第 2 段を行い, 第 3 段で終了するデータ処理の流れは以下のように書ける 2

3 第 1 段の最初において与えるべきデータ x k の順番 ( 添え字番号 k の順番 ) は,,1,,7 を表す 2 進数 () 2, (1) 2, (1) 2, (11) 2, (1) 2, (11) 2, (11) 2, (111) 2 のビットを逆転させたもの ( 以下 ) となっている () 2, (1) 2, (1) 2, (11) 2, (1) 2, (11) 2, (11) 2, (111) 2 = =4 =2 =6 =1 =5 =3 =7 アルゴリズム DFT のアルゴリズムは以下のとおりただし,Y, k, w p j は複素数である 1) ビット逆転により,=2 m 点の入力データ x i ( 対応する時間は t i =it s ) の順序を並べ替え, その結果を変数, 1,, -1 にセットする for i= to -1 データ数 =2 m だけ繰り返す bit := i i を 2 進数で表したときのビットを逆転して変数 bitr にセット bitr := ; C 言語の左右シフト演算子 <<, >>, 論理積演算子 &, 論理和演算子を用いる for k=1 to m bitr := bitr << 1 bitrを左に1ビットだけシフト bitr := bitr (bit & 1) bitの1ビット目をbitrの1ビット目にセット bit := bit >> 1 bitを右に1ビットだけシフト if i<bitr then i < bitr のときのみ,i 番目とbitr 番目のデータを入れ換え Y := i 注 : 既に交換済みのものを再度交換すると元に戻るので i := bitr bitr := Y 2) 偶数列と奇数列への分解の統合を m 段繰り返す次に示すのはアルゴリズムの原型であり,w p j の計算を効率よく行うべく更に改良できる for p=1 to m 第 p 段 :p=1,,m p := 2 p p: 各グループのデータ点数 g := /p ng: グループ数 for i= to ng-1 for j= to p/2-1 k := p i + j l := k+p/2 Y := w j p l l := k -Y k := k +Y この結果, k が DFT の y k ( 対応する周波数は f k =k/(t s )) となるなお,IDFT も全く同じアルゴリズムにより計算できる DFT との相違は回転因子 w を w -1 に置き換え, 最後に (1/) を乗じることのみであるプログラム C 言語のコンパイラによっては, 複素数を扱えないものもあるここでは, 構造体により complex 型を定義する 1 2 /* FFT (Fast Fourier Transform) program */ /* with Complex type supplied */ 3

4 #include <stdio.h> #include <math.h> #define PI #define MA 8192 /* MA = 2^13 =8192 */ typedef struct /* 複素数を扱うための complex 型構造体を定義 */ double re; double im; complex; complex comp_set(double x, double y) /* 複素数の実部と虚部のセット */ complex z; z.re = x; z.im = y; return(z); complex comp_add(complex x, complex y) /* 複素数の加算 */ complex z; z.re = x.re + y.re; z.im = x.im + y.im; return(z); complex comp_sub(complex x, complex y) /* 複素数の減算 */ complex z; z.re = x.re - y.re; z.im = x.im - y.im; return(z); complex comp_mul(complex x, complex y) /* 複素数の乗算 */ complex z; z.re = x.re * y.re - x.im * y.im; z.im = x.re * y.im + x.im * y.re; return(z); complex comp_div(complex x, complex y) /* 複素数の割算 */ complex z; 4

5 double c; c = y.re * y.re + y.im * y.im; z.re = ( x.re * y.re + x.im * y.im)/c; z.im = (-x.re * y.im + x.im * y.re)/c; return(z); void fft_time(complex *x, int n, int M, int idft) complex w,wj,y; float dth; int i,j,k,l,p,bit,bitr,sign,np,np2,ng; for(i=;i<n;i++) /* bit 逆転 : i ==> bitr */ bit =i; bitr=; for(k=;k<m;k++) bitr = bitr << 1; /* bitrを左に 1ビットだけシフト */ bitr = bitr (bit & 1); /* bitの1ビット目を bitrの1ビット目にセット */ bit = bit >> 1; /* bitを右に 1ビットだけシフト */ if(i<bitr) /* i < bitr のときのみ, データを入れ換える */ y = x[i]; x[i] = x[bitr]; x[bitr]= y; if(!idft) /* DFT */ sign=-1; else /* IDFT */ sign= 1; np=1; for(p=1;p<=m;p++) np *= 2; /* p 段目のデータ点数 :np = 2^p */ np2 = np/2; ng =n/np; /* p 段目のグループ数 :ng = n/np */ dth =PI/np2; /* dtheta = 2PI/np */ w = comp_set( cos(dth), sign*sin(dth) ); wj = comp_set( 1.e,.e ); for(j=;j<np2;j++) for(i=;i<ng;i++) k =np*i+j; 5

6 l =k+np2; y =comp_mul( wj, x[l] ); x[l]=comp_sub( x[k], y ); x[k]=comp_add( x[k], y ); wj = comp_mul( wj, w ); if(idft) /* IDFT */ for(i=;i<n;i++) x[i].re /= (float)n; x[i].im /= (float)n; void main(void) complex x[ma]; float time_smpl,sp,time[ma],freq[ma],single_re,single_im; int n,m,i,spectrum,idft=; char fname[3][3]; /* ファイル名用領域 */ FILE *fpin, *fpout; /* ファイルポインターの宣言 */ printf("fast Fourier Transform\n"); printf("* umber of data n=2^m : n, M? "); /* データ点数 :n = 2^M */ scanf("%d%d", &n, &M); /* ファイル名の入力 */ printf("input file name for data ==> "); scanf("%s",fname[]); /* idft= のときは後者 (DFT) を, それ以外のときは前者 (IDFT) を出力 */ printf("output file name for %s complex data ==> ", idft?"idft":"dft"); scanf("%s",fname[1]); printf("output file name for %s spectrum data ==> ", idft?"idft":"dft"); scanf("%s",fname[2]); /* ファイルオープンとエラー処理 */ if( (fpin=fopen(fname[],"r"))==ull ) printf(" ファイルをオープンできません \n"); exit(1); for(i=;i<n;i++) if(!idft) /* data for DFT */ fscanf(fpin, "%g%g",&time[i], &single_re); single_im=.; else /* data for IDFT */ fscanf(fpin, "%g%g%g",&time[i], &single_re, &single_im); x[i]=comp_set(single_re, single_im); 6

7 fclose(fpin); printf("\n* input data were read from file: %s\n",fname[]); /* set of sampling time */ if(!idft) /* DFT: time[]... time */ time_smpl = time[1]-time[]; else /* IDFT: time[]... frequency */ time_smpl = 1./((time[1]-time[])*n); fft_time(x, n, M, idft); /* FFT */ printf("!!! %s finished!!!\n", idft?"idft":"dft"); if(!idft) /* DFT: set frequency */ for(i=;i<n;i++) freq[i]=(float)i/(time_smpl * n); else /* IDFT: set time */ for(i=;i<n;i++) freq[i]=(float)i*time_smpl; /* 変換結果を複素数でファイル出力 */ if( (fpout=fopen(fname[1],"w"))==ull ) printf(" ファイルをオープンできません \n"); exit(1); for(i=;i<n;i++) fprintf(fpout,"%13.6g\t%13.6g\t%13.6g\n",freq[i],x[i].re,x[i].im); fclose(fpout); printf("* complex output were written to file: %s\n",fname[1]); /* 出力スペクトルの種類の選択 */ printf("\ninput spectrum o.:.phase, 1.amplitude, 2.power spectrum? "); scanf("%d",&spectrum); /* 変換結果を各種スペクトルでファイル出力 */ if( (fpout=fopen(fname[2],"w"))==ull ) printf(" ファイルをオープンできません \n"); exit(1); for(i=;i<n;i++) if(spectrum==) sp = atan2( x[i].im, x[i].re ); else if(spectrum==1) sp = sqrt( x[i].re*x[i].re + x[i].im*x[i].im ); else if(spectrum==2) sp = x[i].re*x[i].re + x[i].im*x[i].im; fprintf(fpout,"%13.6g\t%13.6g\n",freq[i],sp); fclose(fpout); printf("* spectrum output were written to file: %s\n",fname[2]); 7

8 実行例時間 t[s] における変位 x[mm] が sin 波による合成 ( 下式 ) により与えられるとする x( t) = a + a sin(2 πν t) + a sin(2 πν t) + a sin(2 πν t) + a sin(2 πν t) ここで, 各振幅と周波数は以下のように設定する a = 2 a1 = 5 ν1 =.8 a2 = 2 ν 2 = 1.35 a3 = 1 ν 3 = 3 a = 1.5 ν = サンプリングタイム T s を.1[s] として =2 1 =124 個のデータを ( エクセルなどにより ) 計算してテストデータとし, ファイル test1.txt に格納するテストデータの時間と変位を以下に図示するファイル test1.txt の中身 ( 時間と変位 ) : : : : 変位 x [mm] 時間 t [sec.] DFT の結果を以下に示す実行開始 Fast Fourier Transform * umber of data n=2^m : n, M? 入力データ点数 124は2の1 乗 input file name for data ==> test1.txt 入力データファイル名 output file name for DFT complex data ==> test1_complex.txt 出力データ ( 複素数 ) ファイル名 output file name for DFT spectrum data ==> test1_spectrum.txt 出力データ ( スペクトル ) ファイル名 * input data were read from file: test1.txt!!! DFT finished!!! * complex output were written to file: test1_complex.txt input spectrum o.:.phase, 1.amplitude, 2.power spectrum? 1 出力スペクトルの選択 * spectrum output were written to file: test1_spectrum.txt おしまい

9 出力スペクトルファイル (test1_spectrum.txt) から振幅スペクトルをエクセルに読み込み, 図示するこの変換におけるサンプリングタイムは T s =.1[s], サンプリング点数は =124, サンプリング区間 ( サンプリングを行った全時間 ) は T=T s =1.24[s] であるしたがって基本周波数は f 1 =1/T=.97656[Hz] であり, ファイルのデータをチェックすると出力周波数は基本周波数の倍数 kf 1 (k=,1, -1) であることが確認されるサンプリング周波数は f s =1/T s =1[Hz] であり, 結果の周波数の範囲は [,f s ) となっていることが確認される振幅スペクトル図には, 中央 (f=f s /2) に関して対称となる特性が現れているサンプリング定理によると, ある波形に含まれる最大周波数を f max とすると, その波形を復元するには 2f max よりも高い周波数 f s でサンプリングする必要があるしたがって計算された周波数の範囲 [,f s ) のうち意味を持つのは,[, f s /2) の範囲であるグラフには, 周波数.8, 1.35, 3, 1 近傍にピークがあり, 入力波が適正にフーリエ変換されたことがわかる振幅スペクトルファイル test1_spectrum.txt の中身 ( 周波数と振幅 ) : : : : : : 周波数 Hz] 周波数 [Hz] 周波数 [Hz] 9

10 確認のため, 更にプログラムの 127 行における変数 idft の初期設定を idft=1 とし,IDFT の計算を行う入力データには先ほどの DFT で得られた結果 ( 複素数 ) のファイルを用いると,DFT での最初の入力波形が復元されるはずである実行開始 Fast Fourier Transform * umber of data n=2^m : n, M? 入力データ点数 124は2の1 乗 input file name for data ==> test1_complex.txt 入力データファイル名 output file name for IDFT complex data ==> test1_idftcomplex.txt 出力データ ( 複素数 ) ファイル名 output file name for IDFT spectrum data ==> test1_idftspectrum.txt 出力データ ( スペクトル ) ファイル名 * input data were read from file: test1_complex.txt!!! IDFT finished!!! * complex output were written to file: test1_idftcomplex.txt input spectrum o.:.phase, 1.amplitude, 2.power spectrum? 1 * spectrum output were written to file: test1_idftspectrum.txt おしまい DFT の結果データ (test1_complex.txt) を入力して IDFT を行った結果, 複素数の実部には最初の入力波形がほぼ復元されており, 複素数の虚部は全てほぼ零に近い値となっており, 逆フーリエ変換が適正に行われたことがわかるファイル test1_idftcomplex.txt の中身 ( 時間と変位 : 複素数の実部と虚部 ) E E E E E E E E E E E E E-6 : : : : : : : : : E E E E-7 変位 [mm] 実部虚部時間 [sec.] 参考文献葉子 : 数値計算の基礎解法と誤差, コロナ社 (27) Heath, Michael T.: Scientific Computing, An Introductory Survey, McGraw-Hill(22) 1

DVIOUT

DVIOUT 第章離散フーリエ変換離散フーリエ変換これまで私たちは連続関数に対するフーリエ変換およびフーリエ積分 ( 逆フーリエ変換 ) について学んできましたこの節ではフーリエ変換を離散化した離散フーリエ変換について学びましょう自然現象 ( 音声 ) などを観測して得られる波 ( 信号値 ; 観測値 ) は通常電気信号による連続的な波として観測機器から出力されますしかしながらコンピュータはこの様な連続的な波を直接扱うことができないため