Document

信号とシステム
信号とシステム
システム科学専攻
林和則
• 第 2 章連続時間信号の変換
• 第 3 章サンプリングと z 変換
• 第 4 章 FFT とその応用
• 第 5 章アナログフィルタとディジタルフィルタ
• 第 6 章適応フィルタ
講義資料：http://www.msys.sys.i.kyoto-u.ac.jp/˜kazunori/ss.html
1
信号とシステム
第 4 章 FFT とその応用
• 実際のスペクトル解析で良く使用される離散フーリエ変換 (DFT)
について説明．
• DFT の高速算法としての FFT アルゴリズムの考え方を説明．
• FFT の使用に際しての具体的な留意事項を説明．
• FFT の実際的な応用として，スペクトル解析，畳み込み，相関関数
の計算，ケプストラム分析について概説．
第 4 章 FFT とその応用
2
信号とシステム
色々なフーリエ変換
フーリエ級数，フーリエ変換，DTFT (離散時間フーリエ変換)，DFT
(離散フーリエ変換) の関係
連続スペクトル
−∞ < t < ∞
離散スペクトル
0≤t≤P
第 4 章 FFT とその応用
連続時間
離散時間
−∞ < f < ∞
−fc < f < fc
フーリエ変換
DTFT
フーリエ級数
DFT
3
信号とシステム
フーリエ級数
x(t), 0 ≤ t ≤ P , −∞ < f < ∞
∞
∑
x(t) =
ck e
j 2π
P kt
k=−∞
∫
1
ck =
P
P
x(t)e−j P
2π
kt
dt
0
周波数 1/P の波を基本波とする調波構造を持った線スペクトル
1.5
1.5
1
0.5
1
0
−0.5
0.5
−1
−1.5
−4
第 4 章 FFT とその応用
−2
0
2
4
6
t
8
10
12
14
16
0
0
0.5
1
1.5
2
2.5
3
3.5
4
4
信号とシステム
フーリエ変換
x(t), −∞ < t < ∞, −∞ < f < ∞
X(jω) = F[x(t)] =
x(t) = F
−1
∫
∞
x(t)e−jωt dt
−∞
1
[X(jω)] =
2π
∫
∞
X(jω)ejωt dω
−∞
周波数軸上の連続関数である連続スペクトル
0.8
70
0.7
60
0.6
50
0.5
40
0.4
0.3
30
0.2
20
0.1
10
0
−0.1
−4
第 4 章 FFT とその応用
−2
0
2
4
6
8
10
12
14
16
0
−100
−80
−60
−40
−20
0
20
40
60
80
100
5
信号とシステム
離散時間フーリエ変換 (DTFT)
x(n), −∞ < n < ∞, −fc < f < fc = 1/2 (−π < ω < π)
∞
∑
X(jω) =
x(n)e−jωn
n=−∞
∫ π
x(n) =
1
2π
X(jω)ejωn dω
−π
|f | > fc では −fc ≤ f ≤ fc の繰り返しの連続スペクトル
10
70
9
60
8
50
7
6
40
5
30
4
3
20
2
10
1
0
第 4 章 FFT とその応用
0
5
10
15
20
25
30
35
0
−40
−20
0
20
40
60
80
100
6
信号とシステム
離散フーリエ変換 (DFT)
x(n), 0 ≤ n ≤ N − 1, −fc < f < fc = 1/2
Xk =
N
−1
∑
x(n)e−j N nk
2π
n=0
N −1
2π
1 ∑
x(n) =
Xk ej N nk
N
k=0
|f | > fc では −fc ≤ f ≤ fc の繰り返しの離散スペクトル
40
12
35
10
30
8
25
20
6
15
4
10
2
5
0
第 4 章 FFT とその応用
0
1
2
3
4
5
6
7
8
9
0
−4
−2
0
2
4
6
8
10
12
14
16
7
信号とシステム
DFT の演算
Xk =
N
−1
∑
x(n)e−j
N −1
1 ∑
j 2π
x(n) =
Xk e N nk
N
2π
N nk
n=0


X0
1

 
 X1  1

 

 
 X2  = 1

 
 .  .
 ..   ..

 
XN −1
1
k=0

e
1
1
−j 2π
N
−j 4π
N
e−j
e
e−j
2π×2
N
..
.
e
−j
2π(N −1)
N
4π×2
N
..
.
e
−j
4π(N −1)
N
···
···
···
..
.
···
1

x(0)






e
  x(1) 

2π(N −1)×2  
−j



N
e
x(2)





..
..



.
.


2π(N −1)(N −1)
−j
N
e
x(N − 1)
−j
2π(N −1)
N
W nk = e−j N nk を変換核とする離散的な直交変換
2π
第 4 章 FFT とその応用
8
信号とシステム
FFT (Fast Fourier Transform) の考え方
• N × N 行列と N 次元ベクトルの積の計算量 (乗算回数) は N 2 の
オーダー．
• FFT はバタフライ演算により計算量を N log2 N のオーダーに削減．
N = 2ν (ν; 正整数), 回転因子 WN = e−j N とし，奇数番目と偶数番目
を分離．
2π
Xk =
N
−1
∑
x(n)WN nk
n=0
N
2
=
−1
∑
n=0
第 4 章 FFT とその応用
−1
∑
N
2
x(2n)WN 2nk +
x(2n + 1)WN (2n+1)k
n=0
9
信号とシステム
N
2
Xk =
−1
∑
x(2n)WN 2nk +
n=0
−1
∑
−1
∑
N
2
x(2n + 1)WN (2n+1)k
n=0
N
2
=
x(2n)(WN 2 )nk + WN k
n=0
N
2
=
x(2n + 1)(WN 2 )nk
n=0
−1
∑
−1
∑
N
2
−1
∑
N
2
x(2n)W N nk + WN k
2
n=0
x(2n + 1)W N nk
2
n=0
= Ek + WN k Ok
• Ek : 偶数番目の標本値列の N/2 点 DFT.
• Ok : 奇数番目の標本値列の N/2 点 DFT.
Ek+N/2 = Ek , Ok+N/2 = Ok より, N 点 DFT を N/2 点 DFT 2 回で
実現．
第 4 章 FFT とその応用
10
信号とシステム
FFT のシグナルフローグラフ
N 点 DFT の分割
E0
x(2)
N/2-point DFT
x(0)
x(N-2)
X0
E1
0
WN
N/2-point DFT
x(3)
x(N-1)
Ek+N/2 = Ek
第 4 章 FFT とその応用
1
WN
E N/2-1
X N/2-1
O0
x(1)
X1
N/2
WN
X N/2
O1
N/2+1
WN
N/2-1
X N/2+1
WN
N-1
WN
O N/2-1
X N-1
Ok+N/2 = Ok
11
信号とシステム
バタフライ演算の修正
WN
N
2
+i
= −WN i を利用して乗算の回数を削減．
Ei
Xi
Ei
Xi
i
WN
Oi
i
N/2+i
WN
(a) butterfly operation
第 4 章 FFT とその応用
X N/2+i
Oi
WN
-1
X N/2+i
(b) modified butterfly operation
12
信号とシステム
分割の繰り返し
X0
x(0)
x(4)
N/4
X1
XN/4−1
(N−4)
x(2)
x(6)
N/4
WN0
WN2
−1
XN/4
−1
XN/4+1
WN2(N/4−1)
−1
x(N−2)
x(1)
x(5)
N/4
x(N−3)
x(3)
x(7)
x(N−1)
第 4 章 FFT とその応用
N/4
WN0
−1
WN2
−1
WN2(N/4−1)
−1
XN/2−1
WN0
WN1
−1
XN/2
−1
XN/2+1
WNN/4−1
−1
WN0
−1
WN1
−1
WNN/4−1
−1
U0
u(0)
XN/2+N/4−1
XN/2+N/4
u(1)
XN/2+N/4+1
XN−1
−1
U1
2 点 FFT
13
信号とシステム
定位置計算
演算結果の格納に元の記憶場所をそのまま使える
X0
x(0)
0
x(4)
W8
X1
-1
0
W8
x(2)
2
0
x(6)
W8
X2
-1
-1
W8
X3
-1
0
W8
x(1)
1
0
x(5)
W8
W8
-1
x(3)
x(7)
第 4 章 FFT とその応用
W8
-1
-1
W8
-1
X6
-1
X7
3
2
0
X5
-1
2
0
W8
W8
X4
-1
-1
W8
14
信号とシステム
ビット逆順配列
ビットを逆順にし，それを小さい順に並べている．
000
100
010
110
001
101
011
111
X0
x(0)
0
x(4)
W8
X1
-1
0
W8
x(2)
2
0
x(6)
W8
X2
-1
-1
W8
X3
-1
0
W8
x(1)
1
0
x(5)
W8
W8
-1
x(3)
x(7)
第 4 章 FFT とその応用
W8
-1
-1
W8
-1
X6
-1
X7
3
2
0
X5
-1
2
0
W8
W8
X4
-1
-1
W8
000
001
010
011
100
101
110
111
15
信号とシステム
FFT の計算量
段数は log2 N − 1 であり，1 段当たりの乗算は N/2.
X0
x(0)
x(4)
N/4
X1
XN/4−1
(N−4)
x(2)
x(6)
N/4
WN0
WN2
−1
XN/4
−1
XN/4+1
WN2(N/4−1)
−1
x(N−2)
x(1)
x(5)
N/4
x(N−3)
x(3)
x(7)
x(N−1)
第 4 章 FFT とその応用
N/4
WN0
−1
WN2
−1
WN2(N/4−1)
−1
XN/2−1
WN0
WN1
−1
XN/2
−1
XN/2+1
WNN/4−1
−1
WN0
−1
WN1
−1
WNN/4−1
−1
XN/2+N/4−1
XN/2+N/4
XN/2+N/4+1
XN−1
16
信号とシステム
逆 FFT
x(n) を Xk で置き換え，指数部の符号を反転して N で割ればよい．
Xk =
N
−1
∑
x(n)e−j N nk
2π
n=0
N −1
1 ∑
j 2π
x(n) =
Xk e N nk
N
k=0
• FFT アルゴリズムの中で回転因子 WN をその複素共役に変える．
• 出力値を N で割る．
FFT を使用する場合には，Xk∗ の DFT を FFT アルゴリズムで計算し,
得られた系列の複素共役をとって N で割る．
(N −1
)∗
1 ∑ ∗ −j 2π nk
x(n) =
Xk e N
N
k=0
第 4 章 FFT とその応用
17
信号とシステム
FFT の留意点: ナイキスト条件
標本化定理連続信号 x(t) のスペクトルが |f | < 1/2T でのみ非零な
らば
∞
∑
x(t) =
k=−∞
sin Tπ (t − kT )
x(kT ) π
T (t − kT )
2.5
70
60
2
50
1.5
40
1
30
0.5
20
0
10
0
−100
0
100
200
300
400
500
−0.5
−5
−4
−3
−2
−1
0
1
2
3
4
5
サンプリング間隔 T が 1/2 max |f | より大きい時，エイリアシング誤差
(折返し誤差) を生じる．
第 4 章 FFT とその応用
18
信号とシステム
アンチエイリアスフィルタ (Anti-Alias Filter)
一般に原信号に含まれる周波数成分はわからないため，あらかじめ適当
な LPF (low-pass filter, 低域通過型フィルタ) に通してから標本化する．
1
1
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
−3
−2
第 4 章 FFT とその応用
−1
0
1
2
3
0
−3
−2
−1
0
1
2
3
19
信号とシステム
FFT の留意点: 波形の切り出し (1)
• 時間範囲が有限でない連続信号の標本値系列は，その一部分を切り
出して FFT を適用．
• 標本値列の長さは通常 28 (= 256) – 211 (= 2048) 程度．
25
25
20
20
15
15
10
10
5
5
0
0
−5
−5
−10
−10
−15
−15
−20
−40
−20
0
20
40
60
80
100
120
140
160
−20
−40
−20
0
20
40
60
80
100
120
140
160
DFT は標本値列を周期系列として扱うことに注意．
第 4 章 FFT とその応用
20
信号とシステム
FFT の留意点: 波形の切り出し (2)
DFT で求めた振幅スペクトル
[dB]
[dB]
0
0
3.0 Hz
3.2 Hz
-10
-10
-20
-20
-30
-30
-40
0
5
10
15
frequency [Hz]
3.0 Hz の正弦波
第 4 章 FFT とその応用
20
25
30
-40
0
5
10
15
20
25
30
frequency [Hz]
3.2 Hz の正弦波
21
信号とシステム
いろいろな窓関数
両端での大きな不連続を防ぐため，窓関数 w(t) （両端で 0 に近い値）
をかける．
f (t) ⇒ w(t)f (t)
方形窓
wR (t) = 1
ハニング窓
wN (t) = 0.5(1 − cos 2π
ハミング窓
wM (t) = 0.54 − 0.46 cos 2π
ブラックマン窓
ガウス窓
第 4 章 FFT とその応用
t
)
P
t
P
t
t
wB (t) = 0.42 − 0.5 cos 2π + 0.08 cos 4π
P
P
(
)
t2
wG (t) = exp − 2 2
P σ
22
信号とシステム
窓関数の具体例
1
25
0.9
20
0.8
15
0.7
10
0.6
5
0.5
0
0.4
−5
0.3
−10
0.2
−15
0.1
0
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
−20
窓関数
ハニング窓
第 4 章 FFT とその応用
0
20
40
60
80
100
120
140
適用例
ハミング窓
ブラックマン窓
23
信号とシステム
周波数領域における窓関数
W (jω): スペクトル窓
fw (t) = w(t)f (t)
Fw (jω) = F[w(t)f (t)] = F[w(t)] ∗ F [f (t)] = W (jω) ∗ F (jω)
0
10
−2
ハニング窓
10
ハミング窓
−4
10
ブラックマン窓
−6
10
−8
10
−10
10
0
第 4 章 FFT とその応用
0.05
0.1
0.15
24
信号とシステム
標本点数の決め方
• アルゴリズムの入手し易さと計算効率から多くの場合 N = 2ν ．
• L が 2ν でない時，0 を加えて 2ν−1 < L ≤ 2ν = N を満たす N に
する．
{x(0), x(1), . . . , x(L − 1), 0, 0, . . . , 0}
1.8
1.6
[3, 2, 1, 0]
1.4
1.2
[3, 2, 1, 0, 0, 0, 0, 0]
1
0.8
[3, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
0.6
0.4
0.2
0
0
0.1
0.2
0.3
0.4
0.5
情報は増えないが，スペクトルの記述精度を向上するような効果がある
(周波数軸上の標本間隔を詰めて標本化関数による補間をしたのと等価)
第 4 章 FFT とその応用
25
信号とシステム
FFT の応用: スペクトル解析
• FFT はもともとスペクトル解析のための手法．
• 有限長で周期的であるとみなす．
• 現実の信号は時変．(例: 音声信号)
8000
6000
4000
2000
0
−2000
−4000
−6000
−8000
1
1.05
1.1
1.15
1.2
1.25
4
x 10
窓関数で切り出した信号を FFT することで，周波数特性の時間変化を
調べる．
第 4 章 FFT とその応用
26
信号とシステム
スペクトログラム
ランニングスペクトル (スペクトログラム) の例
1
6000
0.9
4000
0.8
0.7
2000
Frequency
0.6
0
−2000
0.5
0.4
0.3
−4000
0.2
−6000
0.1
0
−8000
0
1000
2000
第 4 章 FFT とその応用
3000
4000
5000
6000
7000
8000
9000
10000
0
500
1000
1500
2000
2500
Time
3000
3500
4000
4500
27
信号とシステム
FFT の応用: 畳み込み
N 点の標本値列 {x(n)}, {h(n)} の畳み込みの計算
y(n) = x(n) ∗ h(n) =
N
−1
∑
x(m)h(n − m)
n = 0, 1, . . . , N − 1.
m=0
• N 2 の乗算が必要．
FFT を利用．
Yk = Xk Hk
k = 0, 1, . . . , N − 1
{y(n)} = F −1 [Xk Hk ]
•
3
N log2 N の乗算で計算可能
2
→ ただし，通常の畳み込みではない (巡回畳み込み)
第 4 章 FFT とその応用
28
信号とシステム
巡回畳み込み
巡回畳み込み
y(n) =
N
−1
∑
x(m)h((n − m))N ,
((·))N : mod N
m=0
• FFT では {x(n)}, {h(n)} は周期的とみなされる．
• h((n))N = 0 (n < 0) でないので巡回畳込みは因果律を満たさない．
第 4 章 FFT とその応用
29
信号とシステム
FFT の応用: 線形畳み込み (有限長系列)
FFT で線形畳み込みを計算する方法 ({x(n)}, {h(n)} が有限長の場合)：
• {x(n)} の長さを Nx , {h(n)} の長さを Nh とする．
• N ≥ Nx + Nh − 1 とし， {x(n)}, {h(n)} の後ろに 0 を付加．
• n ≥ Nx + Nh では y(n) = 0.
• 周期性を排除して畳み込み計算が可能
第 4 章 FFT とその応用
30
信号とシステム
FFT の応用: 線形畳み込み (有限長と無限長系列)(1)
Overlap-add 法
zero-padding
L
h(n)
Q
x(n)
N
x (l) (n)
y (l) (n)
add
add
add
add
y(n)
第 4 章 FFT とその応用
31
信号とシステム
FFT の応用: 線形畳み込み (有限長と無限長系列)(2)
Overlap-save 法
zero-padding
L
h(n)
Q
x(n)
N
x (l) (n)
y (l) (n)
save
save
save
save
y(n)
第 4 章 FFT とその応用
32
信号とシステム
FFT の応用: インパルス応答と伝達関数
伝達関数 H(jω) はインパルス応答 h(t) のフーリエ変換 F[h(t)]
インパルス応答 → 伝達関数
系にインパルスを印加し，その応答出力 h(t) を標本化して DFT す
ると伝達関数 H(jω) が得られる．
• 標本化速度がナイキスト条件を満たす必要がある．
伝達関数 → インパルス応答
インパルス印加ができない場合，既知信号 x(t) を入力し，
その応答 y(t) を観測する．伝達関数は H(jω) = Y (jω)/X(jω) で
得られ，H(jω) を IDFT するとインパルス応答 h(t) が得られる．
第 4 章 FFT とその応用
33
信号とシステム
FFT の応用: 相関関数
自己相関関数
∫
x(t)x(t − τ )dt
rxx (τ ) =
相互相関関数
∫
rxy (τ ) =
x(t)y(t − τ )dt
ウィーナー-ヒンチン (Wiener-Khinchine) の定理
F[rxx (τ )] = P (jω) = X(jω)X ∗ (jω)
F[rxy (τ )] = X(jω)Y ∗ (jω)
第 4 章 FFT とその応用
34
信号とシステム
FFT の応用: ケプストラム分析
ケプストラム C(τ ): cepstrum, spec-trum からの造語
C(τ ) = F −1 [log |S(jω)|]
ケフレンシ τ : quefrency, fre-que-ncy からの造語
7
0.01
0.009
6
0.008
5
0.007
4
0.006
3
0.005
0.004
2
0.003
1
0.002
0
−1
0.001
0
0.5
1
1.5
2
2.5
0
3
0
1
2
3
4
5
4
対数パワースペクトラム
第 4 章 FFT とその応用
6
4
x 10
x 10
ケプストラム
35
信号とシステム
ケプストラムを用いた音声情報処理
音声波 s(t) は，音源波 v(t)，声道の特性 h(t)，口からの放射特性 r(t)
の畳み込み.
S(jω) = V (jω)H(jω)R(jω)
C(τ ) = F −1 [log |V (jω)|] + F −1 [log |H(jω)|] + F −1 [log |R(jω)|]
V (jω): 調波構造で基本周波数を示す．
H(jω): 典型的には山がいくつか連なった特性．
R(jω): 周波数軸上でなだらかな右上がり．
対数を取ることで，性質の異なる V , H, R を分離できる．
第 4 章 FFT とその応用
36

Download Report