2

計測工学講義第9回目
担当：西野信博
A3-012号室
[email protected]
1
プラズマ実験装置NSTX(Princeton)
目次
• 前回演習の答え
• 第3章信号処理の方法
– スペクトル解析の応用とその例
2
演習
• ある周波数範囲以外の周波数を全く通さない理想的な帯
域フィルター（Band-pass filter ）を用意する。
• このフィルターを通過した白色雑音のパワースペクトル
S(ω)は、以下の式で表させる。ここに、ω1，ω2は定数であ
る。
 a
S (ω ) = 
0
0 < ω1 ≤ ω ≤ ω2
ω < ω1 or ω2 < ω
• Wiener-Khintchineの関係を使って、フィルターを通過した
白色雑音の自己相関関数を求めよ。
3
解答例
• 題意により、以下のパワースペクトルのフーリエ変換を
求めて見る。
 a
S (ω ) = 
0
0 < ω1 ≤ ω ≤ ω2
ω < ω1 or ω2 < ω
−ω1
ω2


iωτ
iωτ
iωτ
C (τ ) = ∫ S (ω )e d ω = a  ∫ e d ω + ∫ e dω 
 −ω2

−∞
ω1
∞
iωτ
−ω1
iωτ
ω2
e 
e 
= a
+a


i
τ
i
τ

 − ω2

ω1
a − iω1τ
a iω2τ
a
− iω2τ
iω1τ




= e
−e
+ e − e  = 2 ( sin ω2τ − sin ω1τ )

τ
iτ
iτ
4
続き
• よって，
 ω2 − ω1   ω2 + ω1 
τ  sin 
τ
sin 
C (τ ) =
τ
 2
  2

4a
sin (α + β ) = sin α cos β + cos α sin β
sin (α − β ) = sin α cos β − cos α sin β
A =α + β, B =α − β
A− B
A+ B
sin A − sin B
∴
= sin
cos
2
2
2
5
第3章信号処理の方法
−スペクトル解析の応用とその例−
• この章では、得られたデータから,測定者にとってより質
の良い信号を得るための処理方法を説明する.
• より良い信号とは何か？ Î余分な物がない信号
＝雑音がない信号のこと
• 計算機やパソコンの普及により,近年における信号処理
技術は,ディジタルデータ処理といって過言でない.
• 本章では、ディジタルデータの信号対雑音比(SN比、S/N）
改善処理について述べる。
6
３章の目次
• ３−１移動平均を用いた平滑化処理
– 単純移動平均
– 重みつき移動平均
– 適応化平滑化法
• ３−２積算平滑化処理
• ３−３装置関数の補正
• ３−４周波数領域における処理
• 時間が許せば、以下を概説する
– 重畳波形の分離
– FFT
– 最大エントロピー法
7
移動平均を用いた平滑化処理
• 移動平均とは、測定されたデータ(例えば,時系列データ)
に対して、各点での値をその周りのデータの平均で置き
換える方法である.
• 平滑化とは、雑音が入った測定値を雑音成分をなるべく
とって滑らかにすること
20
20
15
15
15
10
10
10
＋
5
-5
-5
-10
-10
-10
測定値
信号
49
46
43
40
37
34
31
28
25
22
19
16
13
7
49
46
43
40
37
34
31
28
25
22
19
16
13
7
10
4
1
49
46
43
40
37
34
31
28
25
22
19
16
13
7
10
4
1
-5
10
0
0
0
5
4
＝
5
1
20
雑音
8
アルゴリズムの図
• 単純移動平均
– アルゴリズムの説明
– 等間隔離散点に対する信号xiについて
1 m
x j = ∑ xi + j ,
N = 2m + 1
N i =− m
xj
– が単純移動平均データを与える.
• つまり、注目する点ｊのデータをｊ点を中心とした元の
データの平均値に置き換える.
• では,実際にどのような演算をしているか見てみる.
9
単純移動平均の考え方
1.
平均を取る区間幅を決め(下図では9点），最初の区間を設定する
2.
3.
4.
その区間で平均をとる
平均値を区間の中央点の値とする新しいデータを作る
平均を取る区間をひとつずらして，3から繰り返す
20
15
10
5
49
43
37
46
40
35
37
34
31
28
25
22
19
16
13
10
7
4
1
0
-5
-10
測定値
区間幅9点
16
14
12
10
8
6
4
2
43
41
39
33
31
29
27
25
23
21
19
17
15
13
11
9
7
5
3
1
0
-2
10
-5
-10
信号
49
46
43
40
37
34
31
28
25
22
19
16
13
10
7
4
1
15
10
5
0
13点
-4
-2
18
16
14
12
10
8
6
4
2
0
9
27
29
29
31
33
45
43
41
39
37
35
25
27
43
41
39
37
35
33
31
23
21
19
17
15
13
11
25
23
21
19
17
15
13
20
11
測定値
7
-10
9
9点
5
0
7
10
3
46
43
40
37
34
31
28
25
22
19
16
13
10
7
4
1
15
1
20
3
5
• 以下の測定値は、ガウス分布
に乱数を加えて人為的に作っ
た模擬測定値である。
5点
1
-5
49
46
43
40
37
34
31
28
25
22
19
16
13
10
7
4
1
単純移動平均の適応例
18
16
14
12
10
8
6
4
2
-2
0
-4
5
16
14
12
10
8
6
4
2
0
-2
11
単純移動平均の期待値と分散
• xnをn番目の標本値とし、その平均値、
E ( xn ) = µ
2は、
分散をとすると、標本平均の分散s
σ 2 ( xn ) = σ 2 x
0
N N
1


2
s 0 = 2 E  ∑∑ ( xn − µ )( xm − µ ) 
N
 n =1 m =1

• で表わされた。
[
]
E xn xm
• 但し、は標本間の相関である。
σx
1
1
2
2
完全不相関 s0 = 2 E  ∑ ( xn − µ )  = 2 Nσ x =
N
N
N
2
2
完全相関
1
1 2 2
2
s0 = 2 E  ∑∑ ( x − µ )  = 2 N σ x = σ x 2
N
N
2
12
実際は不完全な相関である
• 上記の中間の時、すなわち、各標本点に部分相関があると
し、各標本が周期T,間隔Δtの時系列データとすると、
σ x2
2
k t
) Cx (k t ) − µ 2  t
s0 =
+ ∑ (1 −
N T
T
2
• となる。ここで、Ckは自己相関関数
• である
C x (k t ) = E  xi x j  , k = j − i


ポイントは,標本値が互いに独立である時に、処理後の
1
雑音の標準偏差はに比例して減少するが、標本値間
N
の相関が大きいほど処理効果が小さくなる。
13
単純移動平均の特徴
• N→大にすると、通常、雑音成分は減少するが、信号波形も
平滑化される欠点がある。
• T＝一定とし、N→∞とすると、
2
τ
s0 = ∫ (1 − )(C x (τ ) − µ 2 )dt
T 0
T
T
2
(k t → τ )
• 従って、Tが大の時、
∞
s0
2
Px (0)
2
2
(C x (τ ) − µ )dτ =
∫
T 0
T
Px (ω )
• ここに、は雑音パワースペクトル密度である。
14
単純移動平均の特徴
• 上式からもわかる様に、平滑化処理は低周波成分に
対して有効ではない。（例、1/fゆらぎ）
• 実用例
この様に、移動平均法は、高周波成分の雑音の中から低周
波の信号成分を抽出するのに有効となる。
15
単純移動平均の流れ図
• 以下に,単純移動平均の流れ図（フローチャート）の例を示す.
N個のデータ
NQ個(奇数)で平均
Begin
End
YES
I=（NQ+1）/2
NO
I>N-NQ/2
J=−NQ/2, SUM=0
I=I+1
SUM=SUM+DATA(I+J)
J=J+1
NO
J>NQ/2
YES
SMOOTH(I)=SUM/NQ
16
説明
• N個のデータが配列DATAに入っている。DATAの添え
字は１から始まる。
Begin
I=（NQ+1）/2
• は、NQ個の平均を取るので、その中心の点は、DATA
（I）である。配列DATAの添え字は0以下はないとした。
J=−NQ/2, SUM=0
• 平均を取るための初期値操作
17
DATAの中のNQ個の平均を取る
• 以下の式でまずNQ個のデータの総和を取っている。
SUM=SUM+DATA(I+J)
J=J+1
J>NQ/2
NO
NQ個数えて終了する
YES
SMOOTH(I)=SUM/NQ
総和をNQで和って、
平均を取る
その後、新しい配列
SMOOTHに入れる
18
平均値を次々求める
• 平均を取るデータ点を進ませる(移動)
あれば、平均を
取る操作へ
I=I+1
NO
I>N-NQ/2
YES
End
まだ、処理すべきデータが
あるかどうか判定する。
19
重みつき移動平均
• 前節の単純移動平均は、平均を取る点数Nを大きくと
ると,信号歪も大きくなりやすい欠点があった.
• 実際は、平滑化の性能を落とさず,簡単な演算処理で,
なるべく信号歪を小さくしたい.
• そこで,考え出されたのが重みつき移動平均で、代表
的な平滑化法の一つである.
• 今,重み係数Wiとすると
o jw =
u
1
u
∑W
i =− L
・∑ x j +・
i Wi
i =− L
i
• ここで、L及びUは重み関数の下限と上限を表してい
る。
20
では、最適な重みは？
• 平滑処理後の不規則雑音の分散を示す式は、
U
U −1
i =− L
i =− L
s0 2 = σ x 2 ∑ W 2 (i t ) + 2 ∑
U
2


−
W
(
i
t
)
W
(
j
t
)×
C
(
k
t
)
µ
∑
 x
j =i +1
•
•
•
•
但し、Wは正規化した重み関数である。
重み関数は、処理対象により、適当な係数が選ばれる。
逆に言うと、処理対象毎に最適値を捜す必要が生じる。
また、真の実時間処理を望むのであれば、非対称（未
来のデータがないため）な重みとなるが、多少でも遅
れを許しうるのであれば対称形重みを用いることがで
きる。
• そこで、汎用性のある重みがいくつか提案されている。
21
多項式適合法
• データ(標本点)が多項式に適合するとして、2乗誤差
を最小とするように係数を決める方法で,多項式適合
による重み係数(Savitzky-Golayの表)がある.
•
•
•
•
•
2次3次多項式：2次多項式と3次多項式は同じ係数
4次5次多項式：4次多項式と5次多項式は同じ係数
・
・
一般的な2n次2n+1次多項式の求め方を以下に述べる.
22
多項式適合法の係数の求め方
• 等間隔データ2m+1個を考える.
• Pn をiについてのｎ次の直交多項式とすると、任意の
多項式Fiは、
Fi = a0 P0 + A1 P1 +・・・・・+an Pn ,n<2m+1
• と表されるÎベクトルの近似を思い出せ！
• ここで、整数の関数で性質の良い直交多項式にGramの
多項式を用いると、
(k )
n
n
+
k
i
Pn,2m(i)=∑ (−1) k ( )(
) (k )
k
k 2m
k =0
n
• •
• 但し、 i ( k ) = i (i − 1)・・・・・(i − k + 1)
23
Gramの多項式の性質
• 以下に、Gramの多項式の例を挙げる
P0,2 m (i ) = 1
i
i
= 1−
P1,2 m (i ) = 1 − 2
m
2m
i
i (i − 1)
+6
P2,2 m (i ) = 1 − 6
2m
2m(2m − 1)
• 直交性より
2m
∑P
i =0
j ,2 m
(i ) Pk ,2 m (i ) = 0, j ≠ k
24
最小２乗適合法を利用
• 観測値（標本の値）をとし、2m+1点において、ｎ次
多項式を最小２乗適合する。
2m
Q = ∑ (a0 P0 + a1 P1 +・・・+ an Pn − Yi ) 2
i =0
• を最小とする様に係数akを求める。よって、
∂Q
=0
∂ak
∴
2m
2m
i =0
j =0
J ≠k
2
(
a
P
∑ k k + ∑ a j Pj Pk − Yi Pk ) = 0
• 直交性を利用して、
2m
2m
i =0
i =0
ak ∑ Pk 2 − ∑ Yi Pk = 0
2m
∴ ak =
∑Y P ,
i =0
2m
i k 2m
(i )
2
P
∑ k ,2m (i)
i =0
25
最小２乗適合法続き
• すると、Yiの推定値は、
Yi

1  2m
Yi = ∑ ak Pk ,2 m (i ) = ∑  ∑ Y j Pk ,2 m ( j )  × Pk ,2 m (i )
k = 0 Ak  j = 0
k =0

2m
• 但し、 A =
k
2m
2m
∑P
j =0
k ,2 m
2
( j)
• 最後に,以下の式が得られる
2m P


k ,2 m ( j )
ˆ
Yi = ∑ Y j × ∑
× Pk ,2 m (i ) 
j =0
 k =0 Ak

2m
26
最小２乗適合法続き
• n=2 の時、すなわち、２次多項式適合の場合は、重
みWiが得られている。
• 今、推定したい観測値の点の指標を0として2ｍ+1点
でこの値を推定したいとすると
• 各標本点の指標iは
i = − m, − m + 1,・・・・・0,1, 2,・・・・, m
• となる。この時、各点の観測値に対する重みは、
W (i ) = {3m(m + 1) − 1 − 5i 2 } , i = −m ∼ m
2
W
(
i
)
=
(4
m
− 1)(2m + 3) / 3
∑
27
実際に利用する場合の注意点
• 汎用性の重み係数は、観測値相互の相関を無視して
決められている。従って観測値の波形歪みをできる
だけ小さくして雑音を除くためには、最適な重み係
数を選ぶ必要がある。
• あらかじめ観測値の波形がわかっている場合は、そ
の波形歪みが重み係数によってどう変わるか、そし
て、雑音の軽減度はどの程度かを見積り、最適な重
み係数を選定する。
• 通常,重みつき平均を取る点数は25点以下の場合が多
い.(データ点数は1024点など、2のべき乗である)
28
Savitzky-Golayの係数表
• 以下に、2次3次多項式適合による平滑化重み係数
を挙げる
• 重みは中心対称であるので、片方のみ書いてある.
• 中心点ｘ０である。
x0
x1
x2
x3
x4
x5
x6
x7
x8
5
17
12
-3
7
7
6
3
-2
9
59
54
39
14
-21
11
89
84
69
44
9
-36
13
25
24
21
16
9
0
-11
15
167
162
147
122
87
42
-13
-78
17
43
42
39
34
27
18
7
-6
-21
29
-10
-2
-4
46
43
0
40
5
37
13点
34
10
31
15
28
-4
25
20
22
-2
19
-10
16
-5
16
14
12
10
8
6
4
2
0
13
10
7
4
43
41
39
37
35
33
31
29
27
25
23
21
19
46
43
40
37
34
31
28
25
22
19
18
16
-5
17
-10
15
13
0
11
5
13
10
10
15
9点
9
20
7
-5
5
49
46
43
40
37
34
31
28
25
22
19
16
13
0
7
1
7
10
5点
3
1
49
46
43
40
37
34
31
28
25
22
19
16
13
10
4
15
4
49
46
43
40
37
34
31
28
25
22
19
16
13
10
7
1
20
1
-5
7
雑音
4
1
信号
4
測定値
1
2次3次多項式適合の重みつき移動平均の例
20
15
10
5
10
5
0
18
16
14
12
10
8
6
4
2
0
30
7
0
20
15
15
10
10
5
5
0
0
-5
-5
-5
-10
-10
-10
39
20
37
-10
35
-10
33
-10
-5
31
-5
-5
29
43
41
39
37
35
33
31
0
29
0
27
5
25
5
27
10
25
10
23
15
23
15
21
20
21
-10
20
19
-10
19
7
46
43
40
37
34
31
28
25
22
19
16
13
10
7
4
1
46
43
40
37
34
31
28
25
22
19
16
13
10
7
4
1
49
46
43
40
37
34
31
28
25
22
19
16
13
10
重みつき
17
-10
17
-5
15
0
15
0
13
0
11
5
13
5
11
10
9
10
5
9
10
7
15
7
15
5
15
3
1
20
4
20
5
43
41
39
37
35
33
31
29
27
25
23
21
19
17
15
13
11
9
7
1
20
3
1
39
37
35
33
31
29
27
25
23
21
19
17
15
13
5
11
10
9
15
7
20
5
-5
5
3
1
0
-5
10
5
49
46
43
40
37
34
31
28
25
22
19
16
13
15
3
7
10
20
1
49
46
43
40
37
34
31
28
25
22
19
16
13
10
雑音
4
1
測定値
4
信号
1
単純移動平均と
2次3次多項式適合重みつき移動平均の比較
• 乱数の大きさを変化させて比較した例１
単純
5点
9点
13点
31
20
20
15
15
10
10
10
5
5
5
0
0
0
39
37
35
33
31
29
27
25
23
21
19
17
15
13
0
-5
-5
-10
-10
-5
-5
-5
-10
-10
-10
39
0
37
5
35
5
33
10
31
10
29
15
27
15
25
15
23
20
21
20
43
41
39
37
35
33
31
29
27
25
23
21
19
-10
19
-10
17
-5
17
-5
15
-10
7
46
43
40
37
34
31
28
25
22
19
16
13
10
7
4
1
46
43
40
37
34
31
28
25
22
19
16
13
10
7
49
46
43
40
37
34
31
28
25
22
19
16
13
10
重みつき
15
0
13
0
11
-5
13
5
11
5
9
0
9
10
7
10
7
5
5
15
4
15
3
1
10
5
20
1
43
41
39
37
35
33
31
29
27
25
23
21
19
17
15
13
11
9
7
5
20
3
15
1
20
3
20
4
20
11
-10
9
-5
7
0
1
5
5
10
3
49
46
43
40
37
34
31
28
25
22
19
16
13
10
7
1
15
1
49
46
43
40
37
34
31
28
25
22
19
16
13
10
7
雑音
4
1
測定値
4
信号
1
単純移動平均と
2次3次多項式適合重みつき移動平均の比較
• 乱数の大きさを変化させて比較した例2
単純
5点
9点
13点
32
重み係数の形による分類
• いろいろな移動平均を重みによって分類すると,例えば下
図のようになる.
W(j)
W(j)
W(j)
ｊ
0
単純移動平均
ｊ
0
2次3次多項式適合
0
ｊ
擬似RCフィルター
33
演習
• 単純移動平均のアルゴリズムを使用して、重み付移動
平均のアルゴリズムを作成したい。
• どの部分を変更すると合理的か？
• 変更部分を具体的に示せ。
34
単純移動平均の流れ図
• 以下に,単純移動平均の流れ図（フローチャート）の例を示す.
N個のデータ
NQ個(奇数)で平均
Begin
End
YES
I=（NQ+1）/2
NO
I>N-NQ/2
J=−NQ/2, SUM=0
I=I+1
SUM=SUM+DATA(I+J)
J=J+1
NO
J>NQ/2
YES
SMOOTH(I)=SUM/NQ
35

Download Report