スライド 1

Fourier変換
例１
sin 3x sin 5x sin 7x sin 9x
f (x)  sin x 



L
3
5
7
9
Plot[{Sin[x],Sin[3x]/3},{x,-Pi,Pi}]
f[x_] := Sin[x] + Sin[3x]/3
Plot[f[x],{x,-Pi,Pi}]
f[x_] := Sin[x]+Sin[3x]/3+Sin[5x]/5+Sin[7x]/7
Plot[f[x],{x,-Pi,Pi}]
例1（続き）
sin(2k - 1)x
f (x)  
2k - 1
k 1
n
f[n_,x_]:= Sum[Sin[(2k-1)*x]/(2k-1),{k,1,n}];
Plot[f[30,x],{x,-Pi,Pi}]
n   では
-p
0
p
2p
3p
4p
例2
sin 2x sin 3x sin 4x sin 5x
f (x)  sin x 



L
2
3
4
5
f[n_,x_]:= Sum[Sin[k*x]/k,{k,1,n}];
Plot[f[50,x],{x,-Pi,Pi}]
n   では
-p
0
p
2p
3p
4p
例3
cos 2x cos 3x cos 4x cos 5x
f (x)  - cos x 

L
2
2
2
2
3
4
5
f[n_,x_]:= Sum[(-1)^k Cos[k*x]/k^2,{k,1,n}];
Plot[f[50,x],{x,-Pi,Pi}]
Euler公式
exp(i )  cos  i sin 
虚数部
exp(i )
sin 

-1
単位円
cos 
1
実数部
単位円を N 等分
n


単位円を N 等分する点： exp i 2p
N

n  0, 1, 2, ..., N - 1
N  12
2
1 13
0 12
11 -1
10
-2
n=-12, 0, 12, 24, … は全て同じ点
n = -1, 11, 23, … は全て同じ点
公式
1
N
nm  1 (n  0)

exp i 2p


N  0 (n  0)

m 0
N -1
等比級数の和の公式からすぐに導ける。あるいは下の図のようにも理解できる。
2
4
1
0
0
11
10
n=1 の場合：これら 12個の点の和 = 0
8
n=4 の場合：これら 3個の点を4回づつ足す=0
Fourier変換
1
nm 

U nm 
exp i 2p

N 
N

N -1
N -1
1
nm 
nm 



U nmU mn   exp i 2p
 exp - i 2p


N m 0
N 
N 


m 0
1 N -1
(n - n)m 

  exp i 2p
   nn
N m 0
N



UU  1
U と U+ は互いに逆行列で

U U  1 も成り立つ
N個のデータ
a2
a1
a0
aN-1
n
-1
0
1
2
3
N-1
N
N+1
n<0 あるいは n>=N の部分は，周期的に繰り返すことにする
Fourier変換
an N  an
N -1
N -1
nm 

exp i 2p
an

N 

n 0
N -1
N -1
1
nm 


an  U nm Am 
exp - i 2p
Am

N 
N m 0

m 0
An N  An
1
Am  U mn an 
N
n 0
特に、 an が実数の場合
A  A-n  AN -n
*
n
変換
逆変換
MathematicaによるFourier変換
関数 Fourier[]
と InverseFourier[] を用いる
a = {0.2,1.0,1.2,3.3,2.1,-1.1,-2.2,-3.1,-1.2}
A = Fourier[a]
a2 = InverseFourier[A]
ファイルからデータを読むには ReadList[]を使う
SetDirectory["D:\\3d"]
a = ReadList["noise.txt"];
a[[1]]
(先頭のデータ）
ListPlot[a];
(読み込んだデータをプロット）
A = Fourier[a];
（Fourier変換を実行）
ListPlot[Abs[A],PlotRange->{0,6}];
Fourier成分 Am の対称性
Abs[Am]
フィルタリング：
ノイズとみなせるなら
ゼロにしてしまう




m
-N
-2 -1 0 1 2 3
An N  An
A  A-n  AN -n
*
n
よって
N-2
N
N+1
N-3 N-1
An  A- n  AN - n  AN  n
たとえば
A1  A-1  AN -1  AN 1
フィルタリング
高周波成分はノイズとみなして捨てる
（ＡのコピーをＢに作る）
B = A;
（ndata = データ数Ｎ）
ndata = Length[A]
k = 3;
Do[B[[m]]=0,{m,k+1,ndata-k+1}];
b = InverseFourier[B];
ListPlot[b];
Bk , Bk 1 , ..., BN -k をゼロにする。
Mathematica では B[[k+1]] = Bk なので，B[[k+1]], B[[k+2]], …, B[[ndata-k+1]]
をゼロにすることになる。
演習
noise2.txt で同様のフィルタリングを行ってみよ。
k の値を上手く設定する必要がある（Abs[A] の値を見
て，どこからどこまでをノイズとみなすかを判断する）。
フィルタリング後に逆変換し，グラフを描いてみる。
連続関数へ
an  f (t )
N 
f (t  T )  f (t )
1
Fm 
T
f (t ) 

T
0
exp( imt ) f (t )dt

2p
m 
m
T
m  0,1,2,3,...
 exp- i t F
m  -
m
変換
m
逆変換
三角関数で表す場合
1 T
Fm   exp( imt ) f (t )dt
T 0
1 T
  cos( mt )  i sin(mt )  f (t )dt
T 0
1
 C m  iS m 
2
2 T
Cm   cos( mt ) f (t )dt  C- m
T 0
2 T
S m   sin(mt ) f (t )dt  - S - m
T 0
逆変換を三角関数で表す
f (t ) 

 exp- i t F
m  -
m
m

 F0   exp- i m t Fm  expi m t F-m 
m 1
1
1
(c - is) (Cm  iS m )  (c  is) (Cm - iS m )
2
2
 Cm cos(mt )  Sm sin(mt )

1
f (t )  C0   Cm cos(mt )  S m sin(mt ) 
2
m 1
実関数、奇関数、偶関数
1
Fm 
T

T /2
-T / 2
exp( imt ) f (t )dt
2 T /2
Cm   cos( mt ) f (t )dt
T -T / 2
f (t ) が実数：
2 T /2
S m   sin(mt ) f (t )dt
T -T / 2
f (t )  f (t )
*
つまり
F  F-m
Cm ,実数
Sm
*
m
f (t ) が偶関数： f (-t )  f (t )
Sm  0
f (t ) が奇関数： f (-t )  - f (t )
Cm  0
例
f(t)
1
奇関数だから
Cm  0
0
-T/2
-1
t
T/2
T
2 T /2
S m   sin(mt ) f (t )dt (１周期分を積分）
T -T / 2
4 1 - cos(mT / 2)
4 T /2
  sin(mt )dt 
T
m
T 0
2p
2
m 
m だから S m 
(1 - cos mp )
T
pm
例（つづき）
cosmp  (-1)
m
を用いると
 4
m: 奇数
2

m
Sm 
(1 - (-1) )   p m
pm
 0 m: 偶数
sin(mt )
 f (t ) 

p m1:奇数 m
4

4  2p
1
6p
1 10p

  sin
t  sin
t  sin
t  ....
p
T
3
T
5
T

グラフで確認してみる
たとえば周期 T=1 とすると（T はもちろん他の値でも良い）
w = 2*Pi;
f[t_]:=(4/Pi)(Sin[w*t]+(1/3)Sin[3w*t]
+(1/5)Sin[5w*t] + (1/7)Sin[7w*t]);
Plot[f[t],{t,-1,1}];
演習
f(t)
1
0
-T/2
-1
T
t
T/2
Cm , Sm を求め、最初の数項の和のグラフを描く
ヒント：奇関数なので Cm=0
-T / 2  t  T / 2
の範囲で
2
f (t )  t
T
2 T /2
2
 S m   sin(mt ) tdt
T -T / 2
T
部分積分せよ

Download Report