講義

計測工学講義
第11回目
担当：西野信博
A3-012号室
[email protected]
プラズマ実験装置NSTX(Princeton)
1
本章の目次
• 前回演習の答え
• ３－１移動平均を用いた平滑化処理
– 単純移動平均
– 重みつき移動平均
– 適応化平滑化法
•
•
•
•
•
３－２積算平滑化処理
３－３装置関数の補正
３－４周波数領域における処理
３－５重畳波形の分離
時間が許せば、以下を概説する
– FFT
– 最大エントロピー法
2
演習
• 一般的に知られていることだが,
• 実空間でコンボルーション（たたみ込み積分）されている
形をフーリエ変換するとどのような演算になるか？

g (t ) 
 a(t ) f (t  t )dt
o
o
o

フーリエ変換
G()  A()  F ()
• 証明は，数学の教科書などを参考にしてください。
3
3－3 装置関数の補正処理
• 前２節で不規則雑音を対象としたSN比改善法を述べた。それ
から、観測時間を一定とした時、不規則雑音の減少に伴い、
波形自体も歪みの増加が起こることがわかっている。この原
因は、今までは数値フィルターの重みであった。実際の測定
においては、類似の歪みが測定装置によりもたらされる。例
えば、分光器に、ある単波長の光を入れたときの出力は、分
光器と光計測器で決まる分解能で決定される。
• これらを表すのに、装置のインパルス応答（入力を単位イン
パルスにしたときの装置の出力）を用い、これが、装置関数
（apparatus function/instrumental function）であった。
a(t )
δ(t)
t=0
時間遅れ
t=0
4
装置関数の補正とは？
• 前回説明した装置関数は、信号の鋭さをなまらせる働き
があった。
• それを補正するとは、所謂、ピークがなまった信号から正
しくピークを得ること
• すなわち、高周波成分を強調することになる。
• 装置関数 a( ) として、得られる信号 g ( ) は、真の信号
f ( ) と a( ) のたたみ込み積分で与えられた。
g ( ) 

 a(
o
) f (  o )d o ,


 a(
o
)d o  1

• 具体的な逆演算の方法を考える。
5
数値計算の考え方
• ディコンボルーションの演算は数値的に不安定であるが、
ここでは、代表的な方法を表す。
• 前式の離算化表示は g ( m ) 
a( n ) f ( m  n )

n
• あるいは
g ( m )   a( m  n ) f ( n )
n
• よって、行列表示すると
 gm    amn  fn 
• m=nでは、連立方程式を解く手法を用いる事ができる。
• 対角要素 aii  非対角要素 aij の時は、反復法が用
いられる。
6
反復法の基礎
• 連立方程式 AX  b
• これを、 X  MX + C の形に変形する。
 x1 
･ 
 
X  ･ 
 
･ 
 xn 
•
•
•
•
•
•
b1 
･ 
 
b  ･ 
 
･ 
bn 
a1n 


ann 
 a11
A  
 an1
次に、初期ベクトル X(0) を適当に決めて、
)
X( （ν
回目の反復値）を使って (  0,1, 2････）
( 1)
X( 1)  MX( )  C
と X
を順次求める。
もし、 X( ) が収束すれば元の方程式の解
収束するかどうかは、Mの形による。
X
となる。
7
行列の分解
•
•
•
•
今、Mのスペクトル半径（固有値の絶対値の最大値）
をρ (M) とすれば、ρ (M)＜1 の時収束する。
さて、A=D+L+U
但し、
0 
 a11


D
,
 0

a
nn 

 0

a21

L


 an1
ann 1
0




0 
 0 a12

U


0
a1n 


an 1n 

0 
8
ヤコービ法
• すると
(L + U  D)X = b
-1
-1
X = -D (L + U)X + D b
• であるから、
(ν+1)
X
-1
(ν)
-1
= -D (L + U)X + D b
• となる。この方法をヤコービ法と言う。
• Dは対角行列なので、D-1は簡単に求まる。成分ごとに書
くと（但し、aii≠0 の様にしておく、下式のkはダンピング
係数）
1
 b
( )
i
a
x

 ij j  a
aii  i
ii
bi
k 
( ) 
    aij x j  
aii  i
 aii
xi( 1)  
xi( 1)
9
ガウス・ザイデル法
• 次に、元の方程式を (D + L)X = -UX + b
• と書き換えると、この時
(ν+1)
(D + L)X
• これは、
(ν)
= -UX + b
X(ν+1) = -D-1 (LX(ν+1) + UX(ν) ) + D-1b
• とも書ける。
• つまり、
( 1)
xi
bi
1
( 1)
( )
  ( aij x j   aij x j ) 
aii j i
aii
j i
• これを、ガウス・ザイデル法と言う。
• （一般にはヤコービ法より収束の仕方が倍早い）
10
Successive Over Relaxation (SOR）
• さらに、変形してGauss-Seidel 法を使って
•
X( )から X( 1) をいったん計算して、
X( ν 1)   X( 1)  (1  )X( )
• とする方法がある。
• これを、SOR(Successive Over Relaxation)法といい、
0＜ω ＜2である。
X( 1)  (D + L)-1 (-UX(ν) + b)  (1  )X( )
• うまくすると、 Gauss-Seidel 法より収束が倍早い
11
Burger-Cittert の擬似ディコンボルーション
• 方法を以下に示す。
g
(k )
nh
( )   f ( k ) ( )a(   )
 0
f ( k 1) ( )  f ( k ) ( )  g ( )  g ( k ) ( )
• を用いる逐次近似法である。（実はJacobi法の k  aii とし
たもの）
k
• この場合の f ( ) の収束性は、フーリエ変換した式
Gk  Fk . A
•
Fk 1  Fk  G  Gk  Fk (1  A)  G
• から F  G 1  (1  A)  (1  A) 2 ･････ (1  A) n 
n


n  
G
Fn  G / 1  (1  A)    F
A
1 A  1
12
実際上の注意
• ｇ(ν )、a(ν )はともに測定誤差を含む事から、高周波成分
がある。すなわち、収束条件を必ずしも満たさない。
• また、ディコンボルーションの演算は高周波域を強調する
• 容易に発散する
• 対策
• 通常、演算前や計算途中に平滑化処理を入れる
13
3-4
周波数領域における処理
• 概要
– 時間空間座標を
x として表すと、通常データは f ( x) の
様に表現できる。ここで、 f ( x) が信号 s ( x) と雑音
n( x) の和、と相加的に表せたとする。
f ( x)  s( x)  n( x)
– この時、フーリエ変換(Fourier)によって周波数領域で表
せば
• 但し、
F ( )  S ( )  N ( )
F ( ) 


f ( x)ei x dx

• などと大文字の関数は小文字の関数のフーリエ変換を表す。
14
フィルター関数によるSN比改善
• この様に表した時、ｓ（ｘ）、n(x)の周波数成分が違う時、
S(ω ）とN(ω ）の分離は比較的容易となる。
• すなわち、信号成分は一定範囲の周波数成分しか含ま
ない場合が多く、雑音成分は一般的にどの周波数成分
も一様な強度の白色雑音で高い周波数成分を含んでい
る。


• 通常のSN比は
S ( )d / N ( )d




• F(ω )にフィルタ関数（重み関数）P(ω )を作用させれば、


•
 S ( ) P( )d /  N () P()d

• となり、異なった値が得られる。

15
例
生データｇ(ｘ)のスペクトル
フィルター関数Ｐ（ω ）
50
1.2
45
40
1
Ｘ
20
Ｐ（ω )
0.8
30
25
0.6
0.4
15
10
0.2
5
0
0
0
20
40
60
80
0
100
20
40
ω
• 生データのスペクトルを
フーリエ変換のプログラ
ム（ＦＦＴ）で求める
• フィルター関数Ｐ（ω）を
かける
• その結果を逆フーリエ
変換する（やはりＦＦＴ）
60
80
100
ω
フィルター通過のスペクトル
50
45
40
35
Ｆ＾（ω )
Ｆ（ω )
35
30
25
20
15
10
5
0
0
20
40
60
ω
80
100
16
フィルター関数の役割
• 従って、雑音成分N(ω )が比較的大きい周波数成分（通
常高周波成分）を滅衰させるフィルター関数をF(ω )に作
用させ、これを変数Xの領域にFourier逆変換すればSN
比の改善ができる。平滑化後の信号は

•
1
i x
Sˆ ( x) 
2
 F ( ) P( )e
d

• となる。
• P(ω )はSN比を大きくし、かつ、S(ω )に大きな影響がない
周波数領域のフィルターが望ましい。
17
たたみ込み積分との関係
• 関数は、フィルター関数P(ω )をフーリエ逆変換して、
1
p ( x) 
2
• を求め、これを用いて、
s '( x) 

 P( ) exp(i x)dw



f ( x ') p( x  x ')dx '

• と表せる。
• この表現は、コンボルーション（たたみ込み積分）で、前に
述べた平滑化の数学的基礎を与える。
• 従って、3.1節、3.2節の平滑化は、前に述べたように低周波
を通すフィルター関数となっていることがわかる。
• たたみ込み積分は、実の世界での処理、一方、フーリエ変
換に代表されるものはスペクトルの世界での処理である
18
使用上の留意点
• また、フィルター関数が周波数0（DC成分）において１に
なっていれば、

F (0) 


f ( x )e
i 0 x

dx   f ( x)dx

• の値が保存されるのでピーク毎の平滑化によってピーク
の強度（積分値）は変わらない
19
フーリエ変換を用いた方法の利点
• 平滑化と同様の手続きでディコンボルーション（分解能の改
善）も行える。
• 前節の結果は、観測データ g ( x)
• 測定装値への入力関数 f ( x)
• 装置関数 a( x) として、

g ( x) 

f ( x ')a( x  x ')dx '

• であった。これをフーリエ変換すると、 G( )  F ( )･A( )
• 我々に必要な情報は f ( x) なので
1
f ( x) 
2
• として、求められる。

 G( ) / A( ) exp(i x)d

20
フィルター関数の選択
• これは、平滑化の計算で P( )  1/ A( )
• となるフィルター関数を用いた事と等価である。
• 但し、前と異なるのは、このフィルターは高周波成分を強調す
る事である。
• 式(42)は、測定装置の分解能関数の周波数分布に逆比例した
フィルター関数を表すので、インバースフィルター（逆フィル
ター）とも呼ぶ。
• もし、ｇ（ｘ）がノイズを含まない理想的な信号なら、
• ω →∞で A(ω )→0であるが、G(ω )はA(ω )より速く0に収束し、
G(ω )/A(ω )は発散しない。
• しかし、実際の観測データはノイズを含むためG(ω )→ (not 0)
となり、発散する場合が多い。
• この様な理由から、1/A(ω )の高周波成分を滅衰させて計算さ
れている。
21
Ｇ(ω)とＡ(ω)、Ｗ(ω)の関係
Ｇ（ω）
1/Ａ（ω）
雑音成分
Ｇ（ω）
ω
1/Ａ（ω）
G( ) / A()が大きくなる
Ｗ(ω)
G( ) / A() の大きさを抑える
G( ) / W () を使用
ω
22
1
f ( x) 
2

 G( ) / A( ) exp(i x)d

• を変更し、部分的に分解能を高める方法もある。
A( )  A1 ( ) A2 ( )
• と表せる時に、A2(ω )に相当する分のみ分解能を改善す
る。すなわち、
1
f1 ( x) 
2





 G() / A () exp(i x)d
2

f ( x ')a1 ( x  x ')dx '
• この様にして、得られる関数は、真の入力分布に比べて
測定系の影響がとして与えられたものであるがに比べて
分解能を上げる事ができる。
フィルター関数の選択つづき
• 以上の様に、装置関数の及びA(ω )の正確な形がわから
ない場合でも、うまくそれなりの形のA(ω )を選ぶことがで
きれば、雑音成分が少ない信号なら分解能改善が可能
となる。
24
フーリエ変換を用いた方法の特徴より広い立場で
• フーリエ変換が必要であるから、移動平均に比べて多少演
算量は増える
• ただし、任意の周波数範囲の強調、あるいは、逆に除去が
可能である。
• 信号(雑音)の性質があらかじめ分かっている、あるいは、予
測できる場合にはきわめて有用である。実際は、スペクトル
解析によって、信号と雑音の区別をつける目安とする。
• 逆に言うと、周波数領域で同じ領域の雑音は除去できない
– 信号と雑音の区別がつかない
• このような場合、一般的な方法がないため、個別対応となる。
• 個別対応とは研究者の腕(能力)の見せ所
– 基礎をよく理解して、応用できる力がある
– 頭が柔軟
25
演習
• 信号処理のアルゴリズムで、移動平均を使用する場合と
フーリエ変換を使用する場合とを比べた場合、どちらの
方法が適用範囲（応用範囲）が広いか。
• 例を挙げて説明せよ。
26

Download Report