乱流工学（スペクトル解析）
機械工学科 飯田 明由
要約

ランダムデータの解析の基礎となるのがフーリエ解析である．フーリエ解析は任意の周期関数を三角関数の級数として表す方法であり，微分方程式の解としてだけではなく，不規則現象の解析である確率論的な解析に拡張されている．工学的には乱流理論への拡張が最も成功した例であり，乱流の中の渦塊の運動を速度変動のスペクトルとして表した．前回までの講義ではKolmogorovの-5/3乗則などスペクトル解析によって得られた研究成果を解説してきたが，この章ではスペクトル解析を行うためのフーリエ解析，特に実際の計測・解析で利用されるFFT（Fast FourierTransform）について解説する．

1. フーリエ級数
19世紀の数学者Fourier(1764−1830)は周期関数を三角関数の級数として表す方法を考案した．
今日フーリエ級数と呼ばれるこの方法は発表当時は学会から正式に認められたものではなかった
が，現在では微分方程式の解法手段としてだけなく，信号波形，特にランダム信号の解析に広く利用されている．（Fourierは熱伝導問題を解くためにフーリエ級数を利用した．）
関数x(t)は区間(-T/2, T/2)を基本とする周期Tの周期関数とする．
1.51x(t)x(t)0.50-0.5-1-1.5-1-0.50t/T

図1 三角系列関数x(t)4
− 2 − tT
 4tx(t ) =  T
 2− 4 t
TT 3T<t <− 
−44

T
 T
− ≤ t ≤ − 4 4TT <t < 4
4-1-0.51この関数を整数個の波を含む正弦関数の過重和で表現するとする．x(t)は奇関数なので，余弦関数の項は考えなくとも良い．
x(t ) = b1 sin2πt4πt6πt2nπt+ b 2 sin+ b3 sin+ K + bn sin+KTTTT

両辺にsin(2nπt/T)を掛けて，積分すると∫T /2x(t ) sin−T / 2T /2T /22nπt2πt2nπt4πt2nπtsinsin
⋅ sin
⋅ sin
dt = b1
dt + b 2dt
TTTT
T−T / 2−T / 2∫∫

+ K + bn2nπt  2 sin dt + KT −T / 2 ∫

T /2

ここで，

∫T /2sinsin−T / 2=02mπt2nπt dt  T⋅ sinTT = 2( m ≠ n)( m = n)で あ る か ら ， m=n の 時 以 外 は 積 分 値 は 0 と な る ． こ の よ う な 関 係 を 持 つ 関 数 列 を 直 交 関 数(orthogonal functions)と呼ぶ．m=n のとき積分値が1 となる場合，正規直交関数(orthonormalfunctions)という．

∫T /2x(t ) sin−T / 2T /2
 T /4 4 2nπt2nπt4 2nπt dt = 2 dt +dt  t  sin 2 − t  sin
TTT TT /4  0  T ∫=∫2T nπ / 22T nπ2Tsinζζdζ−ζ sin ζdζ +22(nπ) 0(nπ) nπ / 2nπ∫∫ 4T(−1) n + 1= (nπ)2
 = 0nπ∫ ζ sin ζdζnπ / 2(n = 1,3,5, K)(n = 2,4,6, K)

したがって，bn =2T

∫

T /2

−T / 2x(t ) sin

2nπtdtT 8T(−1) n + 1
=2
 (nπ)
 = 0(n = 1,3,5, K)(n = 2,4,6, K)

以上の結果から先に与えた周期関数x(t)はx(t ) =14πt
10πt 16πt 1
8  2πt 1− 2 sin+ 2 sin− 2 sin+ Ksin
2π TTTT
753と正弦関数の和として表すことができる．
図2(a)〜（d）は次数nを増やした場合の関数x(t)の近似値である．次数nを増やすと，関数x(t)を精度良く再現できることがわかる．

-2-1.5

1.5x(t)1n=30.5x(t)0.5x(t)x(t)1n=100-0.5-0.5-1-1-1.5-1.5-1-0.500.5-1

1-0.50t/Tt/T(a) n=1(b) n=3

1.5

0.5

1.5x(t)

1x(t)

1n=5

n=70.50.5x(t)x(t)100-0.5-0.5-1-1-1.5-1.5-1-0.500.5

1-1-0.5t/T

00.5

1t/T(c) n=5(d) n=7図2．2 三角級数による三角系列関数の近似
一般に区間[-T/2, T/2]を一周期とする周期関数のフーリエ級数はx(t ) =a0+22an =Tbn = 2T∫∫∞∑  a cos 2nTπt + b sin 2nTπt nnn =1x(t ) cos2nπtdt
Tx(t ) sin2nπtdt
TT /2−T / 2T /2−T / 2となる．
フーリエ級数を用いて関数を表すには項数を非常に多くとる必要があり，一般に20から30項が必要である．しかし，コンピュータが発達した今日では級数の項数を増やすことはそれほど困難ではなく，また，後述するFFTにより高速に計算することが可能である．
フーリエ変換で矩形波を表そうとすると，項数の増加とともに不連続点において振動が発生する．
これをGibbs現象という．

-3-1.41.21x(t)

0.8
0.6
0.4

n=1

n=20.2n=3n=50n=10

n=20-0.2
00.1

0.2

0.3

0.4

0.5t/T

図3 不連続点を持つ関数のフーリエ級数

2. フーリエ変換
基本周波数区間Tが有限値ならば，T/nは第n次高調波の周期であり，その逆数は周波数fnであり，

nT
1
δf =
T
fn =と表すことができる．
三角関数と指数関数の関係
cos θ =(e iθ + e − iθ )2(e iθ − e − iθ )sin θ =2から，a0x(t ) =+2∞∞∑∑

 i 2nπt  i 2nπt An expBn exp −+T  T  n =1
n =1an − ibn 1 T / 2
 i 2nπt =x(t ) exp −dt
 An =2T −T / 2T 
T /2
 i 2nπt  An = an + ibn = 1x(t ) expdt
2T−T / 2
 T ∫∫となるので，周波数を用いてフーリエ級数を表すと

-4-

∞

1
T → −∞
 T=−∞n∑ ∫ x(t) exp(−i2πnt / T )dt  exp(i2πnt / T )x(t ) = lim

T /2−T / 2∞

 T /2δf  x(t ) exp(−i 2πfnt )dt  exp(i 2πfnt )T →∞ −T / 2n = −∞= lim∑ ∫となる．したがって，区間Tを無限大に置き換えれば，∞∫ X ( f ) exp(i2πft )df
X ( f ) = x(t ) exp(−i 2πft )dt∫x(t ) =−∞∞−∞となる．これをフーリエ積分あるいはフーリエ変換と呼ぶ．

3. フーリエ変換の意味
プリズム

白色光

図4

光のスペクトル

プリズムを通過する白色光

暗室内に白色光（太陽光）を入れ，三角プリズムに当てると光の波長ごとに屈折率が異なるため，プリズムの背後に波長の長さ順に7色の光の帯が映し出される．
プリズムは白色光の中に含まれる各波長の光の成分毎に光を分解する役割をしている．分光された光の強さ（明るさ）は白色光に含まれる各波長の光の強さを表している．このようにプリズムは光の波長と強さを分離する装置である．さて，フーリエ変換では信号は周波数と振幅（強度）に分解されていることを思い出そう．すべての光の波長が集まった白色光では各波長の成分の持つ光の色がわからなくなる．これはちょうど，ランダム信号には特徴がないことに似ている．ランダム信号は副速な種々の波長の波の合成と考えられる．不規則波形はフーリエ変換という数学的プリズムを介して，信号をスペクトルに分解することに他ならない．

フーリエ変換スペクトルアナライザー

図5 フーリエ変換による信号の周波数分解-5-フーリエ成分X(f)は周期fの波の振幅を表し，|X(f)|2はその強さ，エネルギーを表している．光の分光との相似性から周波数fの成分波のエネルギーの分布関係|X(f)|2をエネルギースペクトルと定義する．有限区間Tにおける信号や周期関数の場合，|X(f)|2も有限と考えられるので，エネルギーも有限となり，物理的にも理解しやすい．しかし，区間Tが無限の場合，単位時間当たりの平均エネルギーをとって，パワースペクトル密度(Power Spectrum Density Function)P（ｆ）を定義する．
1
1
P ( f ) = lim  X ( f ) 2 = lim  X ( f ) X * ( f )T →∞  T T →∞  T

関数x(t)が確率過程の場合は，|X(f)|2の期待値についてのパワースペクトル1P ( f ) = lim E  X ( f ) 2T →∞  T

を定義することができる．
1T →∞ Tx 2 = lim∫T /2x2 (t )dt =−T / 2∞

∫ P( f )df−∞P(f1)df
スペクトル
P(f)

δf

周波数 f

図6 スペクトルの意味とエネルギーとの関係-6-4. Fast Fourier Transform
フーリエ変換は周期関数を三角関数の級数として表現するだけなく，ランダム信号の特徴を表現することが可能であり，工学的に有効な数学ツールである，しかし，信号を20から30項の級数に展開することは容易ではなく，コンピュータが発達する以前は，電気的なバンドパスフィルターにより信号を（まさにプリズムで分解するように）周波数毎に分析していた．1965年にCooleyとＴｕｋｅｙは計算時間を驚異的に短縮する方法を開発した．
長さＮの時系列信号x0,x1……xN-1の有限余弦フーリエ変換及び有限正弦フーリエ変換は以下のように表すことができる．
N −1

2πkj2
 Ak =xj cosNN

j =0

N −1

2πkj
2xj sin Bk =NN

j =0

∑

∑

(k = 0,1, K N / 2)
(k = 1,2, K N / 2 − 1)

ここで
C 0 = A0 / 2
CN / 2 = AN / 2
Ck = ( Ak − iBk ) / 2

N 0<k <CN − K = ( Ak + iBk ) / 22とすると，Ck =2
N

N −1∑ x exp − i 2πNkj  ( j = 0,1,2K N − 1)
jj =0

逆変換はN −1xj =∑ C exp i 2πNkj  ( j = 0,1,2K N − 1)
kk =0である．
Ｎ＝8のとき，ππ 2π 
W = exp i  = cos + i sin844とおけば,W2=i ,W4=-1, W5=-W, W6=--i,W7=-W3より-7- x 0  1
 x1  1
  
 x 2 1  
 x3 = 1
 x 4 1   x5  1
 x 6  1  
 x 7  1

1
WW2
W3W4W5W6W7111 W
1 − i
1 W3=
1 − 11 − W1 − i1 − W 3

1W2W4W6W8W 10
W 12W 14

1W3W6W9W 12
W 15
W 18
W 21

11i
W3−1 − i− i −W−11i −W3−1i− i −W

1W4W8W 12
W 16
W 20
W 24
W 28

1W5W 10
W 15
W 20
W 25
W 30
W 35

11−1 −W
1i−1 W3−11−1 W−i
1−1 W3

1W6W 12
W 18
W 24
W 30
W 36
W 42

1  C 0 
W 7   C1 W 14 C 2
 W 21 C 3
W 28 C 4
 W 35  C 5 W 42  C 6  W 49  C 7 

11  C 0 − i − W 3  C1 − 1 − i  C 2 − i − W  C 3 − 1  C 4
1
 − i W 3  C 5 −1i  C 6  iW  C 7 

行列の列を入れ替えて，偶数項と奇数項に分解する．
 x0  D0 + D 4  x1   D1 + WD5    x 2  D 2 + iD 6   3 x3 =  D3 + W D 7  x 4  D0 − D 4   
 x 5   D1 − WD5  x 6   D 2 − iD 6    x 7   D3 − W 3D 7 
C 0 
 D0
  D1   = T C 2,
C 4
 D 2  
C 6 
 D3

 C1 
 D 4  D5  = T C 3 C 5 
 D6
 
 
C 7 
 D7

1
1 1 11 i − 1 − i T=1 -1 1 -1 1 -i -1 i となる．フーリエ変換は簡単な加算計算に置き換えられる．
標本数がＮのとき，元のフーリエ変換に必要な演算数はN2のオーダーである．一方，上記の演算方法の場合，Ｎ＝Ｐ・Ｑとおくと，演算回数はＰＮ (logpN)となるので，もとの方法と比較すると演算回数比は P log 2 N となる． Ｎが十分大きければP=4のとき最も計算効率がよくなる．
log 2 P N-8-実際の計算手順：（1）データ数の決定：データ数Ｎを２のべき乗に選ぶ N=2p
データが２のべき乗になっていない場合は，データの一部を削るか，データの末尾に0を加えていく．
（2）ウィンド関数：フーリエ変換は本来，無限に長い区間Ｔに関する積分であるから，有限長のデータではデータの端部でデータが周期的に変化しないための誤差を生じる．データの端部における誤差を無くすために，ウィンド関数を用いる．データの始めと終わりをなだらかに0に近づけることで，無限に長い領域に0というデータがあるように修正する．
一般にはＨａｎｎｉｇ Windowが使われることが多い．（連続信号）
パルス信号の場合はRectangular Windowを用いる．
x( j ) = x(t )Wc(t ) 2πi Wc = 0.51 − cos
 N − 1 

Hanning
Rectangular1µ= T∫ W dt  = 3
T /2

22−T / 2

Wc = 1(3) パワースペクトルの算出：フーリエ係数からパワースペクトルを推定する．
~  k  ∆tP   = [Ak 2 + Bk 2]µT  N（4）平滑化：一般にスペクトル推定値は推定誤差が大きく，激しく振動する．スペクトルの分散を抑えるためには，平滑化処理が必要である．
統計平均 同一条件でM回の平均を取る．
分割平均 データをｌ個の部分に分割して各区間のスペクトルの平均値を取る．
周波数平滑 周波数空間での重み平均を取る．
0.03
1.50.025
P(f) m2/S3

u [m/s]

10.5
0-0.5

0.02
0.015
0.01
0.005

-10-1.5
00.05

00.1

100

t [sec]200300Frequency [Hz]図7 スペクトルの測定例-9-400500

ＦＦＴソフト等
ＦＦＴ解析用のハードウェアとしては，かつては小野測器製のＣＦシリーズが広く使われていたが，現在では，パソコンで直接信号を取得し，ソフトウェアでＦＦＴ処理する方法が一般的である．ＦＦＴのプログラム自体は簡単なので，A/Dコンバータがあれば自作も可能である．

各種のフリーソフトも作られているので，実験データの整理等には，これらのソフトを利用しても良い．
VirtualFFTアナライザは小野測器のＣＦシリーズのインターフェースを模擬した教育用ソフトでFFTアナライザ（波形解析装置）の操作方法や概念を学習することができる．以下のサイトからダウンロードが可能である．
http://rd.vector.co.jp/soft/win95/edu/se056666.html

WINDOWSの音源ボードと組み合わせて利用できるＷＳ及びＷＧ（信号発生器）を使うと音の周波数とスペクトルの関係がわかりやすい．このソフトも書きのサイトからダウンロードできる．
http://www.ne.jp/asahi/fa/efu/soft/ws/ws.html-10-#include<stdio.h>#include<math.h>#define MAXDATA 10000#define NF
4#define DT
0.0004#define N
2048#define WT
1.0/0.875#define PI
3.14156

dummyx[i]=sin(2.0*PI*500.0*dt*(double)i) ;dummyy[i]=0.0;
fprintf(fp,"%e %e\n",(double)i*dt,dummyx[i]);}fclose(fp);

int ipow2(),inv2pow(),atoi(),zero_y(),exit();int sines_2(),fft(),saw_wave();
double atof();strcpy(filename,"fft.d");break;
default:/*
USAGE: fft filename Sampling_Frequency FFT_Length
FillterLength*/
fp=fopen(argv[1],"r");nf=NF;
i = 0;for ( ;; ){if ( fgets(temp,78,fp)==NULL ){fclose(fp);break;}
if ( temp[0] == ';' || temp[0] <= 31 ||temp[1] == 0 ) continue;
if ( temp[0] == 'b' ) break;sscanf(temp,"%lf %lf",&t_tmp[i],&dummyx[i]);
i++;
}
n=i;
l=inv2pow(n);n=ipow2(l);void main(argc,argv)
int argc;
char *argv[];{FILE *fp;
int n,l,i,mf,nf,sample;
char filename[80],temp[256];
double x[MAXDATA], y[MAXDATA];
double sx1[MAXDATA], sy1[MAXDATA], sp1[MAXDATA];
double sx2[MAXDATA], sy2[MAXDATA], sp2[MAXDATA];
double sx[MAXDATA], sy[MAXDATA];
double fsp1[MAXDATA], fsp2[MAXDATA];
double fsx[MAXDATA], fsy[MAXDATA];
doublesumxy[MAXDATA],sump1[MAXDATA],sump2[MAXDATA];
double dummyx[MAXDATA];
double dummyy[MAXDATA];
doubleave,sum,dt,TN,OBT,FRE[MAXDATA],t_tmp[MAXDATA];
switch(argc){
case 0:exit(0);break;

dt=t_tmp[1]-t_tmp[0];break;}

case 1:
fp=fopen("time.d","w");
dt=DT;
l=inv2pow(N);n=ipow2(l);nf=NF;
for(i=0;i<n;i++){for(i=0;i<n;i++){x[i]=dummyx[i];y[i]=0.0;}-11-

/*hanning(x,y,n);*/fft(x,y,l,-1.0);}
int filter(z,fz,n,nf)
double *z,*fz;
int n,nf;{
int i,k,k1,k2;
int mf,m2,m;double sum;m2=nf*2;/*
printf("decomposition FFT data coresspond tooriginal\n");*/for(i=0;i<n;i++){sp1[i] =}(x[i]*x[i]+y[i]*y[i]);for(k=mf=0;k<n/2;k+=m2,mf++){
sum=0;
for(i=0;i<nf;i++){k1=k-i;k2=k+i;

OBT=1.0/(dt*(double)n);fp=fopen(argv[2],"w");for(i=0;i<n;i++){FRE[i]=OBT*(double)i ;

if(k1<=0)
k1= -k1;

fprintf(fp," %e %e %e %e\n",FRE[i],sp1[i]/OBT*WT,x[i],y[i]);}fclose(fp);if(k2>=n/2)
k2= n/2-(k2-n/2);

fft(x,y,l,1.0);
for(i=0;i<n;i++){printf("%e %e\n",dt*(double)i,x[i]);}sum+=(z[k1]+z[k2]);}
fz[mf]=sum/(double)nf;}}
return mf;

int hanning(x,y,n)double *x,*y;
int n;{
int i;
for(i=0;i<n;i++)
{
x[i]
*=
cos(2.0*PI*(double)i/(double)n));y[i]*=cos(2.0*PI*(double)i/(double)n));}}
int fft(x,y,l,f)double *x,*y,f;
int l;{
int i,i0,i1,j,l1,n,ns,n1,k,ipow2();
double s,c,s1,c1,sc,x1,y1,t;0.5*(1.0n=ipow2(l);n1=n/2;
sc=PI;

0.5*(1.0--12-j=0;
for(i=0;i<n-1;i++){if(i<=j){t=x[i];
x[i]=x[j];x[j]=t;}if(f<0.0)
for(i=0;i<n;i++){x[i] /= (double)n;y[i] /= (double)n;}}t=y[i];
y[i]=y[j];y[j]=t;int zero_y(y,n)
double *y;
int n;{
while(n!=0){*y++=0.0;n--;}}}k=n/2;
while(k<=j){j=j-k;
k /= 2;}j=j+k;}

int ipow2(l)int l;{
int n;
n=1;
while (l !=0){n *= 2;l--;}return n;}

ns=1;
while(ns<=n/2){
c1=cos(sc);s1=sin( f * sc);c=1.0;
s=0.0;
for(l1=0;l1<ns;l1++){for(i0=l1;i0<n;i0+=(2*ns)){i1=i0+ns;x1=x[i1]*c-y[i1]*s;y1=y[i1]*c+x[i1]*s;x[i1]=x[i0]-x1;y[i1]=y[i0]-y1;x[i0]=x[i0]+x1;y[i0]=y[i0]+y1;}
t=c1*c-s1*s;
s=s1*c+c1*s;
c=t;}
ns=2*ns;
sc=sc/2.0;

int inv2pow(n)int n;{
int l;
l= -1;
while(n!=0)
{
n/=2;
l++;
}
return l;}-13-