大西
昇（名古屋大学）
信号処理
1.
はじめに
1
1.1 信号とは
1.2 信号処理とは
1.3 信号処理の内容
1.4 信号処理の応用
1.5 信号処理の歴史
1.6 信号の分類
1.7 ディジタル処理の長所と短所
2.
アナログ信号の変換
4
2.1 フーリエ級数展開
2.2 フーリエ変換
2.3 ラプラス変換
3.
不規則信号
14
3.1 不規則信号の記述
3.2 集合平均と時間平均
3.3 特性関数
3.4 キュムラント関数
3.5 ガウス性不規則信号
3.6 定常性とエルゴード性
3.7 相関関数とスペクトル
3.8 不規則信号と線形システム
4.
信号のディジタル化
22
4.1 標本化と量子化
4.2 標本化定理
4.3 エイリアシング
4.4 信号の量子化
4.5 標本化周波数の変換
4.6 D/A 変換
5.
離散時間信号とその変換
29
5.1 離散時間信号の表現
5.2 離散時間フーリエ変換
5.3 離散フーリエ変換
講義信号処理目次.docx
i
大西
昇（名古屋大学）
5.4 z 変換
5.5 逆 z 変換とその求め方
6.
離散時間システム
41
6.1 システムとは
6.2 システムとしての重要な性質
6.3 線形時不変システム
6.4 差分方程式と離散時間システム
6.5 システムの分類
6.6 システム関数
6.7 離散時間システムの周波数特性
6.8 システムの構成
6.9 システムの接続
7.
フィルタ
51
7.1 フィルタとその種類
7.2 周波数特性の性質
7.3 代表的なアナログフィルタ
7.4 ディジタルフィルタの設計法
7.5 直線位相 FIR フィルタ
7.6 低域通過以外の他のフィルタの導出
8.
高速フーリエ変換(FFT)
63
8.1 高速フーリエ変換とは
8.2 FFT アルゴリズムの導出
8.3 ビットリバーサル
8.4 FFT の応用
9.
線形予測法
69
9.1 時系列の解析
9.2 線形予測モデル
9.3 予測係数の決定
9.4 偏自己相関係数
9.5 偏自己相関係数から予測係数
9.6 Levinson-Durbin の方法
9.7 予測次数の最適化（AIC 基準）
9.8 適用例（音声のピッチとフォルマント）
講義信号処理目次.docx
ii
大西
10. 適応信号処理
昇（名古屋大学）
78
10.1 適応システム
10.2 適応信号処理の例
10.3 適応フィルタ
10.4 適応アルゴリズム
11. ケプストラム解析
86
11.1 ケプストラム解析とは
11.2 ケプストラム
11.3 音声のケプストラム解析
11.4 ケプストラム解析によるエコー検出
付録
参考文献
90
93
索引
95
講義信号処理目次.docx
iii
大西
昇（名古屋大学）
講義信号処理目次.docx
iv
大西
1
昇（名古屋大学）
はじめに
信号は，何らかの情報を我々にもたらす．信号から情報を抽出する信号処理の内容，応用，歴
史を説明する．つぎに，現在の主流である，コンピュータでの信号処理すなわちディジタル信号
処理の長所と短所を説明する．
1.1 信号とは
信号は，電圧，音圧，光などの物理量の時間的変化で，情報（抽象概念）を運ぶ．例えば，
音声から音素や単語，レーダの反射波形から対象の位置・姿勢・材質，脳波から脳の活動の有
無，地震波から地震の大きさと位置，株価から企業の力・将来性，顔画像から個人性や喜怒哀
楽の情報を得ることができる．
1.2 信号処理とは
狭義の信号処理は，信号から目的とする情報を取り出すために，最も効果的な処理・解析を
施すことである．一方，広義の信号処理は，情報の抽出に加え，信号の加工・変形・合成をも
含む．信号処理の手法を評価する基準は，処理結果の精度と処理の複雑さ（計算量）である．
1.3 信号処理の内容
信号処理の内容を表 1.1 にまとめる．
表 1.1
信号処理の内容
ろ波（フィルタリング）
雑音除去，信号の強調．
解析（分析）
平均・分散などの統計量，スペクトル，定位．
予測/システム同定
信号を生成するシステムをモデル化（システムの性質を決定）．
認識（信号の分類・識別）
音声認識，文字認識，表情認識，物体識別．
加工/変換
音声変換，音声に感情付与，信号圧縮．
合成（生成）
音声合成，音場再現，モーフィング(morphing)，コンピュータグラフ
ィックス（CG）．
参）morph:２つの異なった形状を関係付けて，なだらかに形状を変形させる．語源は metamorphosis．
1.4 信号処理の応用
信号処理は，計測と制御，通信，音声映像メディア，医療機器などで，幅広く活用されてい
る．計測としては，レーダ，ソナー，地質探査，宇宙探査，地震波が，制御ではシステム同定
や適応制御が挙げられる．通信では，符号化，エコーキャンセラー，等化器が応用例である．
音メディアとしては，音声認識，音声符号化，音響処理，音場再現が，映像メディアとしては，
画質改善1，画像符号化2，映像内容解析(インデキシングやメタデータ付与のため)，物体検出，
物体認識が挙げられる．そして，心電図，脳波，コンピュータ断層撮影画像，血流，筋電図な
どの，医用生体信号の処理の基礎である．
1
ディジタルフラットテレビでの，明るさの制御でまぶしさの低減，動画ボケの補正などがある．
2
ディジタル放送の符号化は MPEG2 を使用．
信号処理 1 章-4 章.docxx
- 1 -
大西
昇（名古屋大学）
1.5 信号処理の歴史
信号処理の歴史を年表式にまとめると次のようになる．
19 世紀初め
フーリエ解析(by J.B.J.Fourier[1768-1830])
1936 年
PCM(pulse code modulation，パルス符号変調)
1939 年
ボコーダ(voice coder，音声符号化装置)(by H.Dudley)
1940 年代
予測・ろ波の（不規則信号を取り扱う）理論(by N.Wiener[1894-1964])
1949 年
標本化定理(by C.E.Shannon[1916-2001]と染谷勲)
1960 年頃
最小 2 乗平均(LMS)アルゴリズム(by B.Widrow)
1963 年
双１次ｚ変換によるディジタルフィルタ設計理論(by J.Kaiser)
1965 年
高速フーリエ変換(by J.W.Cooley and J.W.Tukey)
1969 年
PARCOR(by 板倉文忠)
1971 年
マイクロプロセッサ
1970 年代
線形予測符号化
1980 年代初め
ディジタルシグナルプロセッサ
1980 年代後半
ブラインド信号分離
1.6 信号の分類
信号の分類は，連続か離散，および確定か非確定，の 2 通りある．
1) 連続と離散
連続と離散の分類は，時間軸および振幅軸のそれぞれでなされ，表のように 4 つに分類され
る．気温や電圧などの物理信号は，時間軸および振幅軸の両方において連続的で，アナログ
信号と呼ばれる．一方，時間軸および振幅軸の両方において離散的な信号をディジタル信号
という．コンピュータ内での信号はすべてディジタル信号である，
時間
振幅
連
続
連
続
アナログ信号
離
x(t)
散
多値信号
（連続時間信号）
離
散
標本値系列
x(nT)
ディジタル信号
（離散時間信号 x(n)）
n：整数，T：標本間隔
参）アナログ(analog = analogue) 相似型，計量型（温度はそのアナログである電圧で表現できる）
．
ディジタル(digital)指の，計数型．
2) 確定と非確定
この分類では，確定信号と不規則信号に分けられる．確定信号は，信号の値が時間 t の関数
として完全に確定している信号で，信号の周期性により，さらに次のように細分化される．
周期信号
同一の波形を無限に繰り返す．
概周期信号
周期がほぼ一定，あるいは長い時間繰り返す．
信号処理 1 章-4 章.docxx
- 2 -
大西
非周期信号
昇（名古屋大学）
周期を持たない．
一方，不規則信号は，値が確率的な法則に従って不規則に変動する信号で，統計的処理が必
要である．不規則信号を，時間的に変動する確率現象である確率過程（不規則過程）とも言う．
1.7 ディジタル処理の長所と短所
・長所
精度 演算の係数語長と信号語長を長くすれば，高精度の特性を実現できる．
再現性と安定性（信頼性）
回路素子のバラツキの影響がないので，特性の再現性が高い．
また，温度変化や経年変化による特性の劣化がない．
経済性
集積化により価格が指数関数的に減少する．
柔軟性/汎用性
ハードウェアを変えることなく，係数値の変更で，いろいろな特性を実現
できる（仕様変更に対する柔軟性）
．計算機を用いた汎用性の高い処理が可能である．
・短所
処理速度
処理ステップ数に比例した時間が必要であるので，一標本点に対する処理が標本
間隔内に完了しない場合，実時間処理を実現できない．
誤差 標本化・量子化という離散化による誤差が発生する．
信号処理 1 章-4 章.docxx
- 3 -
大西
2
昇（名古屋大学）
アナログ信号の変換
アナログ信号の変換として，周期関数のフーリエ級数展開，非周期関数のフーリエ変換をま
ず説明する．いずれも，時間関数の周波数領域での表現で，フーリエ級数展開は，周波数軸で
の離散的な表現であるのに対し，フーリエ変換は連続的な表現である．次に，フーリエ変換の
拡張としてのラプラス変換を説明する．この変換は，周波数軸を虚軸にもつ複素平面への変換
である．
2.1 フーリエ級数展開
1) 関数の直交展開
基本周期 T を持つ周期信号 f(t)は，次式のように無限級数の和で表現でき，これをフーリエ
級数展開（Fourier series expansion）という．
0 2 T
f (t ) a 0 {a n cos n 0 t bn sin n 0 t}
n 1
(2.1)
ここで，a0, an, bn はフーリエ係数（Fourier coefficient）で，特に， a0 を直流という．ω0(=2
π/T）を基本角周波数， nω0 の成分を第ｎ次高調波という．フーリエ係数は次式で計算され
る．
1
a0
T
T /2
f (t )dt
T / 2
2
an
T
T /2
f (t ) cos n tdt
0
T / 2
2
bn
T
T /2
f (t ) sin n tdt
(2.2)
0
T / 2
式(2.1)の表現で使用される関数の集合 {1, cos 0 t , sin 0 t , , cos n 0 t , sin n 0 t , } は完全直交系
をなす．すなわち，
1
T
1
T
1
cos m 0 t cos n 0 t dt mn
T / 2
2
T /2
1
T / 2 sin m 0 t sin n0 t dt 2 mn
T /2
T /2
T / 2
1 m n
0 m n
mn
sin m 0 t cos n 0 t dt 0
m, n 1, 2,
ここで，完全とはパーセバルの等式(後述の 3)を参照）が成立することを意味する．
次に，複素フーリエ級数展開を求める．オイラー（Euler）の公式
cos n0t
e jn 0 t e jn 0 t
e jn 0 t e jn 0 t
, sin n0t
2
2j
(2.3)
を，式(2.1)に代入し，整理すると
f (t ) a 0
1
2
1
a0
2
{a
n (e
jn0t
e jn0t ) jbn (e jn0t e jn0t )}
n 1
{(a
(2.4)
n
jbn )e
jn0t
(a n jbn )e
jn0t
}
n 1
n≠0 の場合について
信号処理 1 章-4 章.docxx
- 4 -
大西
昇（名古屋大学）
a n jbn
2
cn
(2.5)
とおき，これに式(2.2)の an, bn を代入すると
T /2
T /2
1 T /2
1 2
2
cn
f
(
t
)
cos
n
tdt
j
f (t ) sin n0tdt
f (t )e jn 0 t dt
0
T / 2
2 T T / 2
T T / 2
T
(2.6)
となる．また，
c n c n *
a n jbn
2
(2.7)
である(∵f(t)は実数値関数)．さらに，式(2.6)で n=0 とすると
1
T
c0
T /2
T / 2
f (t )dt a 0
(2.8)
したがって，式(2.4)に式(2.5),(2.7),(2.8)を代入すると，次式の複素フーリエ級数展開が得
られる．
f (t ) c0
1
cn
T
(cn e jn0t c n e jn0t )
n 1
T /2
T / 2
c e
n
n
jn0t
(2.9)
f (t )e jn0t dt
上式の cn を複素フーリエ係数という．f(t)が偶関数なら cn は実数，奇関数なら cn は純虚数とな
る．
式(2.5)，(2.7)より，フーリエ係数と複素フーリエ係数の関係は次式となる．
a n c n c n , b n j (c n c n )
例
1
0
周期 T=1 の矩形関数 f (t )
(2.10)
t 1/ 4
1/ 4 t T / 2
のフーリエ級数展開を求める．
0 2 T 2 であるので，式(2.9)より
1
cn
1
4
1
4
e j 2 nt 4
e jn / 2 e jn / 2
n
1
e jn 2t dt
sin
j 2 n
n
2
j 2 n 1
4
f (t )
n
n
n
sin
1
2 e j 2nt
2 cos 2nt
n
2 n 1 n
2
sin
上式第 2 項の無限級数和を，項数 N=1, 6, 20 で打ち切った場合の波形を図 2.1 に示す．項数
N の増加につれて，矩形に近づくことが分かる．
信号処理 1 章-4 章.docxx
- 5 -
大西
-1
1
1
1
0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
-0.5
0.5
1
-1
(a) N=1
-0.5
0.5
(b) N=6
図 2.1
1
-1
-0.5
昇（名古屋大学）
0.5
1
(c) N=20
矩形関数のフーリエ級数展開
2) スペクトル
スペクトル (spectrum)とは，振動数（周波数）の関数として，信号の特性（強度，位相）を
表示することである．周期関数の f(t)のスペクトルは，フーリエ係数 cn で表される．
c n c n e jn と極座標表現したとき， c n を振幅スペクトル， n を位相スペクトルといい，そ
れぞれ次式のように計算される．
c n a n2 bn2 2
(2.11)
n tan 1 (bn / a n )
(2.12)
3) パーセバルの等式
フーリエ級数展開では，次式のパーセバルの等式(Parseval’s formula)が成立する．この定理
は信号の 2 乗平均値すなわち電力は，各周波数成分の電力の総和に等しいことを意味する．
1
T
T /2
T / 2
f (t ) dt
2
c
n
2
(2.13)
n
4) ギブス現象
フーリエ級数展開では，図 2.2 のように，不
連続点で振動（リップル：さざ波）が生じる．
これをギブス（Gibbs）現象という．この振動
は，項数 N が有限である限り発生し，リップル
の振幅は，N の値に関わらず一定である．ギブ
ス現象は，不連続点での近似誤差と言える．
2.2 フーリエ変換
図 2.2
ギブス現象(○で囲まれた部分)
1) 定義
フーリエ級数展開の周期 T を無限大にすると，フーリエ変換（Fourier transform）が得られる
(詳細は文献(1)を参照)．
F ( )
f (t )e jt dt
(2.14)
信号処理 1 章-4 章.docxx
- 6 -
大西
ただし，f(t)のフーリエ変換 F ( ) 3が存在するのは，f(t)が絶対可積分（
昇（名古屋大学）
f (t ) dt ）の場
合である．なお，有限継続信号は概ねフーリエ変換可能である．
フーリエ逆変換は
f (t )
1
2
F ( )e jt d
となる．フーリエ変換対を
(2.15)
f (t ) F ( )
と記す4．
一般的に複素数となる F(ω)を
F ( ) F ( ) exp[ j ( )]
(2.16)
と，極座標表現した時， F ( ) を振幅スペクトル， ( ) を位相スペクトルと言う．
2) フーリエ変換の性質
f (t ) F ( )， g (t ) G ( ) とする．
a) f(t)が実数値関数なら
b) 線形性
F * ( ) F ( )
(2.17)
af (t ) bg (t ) aF ( ) bG ( )
(2.18)
c) f ( p(t q))
e jq
F( )
p
p
d)
df (t )
jF ( )
dt
e)
t
f)
（微分は高周波成分を強調）
1
F ( )
j
f ( )d
(2.19)
（積分は高周波成分を抑制）
f ( ) g (t )d F ( )G ( )
g) f (t ) g (t )
1
2
(2.20)
F (u )G ( u ) du
(2.21)
（畳み込み積分）
1
F ( ) * G ( )
2
(2.22)
（複素畳み込み積分）
(2.23)
h) パーセバルの等式（Parseval’s formula）(信号のエネルギーについての等式)
2
f (t ) dt
1
2
2
F ( ) d
(2.24)
3) フーリエ変換の例
a)矩形(方形)波
3
1
f (t )
0
t / 2
t / 2
フーリエ変換を F(jω)と記し，ラプラス変換との関係が分かるようにした本もある．
4
ラプラス変換との関係が分かるように，フーリエ変換を F(jω)と記した本や，逆変換での 1/(2
π)といった係数が現われない，周波数 f を用いたフーリエ変換を記した本もある．
信号処理 1 章-4 章.docxx
- 7 -
大西
F ()
e
2
jt
2
jt
j
e
e
dt
j
2
j
e
j
2
2
昇（名古屋大学）
2 sin
2 sinc( / 2)
(2.25)
2
sin(x) / x で，sinc 関数と言う．
ここで， sinc(x)
F ()
0.5
例
1 2
1
の場合， F () sinc( / 4) で図 2.3 のよ
2
0.4
0.3
0.2
うになる．
0.1
-30
-20
超関数であるデルタ関数 (t ) は次の等式を満足する． 図 2.3
20
30
矩形関数のフーリエ変換
(t )dt 1
(2.26)
(t ) 0
t 0
10
-0.1
b)単位インパルス(デルタ関数)
-10
(t a ) (t )dt ( a)
(2.27)
よって， (t ) のフーリエ変換は，式(2.27)を用いると，次式となる．
F ( )
(t )e jt dt e j 0 1
(2.28)
すなわち，デルタ関数は全ての周波数成分を同じパワーで含む信号である．
c)複素正弦波信号
デルタ関数を導入すると，角周波数がω0 の正弦波のフーリエ変換 F[e j0t ] は次式となる．
F [e j0t ]
e j ( 0 )t dt [
ここで， ( x) lim
2 sin( 0 )t 0
e j ( 0 ) t
] lim
2 ( 0 ).
j ( 0 ) t0
0
(2.29)
sin( x )
0
を使用している．デルタ関数の極限としての定義は文献(2)を参
x
照のこと．したがって，変換対は次のようになる．
e j0t 2 ( 0 )
(2.30)
d)指数関数
指数関数の変換対は次のようになる．
e at u (t )
1
j a
a 0
(2.31)
0 t0
1
t0
ここで，u(t)は単位ステップ(階段)関数で， u (t )
4) フーリエ級数展開とフーリエ変換の関係
信号処理 1 章-4 章.docxx
- 8 -
大西
昇（名古屋大学）
フーリエ級数展開とフーリエ変換は，時間関数の周波数領域での表現である．フーリエ級数
展開は，基本周波数の整数倍となる周波数での離散的な表現（離散スペクトル）であるのに対
し，フーリエ変換は連続的な表現である．
有限継続信号 f(t)を，その継続時間よりも長い周期 T で拡張した周期信号のフーリエ係数
cn は，f(t)のフーリエスペクトル F(ω)のω0(=2π/T）毎の標本値 F (n 0 ) (n 0, 1, 2, ) と 1
対 1 に対応する(図 2.4 参照)．すなわち，次式が成立する．
F (n 0 ) Tc n
n 0, 1, 2,
(2.32)
ここで，cn は各周波数 nω0 の電力で，F(ω)はエネルギー密度である．
上記のことは，図 2.5(a)のように，有限継続信号 f(t)を周期Ｔで拡張することは，f(t)の
フーリエスペクトルを周波数軸上で２π/Ｔ毎にサンプリングすることと等価である．これに
双対5な関係として，有限継続スペクトル F(ω)を周期 W で拡張することは，時間領域で２π/
Ｗ毎にサンプリングすることと等価である（図(b)参照）
F()
図 2.4
フーリエ係数（赤の線, 値は Tcn）とフーリエ変換（実線）との関係
f(t)
t
ω
T
2π/T
(a) 時間領域での周期的拡張
F(ω)
t
ω
W
2π/W
(b) 周波数領域での周期的拡張
図 2.5
5
時間・周波数領域での双対性
双対（duality） 二重性，相似関係．
信号処理 1 章-4 章.docxx
- 9 -
大西
昇（名古屋大学）
2.3 ラプラス変換
1) 定義
ラプラス変換(Laplace transform)は次式で定義される．
F s
0
s j
f (t ) exp( st )dt
(2.33)
これは，f (t )e t u (t ) のフーリエ変換に相当する．したがって，ラプラス変換は，
f (t ) dt
を満足しない関数にも適用できる．式(2.33)の積分が収束する領域σ＞σa を収束領域と言
う．ラプラス変換では，その収束領域の明記が必要である．
ラプラス逆変換は
f (t )
1
2j
j
j
(2.34)
F ( s ) exp(st )ds
で，その積分路は収束領域にとる．ラプラス変換対を f (t ) F ( s) と記す．
2) ラプラス変換の性質
f (t ) F ( s )， g (t ) G ( s ) とする．
af (t ) bg (t ) aF ( s ) bG ( s )
a) 線形性
b) f ( p(t q))
c)
d)
p 0, q 0
(2.36)
df (t )
sF ( s ) f (0)
dt
t
t
f ( )d
0
e)
e sq
F( )
p
p
(2.37)
1
F ( s)
s
(2.38)
f ( ) g (t ) d F ( s )G ( s )
0
(2.35)
（畳み込み積分）
f) e at f (t ) F ( s a)
(2.39)
(2.40)
3) ラプラス変換の例
a)単位ステップ関数
1 t 0
u (t )
0 t 0
u (t ) 1
収束領域は Re( s ) 0
(2.41)
収束領域は Re( s ) Re( a )
(2.42)
b) e at 1
s
sa
導出は次の通りである．
0
e
at
e
st
dt
0
e
( s a )t
e ( s a )t
e ( s a ) 1
dt
(s a)
(s a) 0
信号処理 1 章-4 章.docxx
- 10 -
大西
昇（名古屋大学）
Re( s ) Re( a ) (収束領域)なら，分子の第 1 項はゼロとなり，積分が収束する．
c) t 1
収束領域は Re( s ) 0
s2
d) cos at s
s2 a2
(2.43)
収束領域は Re( s ) 0
(2.44)
導出は次の通りである．
cos at Re(e jat ) であるので，
0
cos ate st dt
0
e ( s ja ) 1
Re(e iat )e st dt Re e ( s ja )t dt Re
0
( s ja)
Re( s ) 0 (収束領域)なら，分子の第 1 項はゼロとなるので，
0
1
s
s ja
Re 2
cos ate st dt Re
2
2
2
s a s a
s ja
e) sin at a
収束領域は Re( s ) 0
(2.45)
収束領域は s 平面の全領域
(2.46)
s2 a2
f) (t ) e s
4) ラプラス逆変換の求め方
X(s)のラプラス逆変換である
x(t )
j
2 j
1
j
X ( s ) exp(st )ds
を求める方法として，部分分数展開による方法と留数計算による方法とを説明する．
a) 部分分数展開
X(s)が実係数有理関数，すなわち
X ( s)
q 0 s m q1 s m 1 q m
s n p1 s n 1 p n
（m≦n）
(2.47)
である場合，部分分数展開で逆変換は得られる．分母の多項式を複素根まで許して因数分解
する．
(i)すべて単根λi(i=1,2,...,n)の時
分母の多項式が次式のように因数分解されるとする（ i (i 1, , n) は全て異なる）と，
s n p1 s n 1 p n ( s 1 )( s 2 ) ( s n )
(2.48)
X(s)は部分分数展開で次式となる．
X ( s) a0
an
a1
a2
s 1 s 2
s n
(2.49)
上式の係数は次式で求まる．
a 0 X (),
a i ( s i ) X ( s ) s
i
i 1, , n
(2.50)
信号処理 1 章-4 章.docxx
- 11 -
大西
昇（名古屋大学）
1
] exp(t )u (t ) であるので，逆変換は次式で与えられる．
s
x(t ) a 0 (t ) {a1 exp(1t ) a 2 exp( 2 t ) a n exp( n t )}u (t )
(2.51)
L1 [1] (t ) ， L1 [
(ii)k 重根αをもち，他は全て単根λi(i=1,2,...,n-k)の時
分母の多項式が ( s ) k を因数に持ち，他は全て異なる単根とする．すなわち，
s n p1 s n 1 p n ( s 1 )( s 2 ) ( s n k )( s ) k
(2.52)
X(s)は部分分数展開で
X (s) a0
an k
bk
a1
a2
b1
b2
2
s 1 s 2
s n k
s
(s )
(s ) k
(2.53)
ここで，係数 bi は次式で求まる．
1
d k i
[( s ) k X ( s)] s
(k i )! ds k i
bi
L1 [
n!
(s )
n 1
] t n exp( t )u (t )
(2.54)
を使用すると，逆変換は次式で与えられる．
x(t ) a 0 (t ) {a1 exp(1t ) a 2 exp( 2 t ) a n k exp( n k t )
t k 1
b1 exp(t ) b2 t exp(t ) bk
exp(t )}u (t )
(k 1)!
(2.55)
b) 留数計算
留数計算による逆変換の求め方を，図 2.6 を用いて説明する(詳細は文献(3)を参照)．
lim X ( s) 0 なら，
R
1
lim
R 2j
X ( s) exp( st )ds
0
t 0
(2.56)
X ( s) exp( st )ds
0
t 0
(2.57)
C2
1
lim
R 2j
C1
となる．ここで，C1 は点 BDA を通る円弧で，C2 は点 AEB を通る円弧である．
まず，t<0 の場合，式(2.56)を利用すると，次式が得られる．
x(t )
1
2j
j
1
lim
R 2j
j
X ( s ) exp( st )ds lim {
R
1
2j
1
X (s) exp(st )ds 2j X (s) exp(st )ds}
A B
C2
(2.58)
X (s) exp(st )ds 0
C'
上式で，閉路 C’は直線 AσB と円弧 C2 で作られる積分路で，絶対収束領域内にあり，X(s)
は正則であるので，一周積分はゼロとなる．
次に，t>0 の場合，式(2.57)を利用すると，次式が得られる．
信号処理 1 章-4 章.docxx
- 12 -
大西
x(t )
1
2j
lim
R
j
1
2j
j
X ( s) exp( st )ds lim {
R
1
2j
昇（名古屋大学）
X (s) exp(st )ds 2j X (s) exp(st )ds}
1
A B
(2.59)
C1
n
X (s) exp(st )ds A
k
k 1
C
ここで，
閉路 C は直線 AσB と円弧 C1 で作られる積分路，
n は積分路 C の内部にある極の数で，Ak はｋ番目の極 ak
の留数である．留数とは，X(s)exp(st)をローラン展開
した時の，1/(s-ak)の係数で，
ak が１次の極ならば，
Ak ( s a k ) X ( s ) exp( st )
(2.60)
s ak
ak がｎ次の極ならば，
d n 1
1
Ak
( s a k ) n X ( s) exp(st )
n 1
(n 1)! ds
s ak
(2.61)
である．
図 2.6
ラプラス変換を求める
ときの積分路
1
X (s)
( s a )( s b) 2
例
Re( s ) a b
の逆
変換を求める．
解法１（部分分数展開）
X ( s)
1 /(a b) 2 1 /(b a) 2 1 /(b a)
sa
sb
( s b) 2
より
1
1
1
t exp(bt )
exp( at )
exp(bt )
2
2
x (t ) ( a b )
(b a )
(b a )
0
(t 0)
(t 0)
解法２（留数計算）
極 s=a の留数は，式(2.60)を用いて，
exp( st )
( s b) 2
sa
1
exp( at )
( a b) 2
(t 0)
極 s=ｂの留数は，式(2.61)を用いて，
t exp( st )( s a) exp( st )
( s a) 2
1
1
t exp(bt )
exp(bt )
(b a)
(b a ) 2
d exp( st )
ds ( s a )
s b
s b
(t 0)
よって，解法１と同じ x(t)が得られる．
信号処理 1 章-4 章.docxx
- 13 -
大西
3
昇（名古屋大学）
不規則信号
気温や音声のように，値の時間的な変化が確定的ではなく，確率的に変動する時間信号が身近
には多い．本章では，値の変化が不確定である不規則信号の記述・解析に関係する特性関数，キ
ュムラント関数，相関関数および，不規則信号としての望ましい性質である定常性とエルゴード
性を説明する．
3.1 不規則信号の記述
1.6 で述べたように，信号は
確定信号と不規則信号とに分
類される．不規則信号は，値
程（stochastic process, random
process）とも言う．
．．．．．
則に変動する信号で，確率過
．．．．．
が確率的な法則に従って不規
x (1) (t )
x ( i ) (t )
(i )
番目に観察された信号を
．．．．．
に，観察のたびに異なる． i
x (t )
．．．．．
不規則信号は，図 3.1 のよう
x ( i ) (t ) と 表 し ， 標 本 信 号
(sample signal)と言う．不規則
合
．．．
集
．．．
信号は，標本信号の無限個の
x ( n ) (t )
(i )
{x (t ), x (t ), x (t ), }
(1)
( 2)
のことで，この集合を信号確
図 3.1
X (t1 )
信号確率集合としての不規則信号
率集合と言い，{x(t)}あるい
は X(t)と表す．本書では X(t)を使用する．
信号確率集合である不規則信号の統計的性質を，確率密度関数で記述する．例えば，n 通りの
時点での， X (t1 ) x1 , X (t 2 ) x 2 , , X (t n ) x n という，信号値の同時生起は，次式の n 次元
確率密度関数で記述される．
p n ( x1 , x 2 , , x n ; t1 , t 2 , , t n )
lim
xi ( i 1, , n ) 0
Pr x1 X (t1 ) x1 x1 , x 2 X (t 2 ) x 2 x 2 , , x n X (t n ) x n x n
x1 x 2 x n
(3.1)
ただし，不規則信号を厳密に記述するには，
n 1,2, ,
p n ( x1 , x 2 , , x n ; t1 , t 2 , , t n )
(3.2)
と，無限次元確率密度関数を必要とする．
3.2 集合平均と時間平均
集合平均(ensemble average)とは，集合の標本関数についての平均である(図 3.1 参照)．例え
信号処理 1 章-4 章.docxx
- 14 -
大西
昇（名古屋大学）
ば，時刻 t1 での平均（1 次モーメント関数）は次式となる．
N
1
N
m1 (t1 ) X (t1 ) lim
N
x
(n)
n 1
(t1 )
(3.3)
一般的に，n 通りの時点でのｎ次モーメント関数は次式となる．
m n (t1 , t 2 , , t n ) X (t1 ) X (t 2 ) X (t n ) lim
N
1
N
N
x
(k )
k 1
(t1 ) x ( k ) (t 2 ) x ( k ) (t n )
(3.4)
集合平均は，確率集合の密度関数を使った期待値演算で表現できる．
X (t1 ) E[ X (t1 )]
x1 p1 ( x1 ; t1 )dx1
X (t1 ) X (t 2 ) E[ X (t1 ) X (t 2 )]
(3.5)
x1 x 2 p 2 ( x1 , x 2 ; t1 , t 2 )dx1 dx2
(3.6)
上式で， p 2 ( x1 , x 2 ; t1 , t 2 ) は，X(t1)=x1 かつ X(t2)=x2 である確率密度である．
時間平均(time average)は，時間軸についての平均である．まず，標本信号 x ( n ) (t ) の時間平
均値は次式となる．
1
T 2T
x ( n ) (t ) lim
T
T
x ( n ) (t )dt
(3.7)
同一標本信号の 2 時点での値の積の時間平均である自己相関関数(3.7 で詳述)は
1
T 2T
x ( n ) (t ) x ( n ) (t ) lim
T
T
x ( n ) (t ) x ( n ) (t )dt
(3.8)
で計算される．２つの不規則信号の相関を調べるために使用される相互相関関数(3.7 で詳述)
は，次式となる．
x ( n ) (t ) ( n ) y (t )
1
lim
T 2T
T
T
x ( n ) (t ) y ( n ) (t )dt
(3.9)
3.3 特性関数
確率密度関数のフーリエ変換を特性関数(characteristic function)という．1 次元確率密度関数
のフーリエ変換を，１次元特性関数と言い
1 (u1 ; t1 )
p1 ( x1 ; t1 )e ju1x1 dx1 E[e ju1 X (t1 ) ]
(3.10)
で表される．一般的に，n 次元特性関数は次式となる．
n (u1 , u 2 , , u n ; t1 , t 2 , , t n )
p n ( x1 , x 2 , , x n ; t1 , t 2 , , t n )e j ( u1x1 u2 x2 un xn ) dx1 dx 2 dx n
(3.11)
E[e j{u1 X ( t1 ) u2 X (t 2 ) un X ( tn )} ]
1 次元特性関数とモーメント関数6には，つぎの関係が成立する．
6
不規則信号のモーメントは，時間の関数であるのでモーメント関数と言う（キュムラント関数も同様）．
信号処理 1 章-4 章.docxx
- 15 -
大西
1 (u1 ; t1 ) E[e ju1 X ( t1 ) ]
d1 (u1 ; t1 )
E[ jX (t1 )]
du1
u 0
より
昇（名古屋大学）
だから
1
m1 (t1 )
1 d1 (u1 ; t1 )
j
du1
u 0
(3.12)
1
同様にして
d 21 (u1; t1 )
2
du1
u
E[ j 2 X (t1 ) 2 ]
だから
1 0
m2 (t1 , t1 )
1 d 2 1 (u1 ; t1 )
2
j2
du1
u 0
(3.13)
1
一般的には，次式のように，モーメント関数が特性関数から計算される．
mn (t1 , t1 ,, t1 )
1 d n 1 (u1 ; t1 )
n
jn
du1
u 0
(3.14)
1
特 性 関 数 を 原 点 u1=0 の 近 傍 で テ イ ラ ー 級 数 展 開 ， す な わ ち マ ク ロ ー リ ン 級 数 展 開
（ f (t ) f (0)
n 1
f ( n ) (0) n
t ）すると
n!
( ju1 ) n
1 (u1 ; t1 ) 1 mn (t1 ,t1 ,, t1 )
n!
n 1
(3.15)
この式から，モーメント関数を与えると特性関数，すなわち確率密度関数が確定することが
分かる．
次に，２次元特性関数は
2 (u1 , u 2 ; t1 , t 2 ) E[e j{u1 X ( t1 ) u2 X ( t2 )} ]
であるので，２次元特性関数とモーメント関数とには，次の関係が成立する．
m1 (t1 )
1 2 (u1 , u 2 ; t1 , t 2 )
j
u1
u u
1
m1 (t 2 )
(3.16)
2 0
1 2 (u1 , u 2 ; t1 , t 2 )
j
u 2
u u
1
m2 (t1 , t1 )
(3.17)
2 0
1 2 2 (u1 , u 2 ; t1 , t 2 )
2
j2
u1
u u
1
m 2 (t1 , t 2 )
(3.18)
2 0
1 2 2 (u1 , u 2 ; t1 , t 2 )
u1u 2
j2
u u
1
(3.19)
2 0
信号処理 1 章-4 章.docxx
- 16 -
大西
m 2 (t 2 , t1 )
1 2 2 (u1 , u 2 ; t1 , t 2 )
u 2 u1
j2
u u
1
m2 (t 2 , t 2 )
(3.20)
2 0
1 2 2 (u1 , u 2 ; t1 , t 2 )
2
j2
u 2
u u
1
昇（名古屋大学）
(3.21)
2 0
２次元特性関数のマクローリン級数展開は，上記のモーメント関数を用いると，次式のよう
に表現される．
ju1
ju
m1 (t 2 ) 2 }
1!
1!
2
2
j u1u 2
j 2 u1u 2
( ju1 )
( ju 2 ) 2
{m2 (t1 , t1 )
m2 (t1 , t 2 )
m2 (t 2 , t1 )
m2 (t 2 , t 2 )
}
2!
2!
2!
2!
2 (u1 , u 2 ; t1 , t 2 ) 1 {m1 (t1 )
(3.22)
3.4 キュムラント関数
キュムラント関数（cumulant function）は次式で定義される．
1 次キュムラント関数
k1 (t1 )
1 d ln 1 (u1 ; t1 )
j
du1
u 0
(3.23)
1
n 次キュムラント関数
k n (t1 , t1 ,, t1 )
1 d n ln 1 (u1 ; t1 )
n
jn
du1
u 0
(3.24)
1
上記のキュムラント関数は，１次元特性関数の自然対数をマクローリン級数展開した時の係
数で
ln 1 (u1 ; t1 ) k n (t1 ,t1 ,, t1 )
n 1
( ju1 ) n
n!
(3.25)
となる．同様に，２次元特性関数の自然対数のマクローリン級数展開は
ju1
ju
k1 (t 2 ) 2 }
1!
1!
2
2
( ju1 )
( ju 2 ) 2
j u1u 2
j 2 u1u 2
{k 2 (t1 , t1 )
k 2 (t1 , t 2 )
k 2 (t 2 , t1 )
k 2 (t 2 , t 2 )
}
2!
2!
2!
2!
ln 2 (u1 , u 2 ; t1 , t 2 ) {k1 (t1 )
(3.26)
式(3.25), (3.26)で，展開の係数がキュムラント関数となっている．
キュムラント関数とモーメント関数も同じ特性関数から計算され，それらには関係が成立す
る．まず， １次元特性関数の場合，
k1 (t1 ) m1 (t1 )
k 2 (t1 , t1 ) m 2 (t1 , t1 ) m1 (t1 )
(3.27)
2
(3.28)
k 3 (t1 , t1 , t1 ) m3 (t1 , t1 , t1 ) 3m1 (t1 )m 2 (t1 , t1 ) 2m1 (t1 ) 3
(3.29)
k 4 (t1 , t1 , t1 , t1 ) m 4 (t1 , t1 , t1 , t1 ) 3m2 (t1 , t1 ) 4m3 (t1 , t1 , t1 )m1 (t1 )
2
12m 2 (t1 , t1 ) m1 (t1 ) 2 6m1 (t1 ) 4
(3.30)
信号処理 1 章-4 章.docxx
- 17 -
大西
昇（名古屋大学）
次に，２次元特性関数の場合，
k1 (t1 ) m1 (t1 ) ， k1 (t 2 ) m1 (t 2 )
(3.31)
k 2 (t1 , t1 ) m 2 (t1 , t1 ) m1 (t1 ) ， k 2 (t1 , t 2 ) m 2 (t1 , t 2 ) m1 (t1 )m1 (t 2 ) ，
2
k 2 (t 2 , t1 ) m 2 (t 2 , t1 ) m1 (t1 )m1 (t 2 ) ， k 2 (t 2 , t 2 ) m 2 (t 2 , t 2 ) m1 (t 2 ) 2
(3.32)
キュムラント関数は，つぎの性質を有する．下記で， k n ( X ) は n 次キュムラントである．
性質１
X, Y が統計的に独立なら， n
すなわち， n
k n ( X Y ) k n ( X ) k n (Y )
(3.33)
kn ( X , Y ) 0
k n ( X ) n k n ( X )
性質２
n
性質３
X が多次元正規分布に従うなら，3 次以上の全てのキュムラントは全てゼロとなる．
β：定数
(3.34)
性質１と３は，モーメント関数にはない．
3.5 ガウス性不規則信号
３次以上のキュムラント関数がすべてゼロである不規則信号を，ガウス性不規則信号 7
（Gaussian random signal）と言う．
簡単な 1 次元特性関数を考える．1 次キュムラント関数 k1 (t i ) および 2 次のキュムラント関数
k 2 (t i , t i ) を与えれば，式(3.25)より
ln 1 (u1 ; t1 ) 0 k1 (t1 )( ju1 ) k 2 (t1 , t1 )
( ju1 ) 2
2!
2
( ju1 ) 2
u1
1 (u1 ; t1 ) expk1 (t1 )( ju1 ) k 2 (t1 , t1 )
expk1 (t1 )( ju1 ) k 2 (t1 , t1 )
2!
2
(3.35)
これをフーリエ逆変換すれば，１次元確率密度関数 p1 ( x1 ; t1 ) は
1
p1 ( x1 ; t1 )
2
1
2
e
1 (u1 ; t1 )e
k 2 ( t1 ,t1 )
u12
2
e
jx1u1
j ( x1 k1 ( t1 )) u1
1
du1
2
du1
e
k1 ( t1 )( ju1 ) k 2 ( t1 ,t1 )
u12
2
e jx1u1 du1
(3.36)
1
2 k 2 (t1 , t1 )
e
( x1 k1 ( t1 )) 2 2 k 2 ( t1 ,t1 )
と，1 次元正規分布の式が得られる．一般に，ガウス性不規則信号のｎ次元の確率密度関数は，
次式で与えられる．
p n ( x1 , , x n ; t1 , , t n )
1
exp ( x1 k1 (t1 ), , x n k1 (t n )) K 21 ( x1 k1 (t1 ), , x n k1 (t n )) T
2
(2 ) n det K 2
1
(3.37)
k 2 (t1 , t1 ) k 2 (t1 , t n )
ここで， K 2 [ k 2 (t i , t j )]
k 2 (t n , t1 ) k 2 (t n , t n )
7
正規不規則信号（normal random signal）とも言う．
信号処理 1 章-4 章.docxx
- 18 -
大西
昇（名古屋大学）
3.6 定常性とエルゴード性
定常(stationary)とは，不規則信号の統計的性質が時刻ｔによって変化しないことである．す
なわち，定常信号では
(3.38)
p1 ( x1; t1 ) p1 ( x1 )
と，確率密度関数が時刻によらず一定で，
p 2 ( x1 , x 2 ; t1 , t 2 ) p n ( x1 , x 2 ; 1 ) 1 t 2 t1
(3.39)
p n ( x1 , x 2 , , x n ; t1 , t 2 , , t n ) p n ( x1 , x 2 , , x n ; 1 , 2 , , t n 1 ) i t i 1 t i
(3.40)
と，全ての同時確率密度関数が時間差τで表現される．これは狭義の定常(強定常)である．条
件を緩め，式(3.38)，(3.39)および E[ X 2 (t )] を満足する場合を，広義の定常(弱定常)と言
う．
定常不規則信号の特徴の一つは，時間平均が意味を持つことである．不規則信号の標本の時
間平均が互いに一致し，集合平均に等しいことをエルゴード性(ergodic)と言い，その性質を有
する不規則信号をエルゴード信号と言う． 例えば，１次モーメントの場合，
x (1) (t ) x ( 2) (t ) x ( n ) (t ) .
x(t ) X (t )
となることである．
定常エルゴード不規則信号では，その平均について次式が成立する．
X (t ) x(t ) lim
T
1
2T
T
T
x(t ) dt
(3.41)
3.7 相関関数とスペクトル
本節で，不規則信号は，定常かつエルゴード性で，平均値はゼロとする．
1) 自己相関関数
自己相関関数(auto-correlation function)は次式のように，τだけ離れた 2 時刻の信号の値に
どれだけの相関があるかを表す．
1
T 2T
xx ( ) x(t ) x(t ) lim
T
T
x(t ) x(t )dt
(3.42)
自己相関関数は，2 次キュムラント関数と一致する．すなわち，
k 2 (t1 , t 2 ) k 2 (t 2 t1 ) k 2 ( ) xx ( ), t 2 t1
(3.43)
自己相関関数は，つぎの性質を持つ．
a） xx (0) x(t ) 2
b） xx ( ) xx ( )
（信号の分散（パワー））
．
（偶関数）
(3.44)
次式を規格（正規）化自己相関関数と言う．
xx ( ) xx ( ) / xx (0).
xx ( ) 1
(3.45)
2) 相互相関関数
信号処理 1 章-4 章.docxx
- 19 -
大西
昇（名古屋大学）
相互相関関数(cross-correlation function)は次式のように，異なる信号 x (t ), y (t ) で，τだけ離
れた 2 時刻での信号値にどれだけの相関があるかを表す．
xy ( ) x(t ) y (t ) lim
T
1
2T
T
T
x (t ) y (t ) dt
(3.46)
相互相関関数については次の関係が成立する(添え字の順番が両辺で異なることに注意)．
xy ( ) yx ( )
(3.47)
規格化相互相関関数は次式となる．
xy ( ) xy ( ) / xx (0) yy (0) .
xy ( ) 1
(3.48)
3) パワースペクトル
不規則信号 X(t)のパワースペクトル xx ( ) は次式で定義される．
X T ( )
xx ( ) lim
2
2T
T
2
(n)
X
(
)
T
lim
T
2T
n 1
N
1
lim
N N
(3.49)
ここで，
X T( n ) ( )
xT( n ) (t )e jt dt ,
x ( n ) (t )
xT( n ) (t )
0
T t T
t T
(3.50)
パワースペクトル xx ( ) と自己相関関数 xx ( ) とには，次のフーリエ変換対の関係（ウィ
ナー・ヒンチンの定理, Wiener-Khintchine’s Theorem）が成立する(証明は文献(6)を参照)．
xx ( )
xx ( )e j d
xx ( )
1
2
xx ( )e j d
(3.51)
パワースペクトルはつぎの性質を有する．
a） xx ( ) xx ( )
(3.52)
b） xx ( ) 0
(3.53)
次に，不規則信号 X(t), Y(t)のクロススペクトル xy ( ) は次式となる．
X T* (ω)YT (ω)
1
xy ( ) lim
lim
T
N
2T
N
*
X T( n ) (ω)YT( n )(ω)
Tlim
2T
n 1
N
(3.54)
ここで，
X
( n)
T
( )
x
(n)
T
(t )e
j t
dt ,
x
(n)
T
x ( n ) (t )
(t )
0
T t T
t T
信号処理 1 章-4 章.docxx
- 20 -
大西
(n)
T
Y
( )
y
(n)
T
(t )e
j t
dt ,
y
(n)
T
y ( n ) (t )
(t )
0
昇（名古屋大学）
T t T
t T
クロススペクトル xy ( ) と相互相関関数 xy ( ) とには，次式の関係がある．
xy ( )
xy ( )e j d
xy ( )
1
2
xy ( )e j d
(3.55)
クロススペクトルはつぎの性質を有する．
a） xy ( ) xy ( )
(3.56)
b） yx ( ) xy ( )
(3.57)
*
*
3.8 不規則信号と線形システム
図 3.2 のように，入力が x(t) ，出力が y(t)
の線形システムがある．システムのインパルス
線形システム
応答を h(t)とすると，出力は
y (t ) h( ) x(t )d
x(t)
y(t)
h(t)
(3.58)
と，x(t) と h(t)の畳み込み積分となる（x(t) ，
図 3.2
線形システム
y(t)， h(t)はすべて実数値関数とする）．入力 x(t)が定常エルゴード性信号であると，出力
y(t)もまた定常エルゴード性信号である．システムの周波数特性を与えるインパルス応答の
フーリエ変換 H ( ) と，パワースペクトルおよびクロススペクトルとには，次の関係が成立
する．
xy ( ) H ( ) xx ( )
(3.59)
yy ( ) H ( ) yx ( )
(3.60)
yy ( ) H ( ) xx ( )
(3.61)
2
すなわち，システムの周波数特性 H ( ) は，式(3.59)もしくは(3.60)で求まる．
信号処理 1 章-4 章.docxx
- 21 -
大西
4
昇（名古屋大学）
信号のディジタル化
コンピュータで信号処理するためには，アナログ信号のディジタル信号への変換，すなわち
標本化と量子化が必要である．まず，標本化と量子化を説明し，標本化において重要な標本化
定理，そして標本化周波数の変換を説明する．さらに，ディジタル信号からアナログ信号への
変換も簡潔に説明する．
4.1 標本化と量子化
アナログ信号をディジタル信号にするためには，図 4.1 のように，時間軸と振幅軸の両方で
の離散化が必要である．時間軸での離散化，すなわち離散的な時点での信号の値を取ることを
標本化と言う．一方，振幅軸での離散化とは，信号の値を有限語長(桁数)で表現すること（信
号値の離散化）で，量子化と言う．図のように，標本化と量子化の両方を行うことを A/D 変換
と言う．標本化と量子化は誤差を発生することがある．前者での誤差はエイリアシング，後者
での誤差は打切り/丸め誤差である．
標
量
本
子
化
化
アナログ信号
ディジタル信号
離散時間信号
A/D 変 換
標本化誤差
量子化誤差
図 4.1 信号のディジタル化
4.2 標本化定理
標本化では，時間間隔 T の整数倍での時点での信号値をとる．時間間隔 T の逆数 1/T を標本
化周波数(sampling frequency)fs，2π/T を標本化角周波数ωs という．標本化の時間間隔 T に関
する重要な定理を，まず説明する．
定理 Ｘ(ω)＝０
｜ω｜≧ωmax というように，帯域制限された信号ｘ(ｔ)を，式(4.1)を
満足する時間間隔Ｔで標本化するなら，
Ｔ ≦ π/ωmax （ωmax / πをナイキスト周波数という）
(4.1)
式(4.2)により，標本値ｘ(ｎＴ)からｘ(ｔ)で完全に復元できる．
x(t )
x(nT )
n
sin
T
T
(t nT )
(4.2)
(t nT )
以下で，上記の標本化定理を導出する．信号 x(t)を時間間隔 T（ T / max ; s 2 max ）
で標本化する．得られた標本値系列 { x ( nT )} を次式で表わす．
信号処理 1 章-4 章.docxx
- 22 -
大西
x s (t )
昇（名古屋大学）
x(nT ) (t nT ) x(t )
n
ここで， T (t )
T
(t )
(4.3)
(t nT ) で，周期 T のインパルス列である． x (t ) は連続時間信号である
s
n
ので， x s (t ) x(t ) T (t ) をフーリエ変換（式(2.23)の複素畳み込み積分を参照）すると，次
式となる．
X s ( )
1
X ( ) * T ( )
2
(4.4)
上式の T ( ) を求める． T (t ) は周期 T の周期関数であるので，フーリエ級数展開でき
T (t )
1
T
e
jn s t
s 2 / T
n
( c n
1
T
T /2
(t )e jn t dt
s
T / 2
1 jns 0
e
)
T
(4.5)
となる．これのフーリエ変換は，複素正弦波のフーリエ変換の式(2.29)と s 2 / T より，
次式となる．
T ( )
T (t )e jt dt
1
T
n
( n
e jn st e jt dt s
n
s
)
(4.6)
上式を式(4.4)に代入すると，次式を得る．
X s ( )
s
2
1
1
X ( ) * T ( )
2
2
[ s
(u n
n
1
(u n s )X ( u )du
T
n
s
)]X ( u )du
(4.7)
X ( n
n
s
)
式(4.7)は，信号 x(t)のフーリエ変換である X(ω)を周期ωs で拡張し，それらの総和を T で
除したもの
1
T
X ( n
n
s
) が，Xs(ω)であることを意味する（図 4.2(b)参照）．
つぎに，標本値系列 { x ( nT )} から連続時間信号 x(t)を復元する．図 4.2(c)のように，低
域通過フィルタ H(ω)を Xs(ω)にかけて，ωが[-ωmax，ωmax]の範囲を切り出す．
X ( ) H ( ) X s ( )
(4.8)
ここで， H(ω)は遮断角周波数ωc=ωmax の理想低域通過フィルタで，
T
H ( )
0
max
max
(4.9)
式(4.8)のフーリエ逆変換は，h(t)と x s (t ) の畳み込みであり（式(2.22)参照）
，H(ω)のフー
リエ逆変換は h(t )
sin max t
T
2 max
であるので，x(t)が次式のように得られる．
max t
2
x(t ) h(t ) * x s (t ) [ x(nT ) (t nT )] * [
n
sin max t
]
max t
sin max (t )
x(nT ) ( nT )
d
max (t )
n
sin max (t nT )
x(nT )
max (t nT )
n
(4.9)
すなわち，時刻ｔでの信号値 x(t)は，標本値 x ( nT ) ( n , , ) の重み付き和として求
まる．実際には，信号の観察時間は有限区間であるため，式(4.9)のように，無限個の標本
信号処理 1 章-4 章.docxx
- 23 -
大西
昇（名古屋大学）
値が得られない．観察時間が有限な信号 x(t)の復元（外挿問題）や，標本化の時間間隔が一
定でない場合については，文献 (7)を参照のこと．
Ｘ（ω）
(a)
ω
ωmax
-ωmax
TXs（ω）
(b)
ω
ωｓ
(c)
H（ω）低域通過フィルタ
ω
Ｘs（ω）H（ω）
(d)
ω
図 4.2
標本化定理の説明
4.3 エイリアシング
標本化周波数がωs≦2ωmax と，標本化定理が満足されない場合，信号の高周波成分が低周波
領域へ，折り返すことで，スペクトルに歪みが発生する．図 4.3 で，朱色の部分がスペクト
ルの折り返しで，その結果，青の実線のようなスペクトルが得られる．このようなスペクト
ルの変形をエイリアシング（aliasing）8と言う．これは不適切な標本化による誤差（スペクト
ル歪み）である．
Ｘs（ω）
ω
ωmax
ωs
図 4.3
エイリアシングの発生
4.4 信号の量子化
量子化とは振幅値を有限語長(桁数)で表現すること（振幅値の離散化）で，そのために，
打ち切り(切り捨て)，または丸め(通常は四捨五入)という操作がなされる．このような量子
8
エリアシング，折り返し歪み，異名現象とも言う．
信号処理 1 章-4 章.docxx
- 24 -
大西
昇（名古屋大学）
化による誤差を量子化誤差という．
1) 数値表現
数値 x の表現には，単極性（0≦x≦R）と両極性（-R≦x≦R），固定小数点表現と浮動小数
点表現がある．まず，数値ｘの両極性の固定小数点表現を次式とする．
Q[ x] b0 ．
b1 b2 bn ,
bi {0,1} (i 0,1, , n)
(4.10)
a) 符号・絶対値表示
b0=0：正の値，b0=1：負の値
表現できる値の範囲：–(1-2-n)～(1-2-n)．数値の刻み(量子化ステップ)： 2-n．
b) １の補数表示
正の数の場合は，符号・絶対値表示と同じ.
負の数(-X)の場合は，X の符号・絶対値表示の各ビットを反転した表示となる（2-2-n から X
を引いた数の表示）
．
表現できる値の範囲：–(1-2-n)～(1-2-n) ．刻み：2-n．
c) ２の補数表示
正の数 X の場合は，符号・絶対値表示と同じ.
負の数(-X)の場合は，X の１の補数表示の最下位に１を加算した表示となる.（２から X を
引いた数の表示）．
表現できる値の範囲：–1～(1-2-n) ．刻み：2-n．
表 4.1 に，n=2 での，数ｘの固定小数点表現を示す．小数点以下の桁数は 2 であるので，
数 1/8 や-1/8 では，小数点以下３桁目のビットは，打ち切りまたは丸めがなされる．なお，
符号・絶対値表示と１の補数表示とでは，x=0 は 2 通りの表現が存在する．
表 4.1
x
数の固定小数点表現（n=2 の場合）
符号・絶対値表示
１の補数表示
２の補数表示
3/4
0.11
0.11
0.11
1/2
0.10
0.10
0.10
1/4
0.01
0.01
0.01
1/8
0.001
0.001
0.001
0
0.00
-1/8
1.001
1.110
1.111
-1/4
1.01
1.10
1.11
-1/2
1.10
1.01
1.10
-3/4
1.11
1.00
1.01
(1.00)
0.00
-1
(1.11)
0.00
1.00
d) 浮動小数点表現
正の数 F はつぎのように表現される．
F M 2E
(4.11)
信号処理 1 章-4 章.docxx
- 25 -
大西
昇（名古屋大学）
ここで，M は仮数部で， M x 0 x1 x n ，E は指数部で， E c 0 c1 c m ．x0, c0 は符号ビットであ
る． 0.5 M 1
である時，正規化された浮動小数点表示と言う．
2) 量子化誤差
量子化後の値 Q[x]と元の値 x との差を，量子化誤差εという．すなわち，
ε= Q[x] – x
以下で，固定小数点表現での量子化誤差をみる．
丸め操作(n+1 桁が 1 なら n 桁に 1 を足す)の場合，
符号絶対値表示
１の補数，２の補数表示
1 n
1
2 2 n ,
2
2
1 n
1
2 2 n ,
2
2
(4.12)
打切り操作(n+1 桁の値にかかわらず切捨てる)の場合
・ 符号・絶対値表示と１の補数表示では
2 n 0 ( x 0)
0 2 n
(4.13)
( x 0)
・ ２の補数表示
2 n 0
(4.14)
となる．量子化誤差についての詳細は文献(8)を参照のこと．
4.5 標本化周波数の変換
システム内に異なる標本化周波数で動作する処理ブロックが存在するマルチレートシステ
ムが存在する．例えば，ディジタルオーディオ機器（fs=48kHz）とコンパクトディスク
（fs=44.1kHz）からなるシステムである．マルチレートシステムでは，標本化周波数が異なる
処理ブロック間で，信号をやり取りするために，標本化周波数の変換が必要となる．
標本化周波数 fs を Ufs/D （U, D は整数）にするには，図 4.4 のように，Ｕ倍のレートでア
ップサンプリング(up-sampling)し，低域通過フィルタを通した信号を，1/D のレートでダウ
ンサンプリング(down-sampling)する.アップサンプリングでは，標本値系列 x (n) に，値がゼ
ロの標本点を(U-1)個追加し，標本化周波数を U 倍にする．すなわち，アップサンプリング後
の系列 v (n) は
n kU
x(n U )
v ( n)
otherwise
0
(4.15)
となる．次に，ダウンサンプリングでは，標本値系列から，Ｄ個毎の値を取り出し，他は捨て，
標本化周波数を 1/D 倍にする．すなわち，アップサンプリング後の系列 w(n ) は次式となる．
w( n) v ( Dn)
(4.16)
D : integer
信号処理 1 章-4 章.docxx
- 26 -
大西
v (n)
w(n )
ダウンサンプラ
低域通過フィルタ
fs
アップサンプラ
x(n)
昇（名古屋大学）
y (n )
Uf s / D
↑U
↓D
図 4.4
標本化周波数の変換
もし，標本化周波数をＵ倍にするだけなら，図 4.4 のダウンサンプラは不要である．逆に，
1/D 倍にするのなら，アップサンプリングは不要となる．なお，図中の低域通過フィルタは，
エイリアシングを回避するためのもので，次の特性を有する．
U
min( D , U )
H ( )
otherwise
0
(4.17)
ここで， は正規化角周波数である．マルチレートシステムの詳細は文献(9)や(10)を参照の
こと．
4.6 D/A 変換
計算機で信号処理した結果に基づいて，モータなどのアナログ機器を制御することがある．
その場合，ディジタル信号からアナログ信号への変換である D/A 変換が必要になる．その変換
を実現する装置を D/A コンバータと言い，それは，図 4.5 のように，復号，再生フィルタおよ
び後置フィルタからなる．復号は，2 進表現されたディジタル信号を周期 T 毎のアナログ振幅
値をもつ離散系列（図の a）)へ変換することで，再生フィルタは，アナログの瞬時値を周期 T
の間保持して，ステップ列（図の b）)として時間的に連続な信号とする．そして，後置フィ
ルタは滑らかな波形とするための，低域通過フィルタである．
後置
再生
復号
フィルタ
フィルタ
a)
c)
b)
T
T
図 4.5
T
D/A コンバータの構成
D/A での問題として，アパーチャ効果がある．これは，図 4.5 の a)のインパルス波形を，b)
信号処理 1 章-4 章.docxx
- 27 -
大西
昇（名古屋大学）
の有限時間の現象で近似することによる高域周波数での特性の劣化で，D/A 変換性能悪化の要
因である．信号の帯域に比べ，サンプリング周波数が十分高ければ（T が十分小さければ），
アパーチャ効果は小さい．D/A 変換の詳細は文献(11)を参照のこと．
信号処理 1 章-4 章.docxx
- 28 -
大西
昇（名古屋大学）
5. 離散時間信号とその変換
時間軸が離散化された離散時間信号の変換としての，離散時間フーリエ変換，離散フーリエ
変換および z 変換を説明する．まず，離散時間フーリエ変換のスペクトルは周期的な連続関数
であるが，離散フーリエ変換のスペクトルは周期的な離散系列となる．そして，z 変換はアナ
ログ信号のラプラス変換に相当し，離散時間信号の複素平面への変換である．
5.1 離散時間信号の表現
離散時間信号を x(n) （ｎは整数）と表わす．標本化周期 T を明示する必要がある場合は，
x(nT)とする．以下に，代表的な離散時間信号を示す．
例１ 単位インパルスδ(n)
1 n0
0 n 0
(n)
(5.1)
例２ 単位ステップ u(n)
1 n 0
u ( n)
0 n 0
(5.2)
δ(n)と u(n)とには，次式の関係が成立する．
( n) u ( n) u ( n 1)
例３ 指数関数信号
(5.3)
x(n) a n u (n)
x(n) e jn0T
例４ 複素正弦波信号
次式を満足する信号を，因果的な信号と言う．
x(n) 0
(5.4)
n0
5.2 離散時間フーリエ変換
離散時間信号 x(n)をフーリエ変換するために，デルタ関数δ(t)を使って，次式のように，
x(n) を連続時間信号 xｓ(t)（添え字の s は標本化の意味）として表現する（文献(12)）．
xs (t )
x(n) (t nT )
(5.5)
n
xs(t)のフーリエ変換 Xs(ω)は次式となる．
X s ( )
[ x(n) (t nT )]e jt dt
n
x(n)e
x(n)[
n
(t nT )e jt dt ]
(5.6)
jnT
n
上式で， T とおき，添え字の s をとると
X ()
x ( n )e
jn
(5.7)
n
信号処理 5 章-7 章.docx
- 29 -
大西
昇（名古屋大学）
となる．式(5.7)を，離散時間信号 x(n)のフーリエ変換（discrete-time Fourier transform）と言う．
上式は， X () のフーリエ級数展開で，変数Ωについて周期 2πの周期関数である． X ( ) の逆
変換は
x ( n)
1
2
X ()e jn d
(5.8)
で，x(n)は周期関数 X ( ) のフーリエ係数に相当する．
離散時間信号のフーリエ変換の例を示す．
a) 矩形信号（パルス）
1 0 n N 1
PN ( n)
0 その他
N 1
X ( )
e jn
n0
1 e jN
1 e j
sin( N )
2 e j ( N 1) / 2
sin( )
2
b) 片側指数関数信号
a n
x(n) a n u (n)
0
X ()
a n e jn
n0
n 0
n 0
(ae
j n
)
n0
1
1 ae j
ae j 1
上記の ae j 1 は，無限級数の和が収束するための条件である．
5.3 離散フーリエ変換
1) 直観的な説明
図 5.1(a)のように，有限時間継続信号 x(t)があり，そのフーリエ変換を X(ω)（ωmax で帯
域制限されている）とする．図(b)のように，x(t)を T<π/ωmax で標本化すると，標本化され
た信号 xs (t)のフーリエ変換 Xs(ω)は，周期がωs=2π/T の周期関数となる．さらに，Xs(ω)
を周波数間隔ω0=ωs /N=2π/(NT)で標本化すると，時間軸では周期 NT の周期離散時間信号
xs(nT)となる．その結果，図(c)のように，時間軸と周波数軸の双方で，離散的かつ周期的と
なる．最後に，図(d)のように，時間軸，周波数軸のそれぞれで，一周期分だけ取り出した x(n)
と X(k)が，離散フーリエ変換の対である．
信号処理 5 章-7 章.docx
- 30 -
大西
昇（名古屋大学）
X ( )
x(t )
max
max
X s ( )
x s (t )
max
max
X s ( k 0 )
x s (nT )
N 0
0
X (k )
x(n)
( N 1) 0
図 5.1
離散フーリエ変換導出の図的説明
(a)アナログ信号とそのスペクトル，(b)離散時間信号とそのスペクトル
(c)周期離散信号とそのスペクトル，(d)離散フーリエ変換の対
2) 離散フーリエ変換対の導出
離 散 時 間 信 号 の フ ー リ エ 変 換 で あ る 式 (5.7) を 用 い る と ， 図 5.1 (b) の
x s (t )
n
x ( n ) (t nT ) のフーリエ変換 X s () は
X s ()
N 1
n
n0
x(n)e jnT
x(n)e
jnT
(5.9)
となる（2 つ目の等号は， x(n)の系列長がＮだから）．図(c)のように，周波数軸を間隔
0 s N 2 (NT) で離散化すると，ω=kω0 での Xs(ω)は次式のように求まる．
X s (k 0 )
上式で， W N e
j
N 1
x ( n )e
jnk0T
n0
2
N
N 1
x ( n )e
jnk
n0
2
T
NT
N 1
x ( n )e
jnk
2
N
(5.10)
n0
とすると，次式となる．
N 1
X (k )
x(n) W
kn
N
0 k N 1
(5.11)
n 0
これが離散時間信号 x(n) (n=0,1,…,N-1)の離散フーリエ変換（discrete Fourier transform；DFT
と略記）である．そして，離散逆フーリエ変換は次式となる．
信号処理 5 章-7 章.docx
- 31 -
大西
x ( n)
1
N
N 1
X (k ) W
kn
N
0 n N 1
昇（名古屋大学）
(5.12)
k 0
(5.11)と(5.12)が，離散フーリエ変換の対となることを以下で確認する．(5.11)に(5.12)を
代入すると，
N 1
X (k )
x ( n) W Nkn
n 0
ここで，
1
N
N 1
W
nr
N
n 0
X (k )
1
N
N 1
1
N
n 0
i 0
1 r mN
0 その他
N 1
N 1
N 1
X (i)
i 0
W
1
X (i ) W Nin W Nkn
N
N 1
N 1
X (i )
i 0
W
( k i ) n
N
(5.13)
n 0
（m は整数）であるので
( k i ) n
N
n 0
1
X (k ) N X (k )
N
(5.14)
したがって，変換対が成立する．
3) 離散フーリエ変換の性質（詳細は文献(13)を参照）
f ( n ) F ( k )，
a) 対称定理
g ( n ) G ( k ) とする．
有限長の系列 f(n)(n=0,1,..,N-1)が実数ならば
Re[ F ( k )] Re[ F ( N k )]
F (k ) F ( N k )
b) 線形性
Im[ F ( k )] Im[ F ( N k )]
(5.15)
arg[ F ( k )] arg[ F ( N k )]
(5.16)
(5.17)
af ( n ) bg ( n ) aF ( k ) bG ( k )
もし，f(n), g(n)の系列長が異なり，それぞれ N f
, N g であるならば，式(5.17)の変換対
の系列長は max(N f , N g ) である．
N 1
c) [
~
f (m)g~(n m)]P (n) F (k )G(k )
N
m0
~
f (n), g~(n) は，
（循環畳み込み定理）
(5.18)
f ( n ), g ( n ) を周期 N で拡張した周期系列である．図 5.2 に，循環畳み込
みの例を示す．
d) 循環推移定理
ここで， f ( n m )
f (n m )
N
N
(5.19)
W N km F ( k )
とは，周期 N の系列
~
f (n) を m だけシフトしたのち，基本周期分（0 か
ら N-1 まで）取り出した信号である．すなわち， f ( n m )
N
~
f ( n m ) PN ( n )
．ここで，PN (n) は
信号処理 5 章-7 章.docx
- 32 -
大西
昇（名古屋大学）
長さ N の矩形信号である．
N 1
e) f (n) g (n)
~
1
~
[ F (m)G (k m)]PN (k )
N m 0
(5.20)
~
~
右辺は，周期系列 F (n), G(n) を畳み込んだ後，その一周期分を取り出したものである．
図 5.2
離散畳み込みと循環畳み込み(x(n)=y(n)=P5(n)の場合)
(a)は，離散畳み込み w(n)
x(k ) y (n k )
k
N 1
(b)～(e)は， u (n) [
~x (k ) ~y (n k )]P
．
N ( n) で，周期系列の畳み込み（循環畳み込み）
k 0
4) 離散フーリエ変換とフーリエスペクトル
連続時間信号の f (t ) とそのフーリエ変換を F () とする． f ( t ), F ( ) をそれぞれ周期 Tp,ωs
（ 2 / T ）で拡張し，周期関数 ~f ( t ),
~
f (t )
f (t nT ),
~
F ( )
p
n
~
f (t )
をT
Tp / N
~
F ( )
をつくる．すなわち，
F ( n )
s
n
で， F~ ( ) を
0 s / N で，それぞれ標本化する．得られた離散系列
~
~
f (nT ), F (k0 ) の間には次式が成立する．
~
1
T f (nT )
N
N 1
~
F (k )W
k 0
0
kn
N
(5.21)
信号処理 5 章-7 章.docx
- 33 -
大西
N 1
~
F (k0 )
~
昇（名古屋大学）
T f (nT )W
kn
N
n 0
(5.22)
このように，アナログ信号 f(t)の標本値 T ~f ( nT ) と，フーリエスペクトル F(ω)の標本値
~
F (k 0 ) とが，離散フーリエ変換の対になる．
5.4 z 変換
1) z 変換の定義
離散時間信号 x(n)のｚ変換（z-transform）は次式で定義される．
X ( z)
両側 z 変換
x ( n) z
n
(5.23)
n
X ( z)
片側 z 変換
x(n)z
n
(5.24)
n 0
ここで，ｚは複素変数である．ｚ変換はラプラス変換と同様に，収束領域でのみ定義される
ので，収束領域の明記が必要である．
x(n)が因果的(x(n)=0 n<0)なら，両側 z 変換＝片側 z 変換となる．因果的な信号の収束領
域は，半径 R（R の値は信号により異なる）の円の外側で， z を含む．すなわち z R と
なる．
X(z)の極と零点を説明する．極(pole)は X(z)の値を発散させる変数 z の値（特異点）で，
X(z)の値をゼロとする変数 z の値である．極は，逆ｚ変換で重要となる．
零点(zero)は
ｚ変換の例を示す．
a) 矩形信号（パルス）
1 0 n N 1
PN ( n)
0 その他
X ( z)
N 1
n
n 0
PN (n) z n z n
1 z N
1 z 1
収束領域は z=0 を除く全平面．
b) 片側指数関数信号
x(n) a n u (n)
X ( z)
a nu(n) z n
n
(az
n0
1 n
)
1
1 az 1
収束領域は az 1 1
(5.25)
c) 両側指数関数信号
信号処理 5 章-7 章.docx
- 34 -
大西
昇（名古屋大学）
x ( n ) a n u ( n ) b n u ( n 1)
X ( z)
収束領域は a z b
1
1
1
1 bz
1 az 1
(5.26)
2) z 変換の性質
f ( n ) F ( z )， g ( n ) G ( z ) とする．
a) f * ( n ) F * ( z *)
F ( z)
証明
(5.27)
f (n) z n
f (n) z *
F ( z*)
だから
n
n
．
n
両辺の共役をとると
F * ( z*)
f * ( n) z
n
n
b) 線形性
(5.28)
af ( n ) bg ( n ) aF ( z ) bG ( z )
c) 時間軸推移
(5.29)
f (n k ) z k F ( z )
d) f ( n ) F (1 / z )
(5.30)
e)
f (k ) g (n k ) F ( z)G( z)
（離散畳み込み演算）
(5.31)
k
2j
f) f ( n ) g ( n ) 1
C1
z
dv
1
F ( )G ( v )
v
v
2j
C2
z dv
F ( v )G ( )
v v
（複素畳み込み積分）
(5.32)
積分路 C 1 , C 2 はそれぞれ，F ( z v ) と G (v ) ，F (v ) と G ( z v ) の収束領域の共通部分にとる．
証明（5.5 の逆ｚ変換を利用）
f (n) g (n) z
n
1
2j
C2
n
2j
1
F (v )v n 1 dv g ( n ) z n
Cf
n
1
2j
Cf
F (v )
g (n)( v )
z
n
n
dv
v
(5.33)
z dv
F ( v )G ( )
v v
g) パーセバルの定理
1
1
f (n) g * (n) 2j F (v)G * ( v *)
n
C
dv
v
(5.34)
積分路 C は F(z)と G(z)の収束領域の共通部分にとる．
証明（5.5 の逆ｚ変換を利用）
左辺＝
1
n 2j
n 1
F ( z) z dzg * (n)
1
2j
n
F ( z) g * (n)z
n
dz
1
z
2j
1
F (v)G * ( v *)
C
dv
v
Q.E.D.
信号処理 5 章-7 章.docx
- 35 -
大西
昇（名古屋大学）
単位円を積分路にとることができ，x(n)=y(n)なら，式(5.34)は次式となる．
f (n)
2
n
1
2
2
F (e j ) d
(5.35)
5.5 逆 z 変換とその求め方
X(z)の逆ｚ変換は次式である．
x(n)
1
2j
X ( z) z
n1
(5.36)
dz
C
ここで，積分路 C は原点周りを反時計方向にまわる閉積分路で，X(z)の収束領域内にとる．
逆ｚ変換の求め方は，べき級数展開，部分分数展開および留数計算の３通りがある．以下で
は，収束領域が， z R （R は任意）である場合，すなわち逆変換 x(n)が因果的な信号である
場合について説明する．なお，収束領域が， z R である場合，留数計算の方法はそのまま使
用できるが，べき級数展開では除算を工夫し，部分分数展開では，非因果的な信号の z 変換対
の表を新たに用意する必要がある．
1) べき級数展開
X(z)が次式のように，べき級数の形で与えられている場合，X(z)の逆変換 x(n)は， z n の
係数である．
X ( z)
x ( n) z
n
n
例１
X (z) a0 a1 z 1 a2 z 2
べき級数の形で与えられているので， x(0)
a0 , x(1) a1, x(2) a2 で，他の x(n)=0．
もし，X(z)が解析的に閉じた表現である場合，べき級数展開を用いる．
例２
X (z) log(1 az1 )
収束領域は z a
log( 1 x ) のべき級数展開を用いると
X ( z)
(1) n 1 a n z n
n
n 1
となるので，
信号処理 5 章-7 章.docx
- 36 -
大西
昇（名古屋大学）
( 1) n 1 a n / n
n 1
x(n)
n 1
0
また，X(z)が有理多項式で与えられる場合，除算を繰り返して，べき級数展開を得る．
例３
X ( z)
収束領域は z a
1
1 az 1
収束領域が z a より，x(n)は右側数列（n<n1 で x(n)=0 となる数列）であるので， z 1 の
べき級数となるように除算を繰り返し，
1
1 az 1 a 2 z 2
1
1 az
x(n) anu(n)
だから，
もし，収束領域が z a なら x(n)は左側数列（n>n2 で x(n)=0 となる数列）であるので，
z のべき級数となるように除算を繰り返し，
z
a 1 z a 2 z 2
az
x(n) anu(n 1)
だから，
2) 部分分数展開（ｚ変換対の表 5.1 を利用）
X(z)が有理多項式，すなわち
X ( z)
q0 q1 z 1 qm z m
1 p1 z 1 pn z n
（m≦n）
(5.37)
である場合，ラプラス逆変換と同様に，部分分数展開で逆変換は得られる．分母の多項式を
複素根まで許して因数分解する．
(a)すべて単根 αi(i=1,2,...,n)の時
X(z)は次式のように部分分数に展開されるとする（ i (i 1,, n) は全て異なる）
．
n
X ( z) a0
i 1
ai
(5.38)
1 i z 1
上式のパラメータは次式で求まる．
a0 X (0), ai (1 i z 1 ) X ( z) z i 1,2, , n
変換対
1 u (n), a n u ( n)
(5.39)
i
1
1 az 1
を用いると，逆変換は次式となる．
信号処理 5 章-7 章.docx
- 37 -
大西
昇（名古屋大学）
n
x(n) (a0
a
i
n
i )u(n)
(5.40)
i 1
(b) k 重根 β をもち，他は全て単根 αi(i=1,2,...,n-k)の時
式(5.37)の分母の多項式が (1 z
1 k
)
を因数に持ち，他は全て異なる単根とすると，X(z)
は次式のように部分分数に展開される．
X ( z ) a0
nk
ai
1 z
i 1
1
k
bi
(1 z
i 1
i
(5.40)
1 i
)
上式のパラメータは次式で求まる．
a0 X (0),
ai (1 i z 1 ) X ( z ) z i i 1,2, , n k
1
d (k i )
1 k
bi
{(
1
z
)
X
(
z
)}
( k i )
k i
1
(k i)! ( ) dz
z
変換対 1 u ( n ),
a n u (n)
i 1,2, , k
(5.41)
1
,
1 az 1
(n 1)(n 2) (n l 1) n
1
a u(n)
を用いると，式(5.40)の逆変換として，
(l 1)!
(1 az 1 ) l
次式を得る．
nk
k
i 1
i 1
x(n) (a0 ai in bi
(n 1) (n (i 1)) n
)u(n)
(i 1)!
(5.42)
3) 留数計算
X(z)の逆変換は次式のように留数定理を用いて計算できる．
x(n)
[ X (z) z
1
X ( z) z n1 dz
2j
C
n1
の極 z i における留数]
(5.43)
zi
上式で，留数計算の対象となる極 zi は，X(z)の収束領域内に取った積分路 C の内側（原点
を含む）に存在する X ( z) z
n 1
の極である．式(5.43)の極 zi における留数 Res ( z i ) は次のよう
に求める．
a) 極 zi が単極の場合， Res ( z i ) ( z z i ) X ( z ) z n 1
(5.44)
z zi
d m1
1
( z z i ) m X ( z ) z n1
b) 極 zi が m 重極の場合， Res( z i )
m 1
(m 1)! dz
z zi
(5.45)
信号処理 5 章-7 章.docx
- 38 -
大西
昇（名古屋大学）
4) 逆ｚ変換の計算例
X ( z)
z
z3
（収束領域は|z| > 3）の逆ｚ変換を，部分分数展開と留数計算の 2 通りで求
める．
a) 部分分数展開（ｚ変換対の表を利用）
z
1
z 3 1 3 z 1
X (z)
表 5.1 の
a n u (n )
←→
1
より，
1 az 1
x ( n ) ( 3) n u ( n )
b) 留数計算
積分路 C は収束領域 z 3 にとる．そして着目する X ( z) z
n 1
の極は，C の内側（原点
を含む）領域に存在するものである．
X ( z ) z n 1
z
zn
z n 1
z3
z3
の極は
n≧0 の時，z = -3 のみ． x(n) ( z 3) X ( z ) z n 1
n=-1 の時，z = -3, 0．
x ( n) ( z 3) X ( z ) z 2
z 3
z 3
(3) n
zX ( z ) z 2
z 0
( 3) 1 3 1 0
n=-2 の時，z = -3 が単極，z=0 が２重極
x ( n ) ( z 3) X ( z ) z 3
z 3
d z 2 X ( z ) z 3
dz
( 3) 2 3 2 0
z0
n =-m の時，z = -3 が単極，z=0 が m 重極
x ( n) ( z 3) X ( z ) z ( m 1)
( 3) m
故に，
z 3
1
d m 1 z m X ( z ) z ( m 1)
( m 1)!
dz m 1
z0
m 1
( 1) ( m 1)! m
3 ( 3) m ( 3) m 0
( m 1)!
x ( n ) ( 3) n u ( n )
信号処理 5 章-7 章.docx
- 39 -
大西
表 5.1
昇（名古屋大学）
ｚ変換対
x(n)
X(z)
収束領域
(n )
1
全平面
u (n )
a n u (n )
na n u (n )
n 2 a n u (n)
sin( n )u ( n)
cos( n ) u ( n )
e n sin( n ) u ( n )
e n cos( n ) u ( n )
( n 1) a n u ( n )
( n 1)( n 2 ) ( n m 1) a n u ( n )
( m 1)!
1
z 1
1 z 1
1
z a
1 az 1
az 1
z a
(1 az 1 ) 2
a(1 az 1 ) z 1
z a
(1 az 1 ) 3
sin( ) z 1
z 1
1 2 cos( ) z 1 z 2
1 cos( ) z 1
z 1
1 2 cos( ) z 1 z 2
e sin() z 1
1 2e
1
cos() z e
2 2
z
1 e cos() z 1
1 2e
1
cos() z e
1
(1 az 1 ) 2
1
(1 az 1 ) m
2 2
z
z e
z e
z a
z a
信号処理 5 章-7 章.docx
- 40 -
大西
昇（名古屋大学）
6. 離散時間システム
入力と出力が離散時間信号であるシステムについて説明する．システムとして望ましい性質と
して，線形性，時不変性，因果性，安定性がある．線形時不変システムの特性はインパルス応答
の z 変換であるシステム関数で記述でき，周波数特性はシステム関数から求まる．システムの構
成方法，複数システムの接続も説明する．
6.1 システムとは
システムとは，特定の働きをするように組み立てられ
x(n)
た装置の組合わせであるが，ここでは，入力信号を変換
システム
y(n)
S
もしくは処理した結果を出力するものとする．そして，
離散時間システムでは，入力x(n)と出力y(n)がともに離
散時間信号である．離散時間システムは
図 6.1 離散時間システム
y(n)=S[x(n)]
と表現される（図 6.1 参照）．
入力がインパルスδ(n)である時の，システムの出力をインパルス応答(impulse response) と
言い，h(n)で表す．
6.2 システムとしての重要な性質
システムは，線形システムと非線形システムに大別される．線形システム(linear system)では，
次式の線形性が成立する．
y 1 ( n ) S [ x 1 ( n )] , y 2 ( n ) S [ x 2 ( n )] とする時，
S [ k 1 x 1 ( n ) k 2 x 2 ( n )] S [ k 1 x 1 ( n )] S [ k 2 x 2 ( n )] k 1 y 1 ( n ) k 2 y 2 ( n )
(6.1)
ここで， k 1 , k 2 は任意の定数である．
例１
y ( n) ax (n ) bx (n 1)
は線形システムで， y ( n )
ax ( n ) bx 2 ( n )
は非線形システムであ
る．
線形の証明
y 1 ( n ) S [ x 1 ( n )] , y 2 ( n ) S [ x 2 ( n )] とする．
k1 x1 ( n ) k 2 x 2 ( n ) の出力 y(n) は，
y (n) S[k1 x1 (n) k 2 x 2 (n)] ak1 x1 (n) k 2 x 2 (n) bk1 x1 (n 1) k 2 x 2 (n 1)
k1 ax1 (n) bx1 (n 1) k 2 ax 2 (n) bx 2 (n 1) k1 y1 (n) k 2 y 2 (n)
非線形の証明
y 1 ( n ) S [ x 1 ( n )] とする．
2 ak1x(n) bk12 x2 (n)
入力 k 1 x1 ( n ) の出力 y(n) は， y(n) ak1 x(n) b k1 x(n)
一方， k1 y1 (n)
で，式(6.1)の関係は成立せず．故に，非線形．
k1 ax1 (n) bx12 (n)
システムの入出力特性が時間とともに変化しない，すなわち次式が成立するシステムを時不
変システム(time-invariant system)と言う．
信号処理 5 章-7 章.docx
- 41 -
大西
y ( n ) S [ x ( n )] の時，
例２
昇（名古屋大学）
y ( n k ) S [ x ( n k )]
(6.2)
y ( n) ax (n ) bx (n 1) は時不変システムで， y ( n ) nx ( n ) は時変システムである．
時不変の証明
S [ x ( n k )] は，与式の n を n-k で置き換えて，
S [ x ( n k )] ax ( n k ) bx ( n k 1)
一方， y ( n k ) は， y (n) ax (n) bx (n 1) を k だけ遅らせたものであるので，
y ( n k ) ax ( n k ) bx ( n k 1) となる．故に，式(6.2)の関係が成立する．
時変の証明
S [ x ( n k )] ( n k ) x ( n k )
一方， y ( n k ) は y ( n ) nx ( n ) を k だけ遅らせたものであるので， y ( n k ) nx ( n k )
故に，式(6.2)の関係は成立せず．
線形性と時不変性が成立するシステムを，線形時不変システム(linear time-invariant system,
LTI)と言う．システムの出力 y(m)が，n≦m である x(n)のみで決まり，未来の入力に依存しな
いシステムを因果的なシステム(causal system)という．
最後に，有界入力に対する出力が有界となるシステムを，安定なシステム(stable system)と言う．
6.3 線形時不変システム
線形時不変システムの出力は，次式のように，入力 x(n)とインパルス応答 h ( n ) S [ ( n )] の
畳み込みで表現できる．
y ( n) x ( n) * h( n)
k
k
x ( k ) h ( n k ) h( k ) x ( n k )
(6.3)
その証明はつぎの通りである． x(n)
x(k ) (n k ) であるので，
k
k
k
k
y ( n) S [
x(k ) (n k )] x(k )S[ (n k )] x(k )h(n k )
上式で，時不変性 h ( n k ) S [ ( n k )] を使用している．Q.E.D.1
インパルス応答 h(n)が因果的なら線形時不変システムは因果的となる（その逆も真）．n<0
について h ( n ) 0 であるので，出力 y(n)は
y ( n)
n
k
k 0
x ( k ) h( n k ) h ( k ) x ( n k )
(6.4)
となり，出力 y(n)は k≦n である x(k)のみで決まる．さらに，入力 x(k)が因果的(x(n)=0
1
n<0)
Q.E.D. quod erat demonstrandum（ラテン語）の略で，証明終わりの意味で使う．
信号処理 5 章-7 章.docx
- 42 -
大西
なら， y ( n )
n
k 0
昇（名古屋大学）
x ( k ) h ( n k ) となる．
線形時不変システムが安定であるための必要十分条件は
h( n)
(6.5)
n
である．証明は以下の通り．
入力は有界であるので， n x(n) X u とする．出力の大きさは次式となる．
h( k ) x ( n k )
y ( n)
k
h(k )
したがって，
h( k ) x ( n k )
n
h( k ) X
n
u
な ら y(n) で あ る （ 十 分 条 件 ）． 必 要 条 件 の 証 明 は ，
n
h(n) なら，出力は有界とならないこと（対偶）を示す．
n
h(n)
となる h(n)
n
を用いて，つぎの有界な入力を構成できる．
h * ( n) / h( n)
h(n) 0
x(n)
h(n) 0
0
ここで，記号*は複素共役である．この有界入力に対する，n=0 での出力は
k
n
h(k)x(k )
y(0)
h(k )
2
h(k )
h(k)
n
となり，有界でない．
6.4 差分方程式と離散時間システム
システムの特性は，連続時間の場合は微分方程式で表わされるが，離散時間の場合は差分
方程式で表わされる．そして，線形時不変な離散時間システムは次式の定係数差分方程式で
記述できる．
y ( n)
M
a x( n i ) b y ( n
i0
例
N
i
j 1
１次の差分方程式
j
(6.6)
j)
y ( n) x(n) by (n 1)
単 位 イ ン パ ル ス x(n) (n) に 対 す る 応 答 は ， 因 果 的 な
y ( n ) b n u ( n 1)
y (n) b n u (n)
と，非因果的な
の 2 通りある．
求め方１
信号処理 5 章-7 章.docx
- 43 -
大西
昇（名古屋大学）
境界条件 y (n) 0 (n 0) として， y ( 0 ) x ( 0 ) 1, y (1) by ( 0 ) b , と順次求めると，
y (n) b n u (n) ．
境界条件 y (n) 0 (n 0) として， y ( n 1) y ( n ) x ( n ) / b より，
y(1) {y(0) x(0)}/ b 1 / b, y(2) {y(1) x(1)}/ b b2 , と順次求めると，
y ( n ) b n u ( n 1)
．
求め方２
Y ( z)
1
収束領域が
1 bz 1
Y (z)
b 1 z
1 b 1 z
z b
収束領域が
として留数計算して求めると， y ( n ) b n u ( n ) ．
z b
として留数計算して求めると， y ( n )
b n u ( n 1)
．
この例から，定係数差分方程式が規定するシステムは，因果的でないシステムをも含むこ
とが分かる．したがって，定係数差分方程式だけでは，出力は一意に決まらないので，シス
テムが因果的か否かを明記する必要がある．
6.5 システムの分類
インパルス応答の持続時間（系列の長さ）に着目した分類と，出力のフィードバックの有無
に基づくシステムの分類を説明する．
1) FIR と IIR
インパルス応答 h ( n ) S [ ( n )] の持続時間（系列長）が有限か無限かで，システムは分類され
る．
インパルス応答の持続時間が有限であるシステムを FIR(finite impulse response)システム，
無限であるのを IIR(infinite impulse response)システムという．
FIR となるインパルス応答の例は矩形信号（パルス）で，IIR の例は片側指数関数信号であ
る．
2) 再帰型と非再帰型システム
再帰型システムは，過去の出力値を計算に利用し（フィードバックがあり）
，FIR システム，
IIR システムのどちらにもなり得る．
非再帰型システムは，出力の表現に過去の出力を用いない（フィードバックがない）．すな
わち，式(6.6)で，すべての bj がゼロである．インパルス応答は有限長（FIR システム）であ
る．
6.6 システム関数
システムの特性を表すシステム関数 H ( z ) は，次式で定義される．
H ( z) Y ( z)
(6.7)
X ( z)
システム関数はインパルス応答 h(n)のｚ変換である．すなわち，
信号処理 5 章-7 章.docx
- 44 -
大西
h( n) z
H ( z)
昇（名古屋大学）
n
(6.8)
n
y(n)
その証明は次の通りである．
x(k )h(n k ) の 両 辺 を
z 変換すると，
k
Y ( z ) X ( z ) H ( z ) ．よって，
ある．
Y ( z)
X ( z) は H ( z ) ，すなわちインパルス応答 h(n)のｚ変換で
Q.E.D.
差分方程式とシステム関数の関係を示す．ｚ変換の推移定理を利用して，差分方程式(6.6)
の両辺をｚ変換すると，
Y ( z)
M
a
m0
[1
m
X ( z)
mz
N
b z
n
n
n 1
N
b z
n
n
n 1
Y ( z)
M
]Y ( z ) [ a m z
(6.9)
m
]X ( z)
m0
Y(z)/X(z)を求めるために整理すると，次式を得る．
M
H ( z) [
a
mz
m
N
] [1
m 0
b z
n
n
(6.10)
]
n 1
システム関数 H(z)の極は H(z)の値を発散させる変数 z の値（特異点）で，零点(zero)は H(z)
の値を零とする変数ｚの値である．N > M として，H(z)が次式のように因数分解できる時，
M
H ( z) C (1 m z 1 )
m 1
零点は z
N
(1
n
z 1 )
(6.11)
n 1
m (m 1,2, , M)とz 0 ．極は z n (n 1,2, , N) となる．
線形時不変システムの安定性の必要十分条件は，式(6.5)である．この条件は，「 H(z)の極
βi が全て単位円内(境界は除く)に存在し（ i 1），かつ H(z)の収束領域が単位円を含むこ
と（システムの因果性を仮定すれば後半の条件は不要）
」と等価である．
簡単な証明を示す． H(z)の全ての極は単極のみとすると，システム関数は次式のように，
部分分数に展開される．
H ( z)
P
ki
1 z
i 1
1
i
信号処理 5 章-7 章.docx
- 45 -
大西
この逆 z 変換は， h(n)
P
k
i
i 1
n
i
となる．全ての極が
昇（名古屋大学）
i 1なら，次式のような展開が得ら
れる．
h( n)
よって，
P
P
P
i 1
i 1
i 1
k i in K in K i
h( n)
n0
n0
P
K i
n
i 1
但し k i K
n
P
K i
n
Q.E.D.
i 1 n 0
6.7 離散時間システムの周波数特性
入力が正弦波 e jn T の時のシステムの出力，すなわちシステムの正弦波応答は次式で与えら
れる．
k
k
h(k )e j (n k )T [ h(k )e jkT ]e jnT H (T )e jnT
y(n) h(n) * e jnT
(6.12)
最後の等号で，
H (T ) H ( z) z e jT
h(k )e
jkT
(6.13)
k
を使用している．式(6.12)は，出力 y(n)は H (T ) と入力 e jn T の積になることを意味する．
すなわち， H (T ) がシステムの周波数特性に相当する．
式(6.13)において， T （ω(rad/s)を標本化周波数 1/T で割っているので，Ωを正規化
角周波数(rad)と言う）とすると，次式となる．
H ()
h(k )e
jk
(6.14)
k
ここで， H () は周期 2πの周期関数であることに注意．
以上をまとめると，離散時間システムの周波数特性は次式で与えられる．
H () H ( z) z e j
h(k )e
jk
(6.15)
k
H ( ) を次のように極座標表現
H () H () e j ()
(6.16)
した時， H () を振幅特性(amplitude characteristic)， ()（実数）を位相特性(phase characteristic)
という．さらに，群遅延(group delay)特性は
( ) d ( ) d
(6.17)
信号処理 5 章-7 章.docx
- 46 -
大西
で表わされる．なお， H()
昇（名古屋大学）
A()e j () とし，負値も取
H ()
りうる実数 A ( ) を振幅特性と定義した本もある（文献
5
4
(14)）．
3
2
以下に周波数特性の例を示す．
1
例1
FIR システム
H ( z ) 1 z 1 z N 1
-3
N 1
H ()
e jk
k 0
-2
1
2
3
Ω
sin(N 2) j ( N 1) / 2
e
sin( 2)
sin(N 2)
A()
,
sin( 2)
-1
(a)振幅特性
() 6
4
() ( N 1) / 2
Ω
2
-3
-2
-1
1
2
3
-2
-4
よって，
-6
-8
sin(N 2)
H ()
,
sin( 2)
(b)位相特性
図 6.2
A() 0
( N 1) / 2
()
A() 0
( N 1) / 2
FIR システムの周波数
特性
図 6.2 に N=5 の場合の， H () と ( ) を示す． A() が負の場合， 1 e j であるので，
位相特性は-πだけ変化させている．
例2
１次 IIR システム
H(z) (1 b) (1 bz1 ), b 1
H ( ) (1 b )
(1 b 2 2b cos ) ,
( ) arctan[ b sin (1 b cos )]
b=0.8 の時の， H () ， ( ) を図 6.3 に示す．
例3
2 次 IIR システム
H ( z)
1
1 2r cos z 1 r 2 z 2
r=0.8, θ=π/4 の時の， H () ， ( ) を図 6.4 に示す． H () で，Ω=π/4 に共振の
ピークが見られる．
信号処理 5 章-7 章.docx
- 47 -
大西
H ()
昇（名古屋大学）
1
4
0.8
3.5
H ()
0.6
3
2.5
0.4
2
1.5
0.2
1
-3
-2
-1
1
0.5
Ω
2
3
(a)振幅特性
-3
-2
0.5
3
()
Ω
0.25
-1
2
(a)振幅特性
0.75
-2
1
Ω
()
-3
-1
1
2
1
3
-0.25
0.5
-0.5
Ω
-0.75
-3
-2
-1
1
2
3
-0.5
(b)位相特性
図 6.3
-1
IIR システムの周波数
(b)位相特性
図 6.4
特性(b=0.8)
IIR システムの周波数
特性(r=0.8, θ=π/4)
6.8 システムの構成
1) 構成要素
離散時間システムを構成する要素は，
つぎの 4 つである（図 6.5 参照）
．
加算器
y ( n ) x1 ( n ) x 2 ( n )
減算器
y ( n ) x1 ( n ) x 2 ( n )
乗算器
y ( n ) ax ( n )
単位時間遅延素子
y ( n ) x ( n 1) ．
本書では，単位時間遅延素子を z-1
図 6.5
離散時間システムの構成要素
で表現するが，Delay の頭文字の D
で表すこともある．
2) 非再帰型システムの構成
次式であらされる非再帰型システムの構成法
には，直接型と転置形の 2 通りがある．
N 1
y(n)
h(k )x(n k )
k 0
直接形構成は，システム関数を式どおりに構成
したもので，図 6.6(a)となる．
転置形構成は，システム関数の次式のような展
開に基づいて構成したもので，直接形構成の，
図 6.6
- 48 -
非再帰型システムの構成
信号処理 5 章-7 章.docx
(a)直接形，(b)転置形
大西
昇（名古屋大学）
入出力端子と信号の向きを反転させたものになっている(図 6.6(b))．
H ( z ) h ( 0 ) z 1 [ h (1) z 1 [ h ( 2 ) z 1 [ h ( N 2 ) z 1 h ( N 1)] ]
3) 再帰型システムの構成
再帰型システムのシステム関数は次式で与えられる．
M
H ( z) [
a m z m ] [1
m 0
N
b z
n
n
]
n 1
直接形構成
M
P( z)
N
b z
a m z m , Q( z)
m0
n
n
とすると，
n 1
(6.18)
H ( z ) P ( z )[1 {1 Q ( z )}]
となる．直接形構成では， P ( z ) と 1 {1 Q ( z )} を縦続に接続する（図 6.7 参照）．
標準形構成
式(6.18)で，P ( z ) と 1 {1 Q ( z )} の順序を換えて，H ( z ) [1 {1 Q ( z )}] P ( z ) とし
ても良い．すると，図 6.8(a)となる．図(a)を見ると，遅延素子を共有化できることに気
づく．共有化した図(b)が標準形構成で，遅延素子数はシステム関数の分母・分子の次数
の最大値に等しい．
図 6.8
図 6.7
再帰型システムの直接形構成
再帰型システムの標準形構成
(a)接続順序を変えた構成，(b)標準形
信号処理 5 章-7 章.docx
- 49 -
大西
昇（名古屋大学）
6.9 システムの接続
H(z)
複数のシステムの接続の形態としては，縦続，並列，フ
Ha(z)
ィードバックの 3 通りがある．説明のために，2 つのシス
Hb(z)
テム Ha(z)と Hb(z)の，前記 3 つの接続を考える（図 6.9 参
(a)縦続
照）
．
縦続接続したシステム H(z)は，
H ( z) H a ( z)H b ( z)
(6.18)
＋
Hb(z)
並列接続したシステム H(z)は
H ( z) H a ( z) H b ( z)
(6.19)
(b)並列
フィードバック接続したシステム H(z)は
H a ( z)
H ( z)
1 H a (z)H b (z)
(6.20)
式(6.20)の導出を以下に示す．H (z)の入力，出力を X, Y
とし，Ha(z)の直前の入力を W とする．
Y H aW
X
＋
－
W
Ha(z)
Y
H(z)
Hb(z)
となる．
W X H bY
H(z)
Ha(z)
(c)フィードバック
図 6.9 システムの接続
(a)
(b)
(a) の両辺に Ha を掛け，(b)の関係を用いると，
HaW Ha X Ha HbY (1 Ha Hb )Y Ha X
よって，式(6.20)が得られる．
信号処理 5 章-7 章.docx
- 50 -
大西
昇（名古屋大学）
7. フィルタ
フィルタは特定の周波数成分のみを通すもので，信号処理の基本である．フィルタの特性記
述に用いる語句，フィルタの種類，ディジタルフィルタの設計法を説明する．設計法の一つに，
アナログフィルタに変数変換を施す間接設計法がある．そのために，代表的なアナログフィル
タと，アナログ低域通過フィルタから他の通過フィルタへの変換も説明する．
7.1 フィルタとその種類
フィルタリングとは，ある周波数帯域の信号成分を通し，他の成分の通過を阻止する操作で，
その装置（システム）をフィルタ(filter)という．以下では，アナログフィルタの周波数特性を
H ( ) で，ディジタルフィルタのそれを H () で表す．
理想的な低域通過フィルタでは，その周波数特性が
1
H IL ()
0
c
c
c : 遮断角周波数
(7.1)
となるが，実際には図 7.1 のようになる．
|H(Ω)|
δp
1
0.8
0.6
Ωc
Ωp
0.4
Ωs
0.2
δs
Ω
0.5
1
通過域
図 7.1
1.5
過渡域
2
2.5
3
阻止域
低域通過フィルタの振幅特性
フィルタの特性に関係する語句を以下で説明する（図 7.1 参照）
．
通過域(pass band) 信号を通す周波数帯域．より厳密には振幅特性
H ( )
が1 p
H ( ) 1
p
となる帯域[0,Ωp]
阻止域(stop band) 信号を通さない周波数帯域．より厳密には振幅特性
H ( )
が0
H ( ) s
と
なる帯域[Ωｓ,π]．
過渡域（遷移域）(transition band) 通過域と阻止域の間の周波数帯域[Ωp,Ωｓ]で，振幅特性は
なだらかに変化する．
遮 断 周 波 数 (cutoff frequency) Ω c
振 幅 特 性 の 最 大 値 H max ( ) で 正 規 化 さ れ た 振 幅 特 性
信号処理 5 章-7 章.docx
- 51 -
大西
H ( )
が 1
H max ()
昇（名古屋大学）
(A=3 dB)となる周波数．
2
減衰量(attenuation) A()
A() 20 log10
振幅の相対的な減衰量は次式で求まる．
H ()
H max ()
(7.2)
[dB]
フィルタは，周波数軸上での通過域のあり方により，図 7.2 のように，低域通過(lowpass)
フィルタ（図(a)）
，高域通過(highpass)フィルタ（図(b)）
，帯域通過(bandpass)フィルタ（図(c)），
帯域阻止(bandstop)フィルタ（図(d)）に分類される．特に，阻止域の周波数帯域が非常に狭い
帯域阻止フィルタをノッチフィルタと言う．なお，ディジタルフィルタの振幅特性 H () は
隅関数で， ( , ] を一周期とする，周期関数である．
｜H(ω)|
｜H(ω)|
ω
ω
ωLC
０
｜H(ω)|
｜H(ω)|
０
(a) 低域通過
ωHC
０
(b) 高域通過
図 7.2
ω
ω
ωLC
ωHC
０
(c) 帯域通過
ωLC
ωHC
(d) 帯域阻止
各種アナログフィルタの振幅特性
7.2 周波数特性の性質
アナログフィルタの周波数特性を H(ω) ( , ) ，ディジタルフィルタの周波数特性を
H(Ω) [0, 2 ) or ( , ] とする．インパルス応答（h(t)あるいは h(n)）が実数であるの
で，
H ( ) H * ( ) ， H ( ) H * ( )
(7.3)
である．したがって，
H ( ) H ( )
，
H ( ) H ( )
．
(7.4)
arg H ( ) arg H ( ) ， arg H ( ) arg H ( ) .
(7.5)
すなわち，振幅特性は偶対称，位相特性は奇対称である．
フィルタの重要な性質の一つに直線位相(linear phase)がある．直線位相とは，位相特性 ( )
が，Ωの一次式となることである．すなわち，群遅延が定数で，どの周波数成分も同じ遅延
時間である．直線位相の場合，図 7.3(b)のように，2 つの正弦波とも同じ遅延時間であるの
で．出力に波形ひずみを生じない．一方，(c)では遅延時間が異なるため，ひずみを生ずる．
表 7.1 に，アナログフィルタとディジタルフィルタの比較をまとめた．
信号処理 5 章-7 章.docx
- 52 -
大西
(a)システムの入力波形
(b)直線位相システムの出力
(c)非直線位相システムの出力
sin(t 10) sin{1.2(t 10) 0.2}
sin(t ) sin(1.2t 0.2)
図 7.3
表 7.1
昇（名古屋大学）
sin(t 10) sin{1.2(t 20) 0.2}
直線位相と非直線位相システムの出力
アナログフィルタとディジタルフィルタの比較
アナログフィルタ
ディジタルフィルタ
システム関数
H(s)
H(z)
周波数の範囲
( , )
[0, 2 ) or
周波数特性
非周期的
周期的
構
R, L, C や演算増幅器
計算機，DSP(Digital Signal Processor)
特性の変更
素子の交換
プログラムの修正
精度・再現性
素子のバラツキの影響あり
語長を長くすれば十分な精度達成可能
処理の遅延
なし（実時間処理）
演算器での遅延の累積が発生
成
( , ]
7.3 代表的なアナログフィルタ
IIR フィルタの間接設計法では，アナログフィルタに対して変数変換を施す．ここでは，代表
的なアナログフィルタを説明する．詳細は文献(17),(18)を参照．
1) バタワースフィルタ
バタワースフィルタ(Butterworth filter, BF)は次式の振幅 2 乗特性をもち，次数 N を大きくす
ることで，過渡域を急峻にできる．
1
2
H B ( ) H B ( ) H B ( )
1
c
2N
(7.6)
ここで， c は遮断角周波数である．HB () HB (s) s j で，H B () H B (s) s j であるので，
式(7.36)において，ω=-js（s=jωより）とすると，次式となる．
H B ( s ) H B ( s )
1
s
1 (1) N
c
(7.7)
2N
上式を満足し，かつ安定（極が左半平面内）な伝達関数 HB(s)を求める．N を偶数とすると，
式(7.7)の極μk は
s 2 N ( 1) N c2 N c2 N e
(7.8)
j ( 2 k 1)
信号処理 5 章-7 章.docx
- 53 -
大西
昇（名古屋大学）
より
k ce
j ( 2 k 1)
k 0 ,1, 2 , , 2 N 1
2N
(7.9)
となる．これらの極は，半径が c の円上で，等間隔 2 2 N に位置する．このうち，左半平面内
に存在する(安定な)極をλk とすると，HB(s)は
H B ( s)
k
(7.10)
N 1
(s
k
)
k 0
上式で， k
cN
とすると，HB(0)=1 となる．
2) チェビシェフフィルタ
チェビシェフフィルタ(Chebyshev filter, CF)は次式の振幅 2 乗特性をもつ．
H C ( )
2
1
1 e 2 Tn
c
(7.11)
2
e はゲインのリプル振幅を与えるパラメータ，TN(x)は N 次チェビシェフ多項式で，
T 1 ( x ) 1,
T 2 ( x ) 2 x 2 1,
T3 ( x) 4 x 3 3 x,
T 3 ( x ) 8 x 4 8 x 2 1,
のように，x の区間[-1,1]において，絶対値が１以下の振動的な多項式である．
3) バタワースフィルタとチェビシェフフィルタの比較
表 7.2 にバタワースフィルタとチェビシェフフィルタの比較をまとめる．
表 7.2
バタワースフィルタとチェビシェフフィルタの比較
バタワースフィルタ BF
チェビシェフフィルタ CF
通過帯域
平坦
等リップル
遮断特性
緩やか
鋭い
阻止帯域
リップルなく，振幅は単調に減少
同左
7.4 ディジタルフィルタの設計法
FIR フィルタと IIR フィルタの設計法を説明する．まず， FIR フィルタと IIR フィルタの比
較を，表 7.3 にまとめる．
表 7.3
FIR フィルタと IIR フィルタの比較
インパルス応答の継続長
FIR フィルタ
IIR フィルタ
有限
無限
-1
システム関数
z の多項式
z-1 の有理関数
安定性
常に安定
注意が必要
信号処理 5 章-7 章.docx
- 54 -
大西
昇（名古屋大学）
直線位相特性
完全に実現可能
近似的に実現可能
急峻な遮断特性の実現
高い次数が必要
低い次数で可能
量子化誤差拡大への対処
不必要
必要
1) FIR フィルタ
理想周波数特性 H I ( ) を近似するディジタル FIR フィルタ H ( ) を求める．
FIR フィルタの周波数特性は
N 1
H () H ( z ) z e j
h(k )e
jk
(7.12)
k 0
で与えられる．理想特性 H I ( ) を実現するように，h(k) (k=0,1,…,N-1 で実数値)を決める．
その方法には，フーリエ級数法（窓関数法）と周波数サンプリング法がある．以下で，これ
らを説明する．
a)フーリエ級数法（窓関数法）
H ( ) の理想特性 H I ( ) からのズレの 2 乗和は次式となる．
Je
2
H I () H () d
(7.13)
フーリエ級数法は上式を最小にするように，h(k) (k=0, 1, ..., N-1)を決める．Je を h(k)で
微分すると
J e
h ( k )
e
jk
( H I * ( ) H * ( )) d
e
jk
(7.14)
( H I ( ) H ( )) d 0
H* () H() (h(n)は実数)であるので，上式の第 1 項と第 2 項は等しく，
e jk (H I * () H * ())d e jk (H I () H ())d
である．したがって，(7.14)は次式となる．
e
jk
(H I () H ())d 0
これに， H ()
N 1
h(k )e
jk
(7.15)
を代入すると，次式を得る．
k 0
N 1
N 1
e jk ( H I () h(l )e jl )d e jk H I ()d
l 0
e
jk
h(l )e
l 0
j (l k )
d
(7.16)
H I ()d 2h(k ) 0
だから
h(k )
1
2
e
jk
H I ( ) d
(7.17)
信号処理 5 章-7 章.docx
- 55 -
大西
昇（名古屋大学）
式(7.17)は，H I ( ) をフーリエ級数展開した際のフーリエ係数が h(k)であることを意味する．
すなわち，次の関係が成立する．
H I ( )
h ( n )e
jn
I
(7.18)
n
1
hI ( n )
2
H
I
( ) e jn d
しかし，式(7.18)の h (n)は無限個あるので，次式のように，時間軸で窓関数 wt(n)を乗じ
て（このことから窓関数法とも呼ばれる）
，有限項数 N（=2K+1）の FIR にする．
h(n) hI (n)wt (n)
(7.19)
ここで，
1 n K
wt (n)
0 n K
(7.20)
である．さらに，因果的なフィルタにするために，K だけ時間シフトする．したがって，求め
るべきフィルタ係数は次式となる．
h (n ) hI ( n K ) wt (n K )
1
2
e
j (n K )
H I ( ) d
n 0, 1, , 2 K
(7.21)
式(7.20)の方形窓では，遮断周波数近傍でリップルが生ずるので，次式のハミング窓関数
(Hamming window function)を用いるとよい．その他の窓関数は付録１を参照．
2n
0.54 0.46 cos
wHam (n)
2K
0
n K
(7.22)
n K
b)周波数サンプリング法
周波数サンプリング法は，長さ N=2K+1 の系列 h(n)の離散フーリエ変換が，長さ N の離散系
列 H(k)となること（離散フーリエ変換対）を利用する（n, k は整数）
．すなわち，
K
H (k )
j
h ( n )e
2
kn
N
K k K
(7.23)
k K
理想の特性 H I ( ) を，周波数間隔 2 / N で標本化して
H I ( k ) H I ( )
k
2
N
(7.24)
K k K
を得る．h(n)(-K≦n≦K)と H(k) (-K≦k≦K)は離散フーリエ変換対であるので，長さ N のフィ
ルタ係数 h(n)（n=0,1,..,N-1）は， H I (k ) の離散逆フーリエ変換で求まる．すなわち，
h I (n)
1
N
K
k K
H I ( k )e
j
2
nk
N
K n K
(7.25)
ただし，因果性を満足するために，K だけ時間シフトする．すなわち，次式が求めるべき h(n)
である．
信号処理 5 章-7 章.docx
- 56 -
大西
h(n) h I (n K )
1
N
K
H
k K
I
( k )e
j
2
k (n K )
N
昇（名古屋大学）
(7.26)
0 n N 1
2) IIR フィルタ
IIR フィルタの設計法は，ディジタルフィルタを直接設計する方法（直接設計法）と，知見
の蓄積があるアナログフィルタの伝達関数に変換を施して設計する間接設計法とに大別され
る．本節では，間接設計法の，標準 z 変換法と双 1 次変換法を説明する．直接設計法について
の詳細は文献(15),(16)を参照．
アナログフィルタの伝達関数 Ha(s)（下付添え字の a はアナログを意味する）からディジタ
ルフィルタ H(z)を求める．
n
ss
H a ( s)
ki
i 1
i
として説明する．
a) 標準 z 変換法（インパルス不変法）
Ha(s)をラプラス逆変換して，インパルス応答 ha(t)を求める．
n
ha (t )
k e
i
si t
(7.27)
i 1
インパルス応答を標本化周期 T でサンプリングして ha(nT)(n=0,1,..)を得る．
n
ha (nT )
k e
i
si nT
(7.28)
i 1
これを z 変換して，次式のように H(z)を求める(フィルタは再帰型で実現)．
H ( z)
ha (kT ) z k
k 0
n
k i e si kT z k
k 0 i 1
n
ki
i 1
k 0
(e siT z 1 ) k
n
1 e
i 1
ki
siT
z 1
(7.29)
sT
Ha(s)と H(z)を比較すると，Ha(s)の極 si は，H(z)の極 zi e i になることが分かる．なお，
後述の式(7.30)の関係より，周波数特性を一致させるためには，式(7.29)を T 倍しなければ
ならない．
標準 z 変換法の短所を述べる．インパルス応答を周期 T でサンプリングするので，ディジ
タルフィルタとアナログフィルタの周波数特性の間に次式が成立する（式(4.7)参照）．
H ()
1
T
H
n
a
(
2n
)
T
( T )
(7.30)
信号処理 5 章-7 章.docx
- 57 -
大西
H a ( ) 0 ,
したがって，アナログフィルタが
昇（名古屋大学）
を満足しない場合に，ディジタル
T
フィルタの周波数特性にエイリアシングを生じる．
b) 双 1 次変換法
Ha(s)から
s
2 1 z 1
T 1 z 1
(7.31)
という双 1 次変換により， H(z)を求める．すなわち，
H (z) H a (
2 1 z 1
)
T 1 z 1
(7.32)
式(7.25)の双 1 次変換で， s j , z
2 tan 1
T
or
2
e
j
を代入すると
2
tan
T
2
(7.33)
上式をプロットすると図 7.4 となり， ( , )と ( , ) は 1 対 1 に対応することが
分かる．このように，エイリアシングを生じないことが双 1 次変換法の長所である．一方，
短所は，図 7.4 のように，アナログの角周波数 ω とディジタルの正規化角周波数 Ω との間
の非線形歪みである（直線関係ではない）
．
ω アナログフィルタの角周波数
T=1
2 tan( / 2)
4
3
歪み
Non-linear warping
2
1
π/８
0.5
1
1.5
2
Ω：ディジタルフィルタの角周波数
図 7.4
双 1 次変換でのωとΩ
以下で，双１次変換によるディジタルフィルタの設計手順を，バタワースフィルタを例に
して，説明する．ディジタルフィルタの設計仕様として，遮断周波数をΩc[rad]，阻止域端
周波数Ωs[rad]，阻止域減衰量 AS[dB]とし，T=1 とする．バタワースフィルタの振幅特性は
H a ( )
手順１
1
1
c
2N
である．
式(7.33)（ここでは，T=1 とする）により，アナログフィルタ Ha(ω)のωc, ωs
信号処理 5 章-7 章.docx
- 58 -
大西 昇（名古屋大学）
を求める（プリワーピング）．
手順２
As 20 log10 H a (s )
手順３
バタワース多項式（式(7.10)の分母で，例を表 7.4 に示す）の s を s/ ωc で置換
を満足する次数 N（整数）を決める．
え，Ha(s)を得る．
手順４
Ha (s) のｓに，式(7.31)（T=1 とする）を代入し，H(z)を得る
表 7.4 バタワース多項式（ωc=1 rad/s）
次数 N
バタワース多項式(Ha(s)の分母多項式)
１
(s+1)
２
(s2+√2s+1)
３
(s2+s+1)(s+1)
４
(s2+0.7653s+1)(s2+1.848s+1)
５
(s+1)(s2+0.6180s+1)(s2+1.618s+1)
６
(s2+0.5176s+1)(s2+√2s+1)(s2+1.932s+1)
c) 直接設計法
理想振幅特性を近似する N 次ディジタルバタワース低域フィルタの振幅特性は次式となる
（付録２参照）
．
1
H ( )
tan( / 2)
1
tan( c / 2)
2N
c : 遮断角周波数
(7.34)
上式から
H ( )
2
tan 2 N ( c / 2)
tan 2 N ( c / 2) tan 2 N ( / 2)
tan 2 N ( c / 2)
tan
2N
1 z 1
( c / 2) 1
1
1 z
2N
N
(7.35)
導出は省略するが，この式の分母多項式の根を求めることで，次式のシステム関数 H(z)が
得られる．
K (1 z 1 ) N
H ( z)
(1 pk z 1)
(7.36)
pk 1
システム関数 H(z)の零点は (1 z
1 N
) 0 の根で，安定な極 pk は
信号処理 5 章-7 章.docx
- 59 -
大西 昇（名古屋大学）
pk
1 qk
1 qk
(7.37)
ここで，
2k 1
c
tan 2 exp( j 2 N ), N even, k 0,...,2 N 1
qk
k
tan c exp( j ), N odd, k 0,...,2 N 1
2
N
(7.38)
設計仕様を，遮断角周波数Ωc [rad]，阻止域端角周波数Ωs [rad]，阻止域減衰量 As[dB]と
する．まず，次数 N は
As 20log10 H (s )
(7.39)
を満足する整数である．次にシステム関数の極 pk を，式(7.37)で求める．最後に， H (0) 1
となるように，K の値を決める．
7.5 直線位相 FIR フィルタ
直線位相では，入力波形が歪むことなく出力され，データ伝送系の重要な性質である．直線
位相には，Ⅰ型直線位相とⅡ型直線位相とがあり，どちらも FIR フィルタで実現できる．表
7.5 に，直線位相特性と，それを実現する FIR フィルタのインパルス応答が満足すべき性質を
まとめる．
Ⅰ型直線位相 FIR フィルタのシステム関数は，
N が奇数の場合
H ( ) e j ( N 1) / 2
( N 1) / 2
a
k 0
k
cos( k ),
(7.40)
a 0 h (( N 1) / 2), a k 2 h (( N 1) / 2 k )
N が偶数の場合
N/2
H () e j ( N 1) / 2 bk cos((k 1 / 2)),
k 1
bk 2h( N / 2 k )
(7.41)
となる．
表 7.5 直線位相特性（系列長 N の場合）
位相特性φ(Ω)
インパルス応答 h(k)の性質
Ⅰ型直線位相
-(N-1)Ω/2
偶対称；h(n)=h(N-1-n)
Ⅱ型直線位相
[π-(N-1)Ω]/2
奇対称；h(n)=-h(N-1-n)
信号処理 5 章-7 章.docx
- 60 -
大西 昇（名古屋大学）
7.6 低域通過以外の他のフィルタの導出
以上では，低域通過フィルタの設計だけを説明した．高域通過，帯域通過，帯域阻止などの
フィルタは，プロトタイプ低域通過フィルタに変数変換（周波数変換）を施して求める．
まず，アナログフィルタの場合を説明する．プロトタイプ低域通過フィルタを Hp(s)(遮断
周波数 ωc=1 rad./s)とし，変換後のフィルタを H(s')とすると，以下のような周波数変換に
より，所望のフィルタが求まる（図 7.5 参照）
．
低域通過フィルタ(遮断周波数が ωLc)
高域通過フィルタ(遮断周波数が ωHc)
H LP ( s ) H p ( s )
s
H HP ( s ) H p ( s)
s
(7.42)
s
Lc
(7.43)
Hc
s
帯域通過フィルタ(遮断周波数が ωLc と ωHc；ωLc<ωHc)
H BP ( s ) H p ( s )
s
s2 Lc Hc
s ( Hc Lc )
帯域阻止フィルタ(遮断周波数が ωLc と ωHc；ωLc<ωHc)
H BS ( s ) H p ( s )
s
s ( Hc Lc )
s2 Lc Hc
｜H(ω)|
｜H(ω)|
ωLC
ω
０
(a) 低域通過
ωHC
(b) 高域通過
図 7.5
ω
ω
０
ωLC
(7.45)
｜H(ω)|
｜H(ω)|
ω
０
(7.44)
ωHC
０
(c) 帯域通過
ωLC
ωHC
(d) 帯域阻止
各種フィルタの振幅特性
次にディジタルフィルタでは，プロトタイプ低域通過フィルタを Hp(z)(遮断周波数θc=1
rad.)とし，変換後のフィルタを H(z')とすると， H ( z ) H p ( z )
z 1 f ( z 1 )
で，所望のフィルタが
求める．周波数変換 z 1 f ( z 1 ) は次のようになる（文献(19)参照）．
低域通過フィルタ (遮断周波数が ωLc）
z 1
z1
Lc
Lc
, sin c
sin c
1
1 z
2
2
(7.46)
高域通過フィルタ(遮断周波数が ωHc）
z 1
z1
Lc
Lc
, cos c
cos c
1
1 z
2
2
(7.47)
帯域通過フィルタ(遮断周波数が ωLc と ωHc；ωLc<ωHc）
信号処理 5 章-7 章.docx
- 61 -
大西 昇（名古屋大学）
2 k 1 k 1
z
k 1
k 1 ,
z 1
k 1 2 2 k 1
z
z 1
k 1
k 1
Lc
Lc
cos Hc
cos Hc
,
2
2
z 2
(7.48)
Lc
k cot Hc
tan c
2
2
帯域阻止フィルタ(遮断周波数が ωLc と ωHc；ωLc<ωHc）
2
1 k
z 1
1 k
1 k ,
z 1
1 k 2
2
z
z 1 1
1 k
1 k
Lc
Lc
cos Hc
cos Hc
,
2
2
z 2
(7.49)
Lc
k tan Hc
tan c
2
2
信号処理 5 章-7 章.docx
- 62 -
大西
昇（名古屋大学）
8. 高速フーリエ変換(FFT)
離散フーリエ変換を効率的に計算できるアルゴリズムが，高速フーリエ変換(FFT)である．FFT
アルゴリズムの導出と原理を説明し，畳み込み，相関関数，スペクトルの計算など，信号処理へ
の FFT の応用を示す．
8.1 高速フーリエ変換とは
離散フーリエ変換を高速に実行するアルゴリズムが，高速フーリエ変換 (Fast Fourier
Transform, FFT)で，1965 年に Cooley & Tukey によって，発表された．
次式の離散フーリエ変換を計算するためには，N2 回の乗算が必要である．
N 1
X (k )
x ( n) W
0 k N 1
kn
N
(8.1)
n 0
しかし，FFT では
データ数 N が N p1 p 2 p r
データ数 N が N p r
の時
乗算回数 M は
の時
乗算回数 M は
M N ( p1 p 2 p r )
M pN log p N
となる．特に， N p r の場合，直接計算と比べた乗算回数の比率は
M
N
2
p
log 10
log 10 N
p
N
p
で，p=3 の時，
log 10 p
なる p=2 では，
p
log 10 p
(8.2)
6.29 で最良となる（文献(20)を参照）
．しかし，アルゴリズムが簡素に
6.64 と p=2 と僅かな差である．したがって，FFT と言えば，データ数 N
は２のべき乗（p=2）とするのが一般的である．
8.2 FFT アルゴリズムの導出
データ数 N=2m とする．離散フーリエ変換の式
X (k )
N 1
n0
k 0,1, , N 1
x(n) W Nkn
WN e
j
2
N
(8.3)
を，次のように，x(n)の偶数番目と奇数番目とに分ける．
N / 2 1
X (k )
N / 2 1
x(2r ) W N2 rk
r 0
x(2r 1) W N( 2 r 1) k
r 0
N / 2 1
N / 2 1
x(2r )(W N2 ) rk W Nk
r 0
x(2r 1) (W
2 rk
N)
(8.4)
r 0
2
j 2 /( N / 2)
WN / 2 を用いると，上式は，
y ( n) x( 2n), z ( n) x ( 2n 1) ( n 0,1, , N / 2 1) とし，WN e
X (k )
N / 2 1
y (n) (W
n0
N /2
) nk WNk
N / 2 1
z (n) (W
n0
N/2
) nk Y (k ) WNk Z ( k )
k 0,1, , N 1
(8.5)
となる．すなわち，X(k)は，偶数番目のデータ x(2n)の DFT である Y(k)と，奇数番目の DFT で
ある Z(k)とから求まる．Y(k)と Z(k)は周期 N/2 の周期系列であるので，
N
N
(8.6)
) Z (k )
k 0,1, ,
1
2
2
となる．したがって，式(8.5)により，すべての k (k 0,1, , N 1) について X(k)を計算できる．
Y (k
N
) Y (k ),
2
Z (k
信号処理 8 章-10 章.docx
- 63 -
大西
昇（名古屋大学）
Y(k)と Z(k) (k 0,1, , N / 2 1) それぞれのデータ系列について，上記と同様のことを繰り返
す．ここでは，Y(k)についてのみ説明する．次式のように，y(n)を偶数番目と奇数番目とに分
け，変形すると，次式を得る．
N / 2 1
Y (k )
N / 4 1
y (n) (W N / 2 ) nk
n 0
N / 4 1
y(2r 1) W
r 0
N / 4 1
y(2r )(W
2
rk
N /2)
y(2r )(W
N /4)
W Nk / 2
r 0
N / 4 1
N / 4 1
y (2r ) W N2 rk/ 2
( 2 r 1) k
N /2
r 0
y(2r 1) (W
2
rk
N /2)
y(2r 1) (W
N /4)
(8.7)
r 0
N / 4 1
rk
W Nk / 2
r 0
W N2 / 2 W N / 4
rk
r 0
すなわち，n が偶数の y(n) (y(n)= x(2k)であるので x(4k))の DFT と，n が奇数の y(n) (x(4k
＋2))の DFT とから，Y(k)は計算される（n の初期値はゼロであることに注意）
．そして，分割
された系列の分割を長さが２になるまで繰り返す．
このように，FFT は，長さ N の系列を，長さ N/2 の系列への分割（分割統治）の繰り返しと，
W Nk
e
j
2
k
N
の周期性を利用したアルゴリズムである．
以下で，N=4 の場合について，FFT アルゴリズムの導出を説明する．
1
X (k )
1
x(2r )(W42 ) rk W4k
r 0
x(2r 1) (W
k 0,1, 2, 3
2 rk
4 )
r 0
より， Y (k )
1
r 0
1
x(2r )(W42 ) rk , Z (k )
r 0
x(2r 1) (W42 ) rk
とすると
X (k ) Y (k ) W4k Z (k ) k 0,1, 2, 3
また，
Y ( k 2) Y ( k ),
で， W4 e
j 2 / 4
Z ( k 2) Z ( k )
k 0,1
あるので
k 0
X (0) Y (0) W40 Z (0) Y (0) Z (0)
k 1
X (1) Y (1) W41 Z (1) Y (1) jZ (1)
k 2
X (2) Y (0) W42 Z (0) Y (0) W40 Z (0) Y (0) Z (0)
k 3
X (3) Y (1) W43 Z (1) Y (1) W41 Z (1) Y (1) jZ (1)
次に，Y(k), Z(k) (k=0,1)は
1
Y (k )
1
x(2r )(W42 ) rk , Z (k )
r 0
x(2r 1) (W
2 rk
4 )
k 0,1
r 0
で， W2 e j 2 / 2 あるので，
k 0 Y (0) x(0)(W21 ) 0 x(2)(W21 ) 0 x(0) x(2)
Z (0) x (1)(W21 ) 0 x(3)(W21 ) 0 x(1) x(3)
k 1 Y (1) x(0)(W21 ) 0 x(2)(W21 )1 x(0) x(2)
Z (1) x(1)(W21 ) 0 x (3)(W21 )1 x(1) x(3)
信号処理 8 章-10 章.docx
- 64 -
大西
昇（名古屋大学）
以上の計算処理は図 8.1 のように表わされる．図で，入力 x(n)の順序がビットリバーサ
ル(次節で説明)になっていることに注意．また，図は，図 8.2 のような構造の繰り返しで
構成されている．この基本となる計算は，羽を広げた蝶の形に似ているので，バタフライ
計算（Butterfly；蝶，X 型の支柱）と呼ばれる．
x(0)
ビットリバーサル
x(2)
Y(0)
X(0)
Y(1)
X(1)
-1
W40
Z(0)
x(1)
-1
X(2)
-1
X(3)
-1
図 8.2
W41
-1
x(3)
Z(1)
FFT 計算のフロー（N=4）
図 8.1
W41 e
j
2
4
バタフライ演算
j
8.3 ビットリバーサル
ビットリバーサル(bit-reversal)とは，データ番号を表わすビットの並びを，次式のように，逆
順にすることである．
bk bk 1 b2b1
b1b2 bk 1bk
(8.8)
表 8.1 に，元のデータ番号 n のビットを逆順にした番号 n’（N=8 の場合）を示す．
表 8.1
n
ビットリバーサルの説明（N=8）
0
1
2
3
4
5
6
7
n のビット
000
001
010
011
100
101
110
111
ビット逆順
000
100
010
110
001
101
011
111
0
4
2
6
1
5
3
7
n’
8.4 FFT の応用
1) 畳み込みの計算
a) 循環畳み込み
周期Ｎの 2 つの周期系列 ~x (n) , ~y (n) の循環畳み込み(5.3c)参照)は次式で計算される．
N 1
u ( n) ~
x (k ) ~
y (n k ) PN (n) n 0,1, , N 1
k 0
上式のｚ変換 U(k)は， x(n)のｚ変換
である X(k)と， y(n)のｚ変換である
x(n)
(8.9)
ＤＦＴ
~
x ( n) * ~
y ( n)
X (k )
Y(k)との積，すなわち，U(k)=X(k)Y(k)
となる．したがって，U(k)は図 8.3 の
ように計算できる． 離散フーリエ変換
×
y (n)
ＤＦＴ
ＩＤＦＴ
X (k )Y (k )
Y (k )
信号処理 8 章-10 章.docx
- 65 - 図 8.3
循環畳み込みの FFT による計算
大西
昇（名古屋大学）
(DFT)を 2 回実行し，その後，離散逆フーリエ変換(IDFT)と，一見複雑に見える．しかし，
DFT と IDFT の計算に FFT を使用すると，計算量は O(Ｎ２)から O(Ｎlog２Ｎ)となり，計算量
をかなり少なくできる．ここで，O(x)（オーダｘ）は，計算量が x の定数倍で押さえられる
ことを意味する．
b) 離散畳み込み
2 つの系列 x(n), y(n)の離散畳み込みは（系列の長さを無限とすると）
w(n)
x(k ) y (n k )
(8.10)
k
となる．
この離散畳み込み w(n)と，
式(8.9)の循環畳み込み u(n)との間には次式が成立する．
u ( n)
w(n kN )
0 n N 1
(8.11)
k
式(8.10)で系列 x(n), y(n)の長さを Lx, Ly とすると，w(n)の長さは Lx+Ly-1 となる．周期 N
が， N L x L y 1 を満足するなら，その区間内で，u(n)=w(n)となる．すなわち，系列 x(n),
y(n)にゼロを追加して，共に長さ N の系列とした後，循環畳み込み u(n)を計算し，0 n N 1
を取り出せば，離散畳み込み w(n)が得られる．FFT を循環畳み込みの計算に利用することで，
離散畳み込みの計算量を低減できる．
2) 相関関数の計算(文献(21))
ウィナー・ヒンチンの定理(3.7 参照)に
x (n)
よれば，x(t)の自己相関関数のフーリエ変
2
×
ＤＦＴ
換 は ， X( ω ) の パ ワ ー ス ペ ク ト ル
xx ( ) X ( ) になる．これを離散時間
X (k ) X * (k )
X (k )
共役
図 8.4
X * (k )
ＩＤＦＴ
N rxx (k )
自己相関関数の FFT による計算
で考える．離散系列 x(n)(n=0,1,..,N-1)
の自己相関関数は次式で計算される．
rxx (k ) E x(n) x(n k )
1
N
N 1
x ( n) x ( n k )
(8.12)
n0
上式で，x(n)を周期 N の周期系列に拡張した ~
x ( n ) を用いて循環型相関関数を計算する．
rxx (k ) の離散時間フーリエ変換は，
R xx (k )
1
X (k ) X * (k )
N
(8.13)
となる．すなわち，x(n)の離散時間フーリエ変換 X(k)のパワーの 1/N 倍になる．そこで，
図 8.4 のように，FFT を用いて，x(n)の DFT を求め，|X(k)|2 を計算した後，それを IDFT す
ることで，自己相関関数 rxx(k)の N 倍が得られる．
x(n)と y(n)の相互相関関数の場合，図において，y(n)の DFT を追加し，Y(k)の共役を求
めるような変更をすればよい．
信号処理 8 章-10 章.docx
- 66 -
大西
昇（名古屋大学）
3) 確定信号のスペクトル推定
確定信号 x(t)のスペクトル X(ω)を，FFT を用いて計算することを考える．図 8.5 のよう
に，持続時間が有限の時間制限された信号 x(t)を周期 Tp で拡張して周期信号を作成する．
つぎに，図(b)の周期信号 xc(t)を，時間間隔 T で標本化する．この離散フーリエ変換は，図
(c)のように，離散周期系列 Xc(k)となる．すなわち，次式が成立する．
1
T
X c (k )
X (k m
m
s
(8.14)
)
信号 x(t)のスペクトル X(ω)が，ωs/2 以上の成分がゼロであるように，帯域制限されてい
るならば，スペクトル X(ω)の離散的な周波数での推定 Xˆ (k ) が次式で得られる．
Xˆ (k ) TX c (k )
(8.15)
もし，帯域制限が完全でなく，スペクトルの折り返し(エイリアシング)がある場合は，標
本間隔 T を小さくして，折り返しがなくなるまでωs を大きくする．なお，T を小さくする
際，標本数 N（=Tp /T）を増やして，周期 Tp が変わらないようにする(その結果，時間信号
x(t)の重なりを避けられる)．
Tp
x(t)
ωs
xc(t)
Xc(k)
t
t
T
(a)有限持続信号
s
(b)周期信号と離散時間信号
図 8.5
ω
2
2
Tp
(c)離散時間周期信号のスペクトル
有限持続信号のスペクトル
しかし， x(t)が時間制限されていない場合が多い.その場合は， x(t)に窓関数(window
function)をかけて時間制限する.すなわち
x (t ) w(t ) x(t )
(8.16)
このフーリエ変換は
X ( ) W ( ) * X ( )
(8.17)
と，周波数領域での畳み込みになる．したがって，得られるスペクトルは，元のスペクトル
X(ω)を，窓関数のスペクトル W(ω)によって平滑化(smoothing)したものになる．
4) 不規則信号のパワースペクトル推定
平均値がゼロの定常不規則信号 x(n) (n=0,1,..,N-1)の自己相関関数 rxx(k)(-(N-1)≦k≦
(N-1))の離散時間フーリエ変換(5.2 参照)は次式となる．
R xx ()
N 1
k
xx
k ( N 1)
rxx (k )e jk
r
(k )e jk
(8.18)
上式の 2 つ目の等号では，rxx(k)の系列長が有限であることを用いている．式(8.18)はつぎ
信号処理 8 章-10 章.docx
- 67 -
大西
昇（名古屋大学）
のように書き換えることができる．
R xx ()
1
2
X ()
N
( )
(8.19)
すなわち，自己相関関数 rxx(k)のフーリエ変換 R xx () を N 倍すると，信号 x(n)のパワース
ペクトル X () が得られる． R xx () をピリオドグラム(periodgram)と言う．式(8.18)のフ
2
ーリエ変換は FFT で計算でき，さらに， 2)で述べたように,自己相関関数の計算にも FFT を
使用できる．
信号処理 8 章-10 章.docx
- 68 -
大西
昇（名古屋大学）
9. 線形予測法
不規則信号の解析方法として線形予測法がある．これは，入力が白色雑音である動的システ
ムの出力として，不規則信号をモデル化するもので，音声信号の処理などに用いられる．予
測モデルと予測係数および最適予測次数の求め方，および予測係数と数学的に等価な偏自己
相関係数を説明する．
9.1 時系列の解析
離散時間不規則信号のサンプル値系列を時系列(time series)と言う．時系列の解析方法とし
ては次の２つがある．
方法１
時系列の統計的性質を相関関数やスペクトルで記述する．
方法２
時系列 s(n)を動的なシステムの出力としてモデル化する（図 9.1 参照）
．
図 9.1
モデルによる信号の表現
方法２の代表的なものが線形予測法(linear prediction method)である．本方法では，観測信号
s(n)は次式でモデル化される．
q
p
s ( n)
a k s(n k )
k 1
b e( n l )
(9.1)
l
l 0
ここで，e(n)はシステムの入力で，白色雑音を仮定する．上式の両辺を z 変換すると
S ( z) H ( z)E( z)
(9.2)
となる．ここで，H(z)はシステム関数で，次式で与えられる．
q
b z
l
l
l 0
p
H ( z)
1
a
k
(9.3)
z
k
k 1
9.2 線形予測モデル
式（9.1）は，極零モデル(pole-zero model)（自己回帰－移動平均モデル）という．これの特
別な場合として，つぎの２つがある．
全零モデル(all-zero model)（移動平均(moving average)モデル；非再帰型；FIR）：ak=0,
k=1,2,..,p
全極モデル(all-pole model)（自己回帰(auto-regressive)モデル；再帰型）
：bl=0, l=1,2, ..,q
以下では，次式で表される全極モデルに限定して，説明する．
p
s ( n)
a
k s ( n k ) e( n )
(9.4)
k 1
信号処理 8 章-10 章.docx
- 69 -
大西
昇（名古屋大学）
上式右辺の第 1 項を，s(n)の予測（推定値） sˆ( n) と呼ぶ．すなわち，予測 sˆ( n) は
p
sˆ(n) a k s(n k ) a T s (n)
(9.5)
k 1
こ こ で ， ak は 予 測 係 数 (predictor coefficients) ， p は 予 測 次 数 で ， a (a1 , a 2 , , a p ) T
, s ( n) ( s ( n 1), s (n 2), , s ( n p )) T ある．すると， e(n)は，s(n)と sˆ( n) の差であり，
これを予測残差(residual error)と呼ぶ．予測残差の性質（詳細は文献(23)を参照）として，つ
ぎの 2 点を仮定する．
・予測残差 e(n)と信号 s(n)は無相関である．
・次数ｐが十分大きい時，予測残差の自己相関関数はインパルスδ(n)（白色雑音）になる．
9.3 予測係数の決定
式(9.4)の未知パラメータは予測次数 p(その最適化は 9.7 で説明)と予測係数である．
まず，
予測係数の求め方を説明する．
予測残差 e(n)の２乗の総和 E は次式となる．
E
e 2 ( n)
n
s(n) sˆ(n)2
n
s ( n)
n
2
a k s(n k )
s ( n) a T s ( n )
k 1
n
p
2
(9.6)
E を最小にする ak (k=1,2,..,p)を求める（これを最小２乗法という）．E は ak の２次形式で
あるので，求めるべき予測係数 ak (k=1,2,..,p)は
E
E en
2 {s (n)
a m
e(n) a m
n
2[ s (n) s (n m)
n
p
a
k 1
k
s (n k )}s (n m)
(9.7)
p
a { s(n k )s(n m)}] 0
k 1
k
1 m p
n
で構成される連立方程式（正規方程式(normal equation)あるいは，ユール・ウォーカー
(Yule-Walker)方程式と呼ばれる）の解である．
式(9.7)を，s(n)の自己相関関数 rss ( k )
s(n k )s(n) を使って書き直すと，
n
p
rss (m)
a r
k ss ( m k )
0
1 m p
(9.8)
k 1
となる．ここで，自己相関関数では次式が成立することに注意．
rss ( m) rss (m)
(9.9)
式(9.8)を行列表現すれば，
rss ( p 1) a1
rss (1)
rss (0)
rss (1)
r (1)
(
0
)
(
2
)
r
r
p
a
ss
ss
ss
2
rss (2)
r ( p 1) r ( p 2)
rss (0) a p
rss ( p )
ss
ss
(9.10)
上式で，左辺の行列は，対角要素および副対角要素が同じで，このような行列を Toeplitz
行列（方程式の解を高速に得るアルゴリズムが存在する）と言う．
信号処理 8 章-10 章.docx
- 70 -
大西
昇（名古屋大学）
9.4 偏自己相関係数
予測係数は，式(9.10)の連立方程式の解として
求めることができる．しかし，予測係数は，次数
p が変わると，新たに計算しなおさなければなら
ない．すなわち，次数(p-1)での予測係数に基づ
いて，次数 p での係数を計算することができない
(その結果，計算量が大となる)．
本 節 で 説 明 す る 偏 自 己 相 関 係 数 (Partial
Autocorrelation Coefficient : PARCOR)は，予測係
数と数学的には等価であるが，予測次数に影響さ
れないという，予測係数にない優れた利点がある．
まず，偏自己相関係数の求め方を説明する．
e (fn 1)
( 1)
ebn
図 9.2
偏自己相関係数の求め方
図 9.2 のように，時間軸にそった前向き予測と,
時間軸を逆行する後向き予測を考える．前向きの予測 sˆ(n) を
1
sˆ(n) i( 1) s (n i )
(9.11)
i 1
後向きの予測 sˆ(n ) を
1
sˆ(n ) i( 1) s (n i )
(9.12)
i 1
( 1)
と表す．式(9.11)，(9.12) による予測残差をそれぞれ， e fn
( 1)
， ebn
とすると，τ次の偏自
( 1)
己相関係数は次式のように，(τ-1)次の前向予測残差 e (fn 1) と後向予測残差 ebn
の相関係数
k として求まる．
k
ここで，
( 1)
e (fn 1) ebn
e ( 1) 2 e ( 1) 2
bn
fn
(9.13)
1/ 2
( )
( )
は時間平均を表す．そして，τ 次の予測残差 e fn ， ebn
は，(τ-1)次の残差と
τ次の偏自己相関係数を使って，次式で計算される．
( 1)
e (fn ) e (fn 1) ebn
(9.14)
( )
ebn
eb(( n1)1) e (f( n 1)1)
(9.15)
次に，(9.13)，(9.14)および(9.15)の導出を簡単に説明する(詳細は文献(24))． 前向きの
予測と残差を用いると，観測信号 s(n)は次式となる．
信号処理 8 章-10 章.docx
- 71 -
大西
昇（名古屋大学）
1
s (n) i( 1) s (n i ) e (fn 1)
(9.16)
i 1
同様に，後向きの予測 s(n-τ)は
1
( 1)
s (n ) i( 1) s (n i ) ebn
(9.17)
i 1
となる．定常信号では時間方向に対して相関特性は不変であるので
i( 1) (i 1) , 1 i 1
(9.18)
( 1) 2
{e (fn 1) }2 {ebn
}
(9.19)
( )
が成立する．τ次の前向き予測残差 e fn は，次のように展開できる．
e (fn ) s(n)
s ( n)
i( ) s(n i) s(n)
i 1
1
s ( n)
1
i 1
( )
i
s (n i ) ( ) s(n )
( 1)
s (n i ) ( ) { i( 1) s (n i ) ebn
}
( 1)
i
i 1
1
( )
i
i 1
1
(9.20)
i 1
( 1)
s(n i ) ( ) ebn
ここで，
i( 1) i( ) ( ) i( 1) ,
式(9.20)の残差の総和 E
1 i 1
{e }
n
( ) 2
fn
(9.21)
を最小にする i( 1) , i( ) (1 i 1) は次式の解であ
る．
E
E
i( 1)
( )
1 i 1
0,
(9.22)
0
(9.23)
( 1)
仮定より， ebn
, s ( n) は無相関であるので，式(9.22)は
E
(j 1)
2{s(n)
n
1
n
i 1
2{s(n)
1
i 1
( 1)
i
( 1)
i
( 1)
s (n i ) ( ) ebn
}s (n j )
(9.24)
s (n i )}s (n j ) 0
となる．上式は(τ-1)次の正規方程式そのものであるので
i( 1) i( 1) ,
1 i 1
(9.25)
となる．これを，式(9.20)に代入すると次式となる．
信号処理 8 章-10 章.docx
- 72 -
大西
e (fn ) s(n)
s ( n)
1
i 1
1
( 1)
i
i 1
( 1)
i
昇（名古屋大学）
( 1)
s (n i ) ( ) ebn
( 1)
s(n i ) ( ) ebn
(9.26)
( 1)
e (fn 1) ( ) ebn
( ) 2
} と，総和を時間平均に置き換える），
上式を式(9.23)に代入すると（ただし， E {ebn
E
( )
{e (fn ) }2
( 1)
( 1)
2ebn
{e (fn 1) ( ) ebn
}0
( )
(9.27)
( 1) 2
{e (fn 1) }2 {ebn
} （式(9.19)）であるので，上式より
( )
( 1)
e (fn 1) ebn
[e
( 1) 2
fn
e
(9.28)
( 1) 2 1 2
bn
]
したがって，式(9.26)のτ次の前向き予測残差は，次式のように，(τ-1)次の予測残差と偏
自己相関係数で求まる．
( 1)
( 1)
e (fn ) e (fn 1) ( ) ebn
e (fn 1) k ebn
(9.29)
同様にして，後向予測残差について次式が得られる．
1( ) k
(9.30)
( )
ebn
eb(( n1)1) k e (f( n 1)1)
(9.31)
よって，式(9.13)，(9.14)および(9.15)が導出された．
9.5 偏自己相関係数から予測係数
偏自己相関係数から予測係数を求める方法を説明する．式(9.21)の
i( 1) i( ) ( ) i( 1) ,
( 1)
に， i
1 i 1
i( 1) (1 i 1) ， ( ) k ， i( 1) (i 1) (1 i 1)
を代入
すると，次式となる．
i( ) i( 1) k i( 1) i( 1) k (i 1) ,
1 i 1
(9.32)
すなわち，τ次の予測係数 i( ) は，(τ-1)次の予測係数 i( 1) とτ次の PARCOR 係数 k から計
算できる．
( ) 2
次に，残差の総和 E {ebn
}
について成立する漸化式を導出する．
Eτに式(9.15)を代入すると，次式となる．
信号処理 8 章-10 章.docx
- 73 -
大西
昇（名古屋大学）
( ) 2
E {ebn
} {eb(( n1)1) k e (f(n 1)1) }2
(9.33)
{eb(( n1)1) }2 2k eb(( n1)1) e (f(n 1)1) k {e (f(n 1)1) }2
2
上式に，式(9.28)の関係
e (f(n 1)1) eb(( n1)1) [e (f(n 1)1) eb(( n1)1) ]1 2
2
2
を代入すると，式(9.33)は次式となる．
( ) 2
{ebn
} {eb(( n1)1) }2 2k2 [e (f(n 1)1) eb(( n1)1) ]1 2 k {e (f(n 1)1) }2
2
2
2
(9.34)
( 1) 2
{e (fn 1) }2 {ebn
} の関係を，上式右辺の第２項に適用すると，
( ) 2
{ebn
} {eb(( n1)1) }2 2k2 [eb(( n1)1) eb(( n1)1) ]1 2 k {eb(( n1)1) }2
2
( 1) 2
b ( n 1)
{e
( 1) 2
b ( n 1)
} k {e
2
2
( 1) 2
b ( n 1)
} {e
2
(9.35)
} {1 k }
2
( ) 2
} は，(τ-1)次の残差の総和 E 1 {eb(( n1)1) }2
となる．上式より，τ次の残差の総和 E {ebn
とτ次の PARCOR 係数 k から計算できることが分かる．また，式(9.13)の右辺の分子，残差の
積の時間平均は
( 1)
e (fn 1) ebn
s ( n)
s ( n)
rss ( )
1
1
i 1
i 1
1
( 1)
i
i 1
( 1)
i
s (n i ) s (n )
1
i 1
( 1)
i
s (n i )
s (n i ) s (n )
(9.36)
r ( i )
( 1)
i
ss
となり，信号の自己相関関数 rss (i ) と予測係数から計算できる．
9.6 Levinson-Durbin の方法
Levinson-Durbin の方法は，予測次数 i を順次大きくしながら，偏自己相関係数 ki と予測
係数 aj(i) (j=1,2,..,i)を計算する方法である．自己相関関数を rss (i )
n
s ( n ) s ( n i ) ，i
(i )
次の偏自己相関係数を ki，予測係数を a j として，Levinson-Durbin の方法を以下に示す．
Step 1
Step 2
W0 rss (1),
E 0 rss (0),
k i Wi 1 Ei 1 ,
a i(i ) k i ,
i 1
E i E i 1 (1 k i2 )
a (ji ) a (ji 1) k i ai(i j1) (1 j i 1)
信号処理 8 章-10 章.docx
- 74 -
大西
Step 3
昇（名古屋大学）
i = p （p はあらかじめ決めた最高次数）なら終了．さもなければ(i < p)
Wi rss (i 1)
i
a
k 1
r (i 1 k )
(i )
k ss
i = i +1
Step 2 へ
9.7 予測次数の最適化（AIC 基準）
予測次数を大きくすれば，予測残差を小さくすることができる．しかし，これは手元にあ
る標本系列への過剰な当てはめで，不規則信号の他の標本系列に対しては残差が多くなるこ
とが予想される．そこで，線形予測モデルの最適次数を求める基準が必要である．
以下で，よく使用されている AIC 基準 (Akaike information criterion)を説明する．最適次数
は，次式の AIC 基準を最小にする p である．
AIC ( p ) log e
Ep
E0
2p
Ne
(9.37)
上式で， p は次数，Ne は実効データ数(=c N ; N はデータ数，c1は補正係数)，Ep は予測
残差の２乗和で
e
2
(n) ，そして E0 は系列のパワー rss (0)
n
n
s 2 (n) である．式(9.37)の
右辺において，次数 p の増大とともに第 1 項は減少するのに対し，第 2 項は増加する．した
がって,最小点 p は必ず存在する．なお，極小値が複数存在するので，次数の上限を 3 N と
して，AIC の最小値を求めると良い（詳細は文献(24)参照）
．
9.8 適用例（音声のピッチとフォルマント）
図 9.3 (a)に母音生成のモデルを示す．図で，音源は，声帯の振動により生ずる断続波で，
三角波の繰り返しで近似される．声道とは，舌，歯，唇，顎などによって形作られた音響的
な管で，声道の形を種々変化させ特定の形にすることを調音という．調音された音波は唇か
ら空間に放射され，音声として知覚される．音源の繰り返し周期をピッチ(pitch)周期と言い，
声道の伝達関数のピーク周波数（共振周波数）をフォルマント(formant)（低い方から第１，
第２フォルマント）と言う．ピッチ周波数（基本周波数）で音の高低（男性の低い声は 50Hz，
女性の高い声は 500Hz）が決まる．一方，フォルマントは，母音を特徴づける優勢な周波数
成分であるので，母音の識別に利用できる．
つぎに音声信号 s(n)から，ピッチ周期とフォルマントを求める方法を簡単に説明する．線
形予測モデルでは，時刻 n の信号 s(n)は，過去の観測信号 s(n-k) (k=1,...,p)の重み付の和
で表される予測値 sˆ(n) と予測残差 e(n)で，次式のように表される．
p
s (n) sˆ(n) e(n)
a
k s ( n k ) e( n )
k 1
上式を，e(n)が入力，s(n)が出力となるシステムとして表現すると，図 9.3 (b)のようになる．
点線で囲った部分が声道に対応する．音声信号 s(n)から，AIC 基準を用い，Levinson-Durbin
数 c はデータウィンドウと矩形窓のエネルギー比で，Hamming 窓の場合 c=0.4 となる．
1
信号処理 8 章-10 章.docx
- 75 -
大西
昇（名古屋大学）
の方法で，最適予測次数と予測係数を求める．音源 e(n)は三角波の繰り返しで，周期関数で
あるので，予測残差 e(n)の自己相関関数をもとめ，そのピークの間隔から，ピッチ周期が求
まる．つぎに，フォルマント（共振周波数）はシステム関数
H ( z)
S ( z)
1
E ( z ) 1 P( z )
p
where P( z )
a
kz
k
(9.38)
k 1
の極に対応するので，H(z)の振幅特性において，ピークを与える周波数をフォルマントとす
ればよい．
音声
音源
(a)
声帯
声
道
H ( )
ω
s (n)
e(n)
(b)
＋
P( z ) k 1 ak z k
P
sˆ(n)
図 9.3
母音の生成モデル(a)と線形予測器を用いたモデル(b)
図 9.4 は，
「あ」の音声波形(男性)である（サンプリング周波数８kHz，量子化ビット数 16bits）
．
予測モデルの次数が 20 の時の結果を以下に示す．予測残差の自己相関関数は図 9.5 のようにな
り，ピークの間隔から基本周波数は 116.4Hz となる．予測係数から(9.37)のシステム関数 H(z)
を求め，振幅スペクトルをプロットしたのが図 9.6 である．約 800Hz と 1.1kＨｚおよび 1.9kHz
にピーク（フォルマント）が見られる2．
2
「あ」の場合，第１フォルマントは 0.5-0.8kHz，第２は 1-1.3kHz，第３は 2.2-2.9kHz である．
信号処理 8 章-10 章.docx
- 76 -
大西
図 9.4
昇（名古屋大学）
「あ」の音声波形(男性)
図 9.5
予測残差の自己相関関数
10
5
0
-5
-10
-15
-20
-25
0
1000
図 9.6
2000
3000
4000
5000
システム関数の振幅特性
横軸は周波数，縦軸は振幅スペクトル(dB)
信号処理 8 章-10 章.docx
- 77 -
大西
昇（名古屋大学）
10. 適応信号処理
システムの入力や条件・環境の変化に応じて，システムの特性を適応的に変化できる機能を備え
た信号処理を適応信号処理(adaptive signal processing)と言う．本章では，適応信号処理の例，FIR
型適応フィルタ，代表的な適応アルゴリズムを説明する．
10.1 適応システム
従来のシステムでは，入力信号や環境の変化
入
力
応 答
可変システム
に対してシステムのパラメータ値は固定であ
る．したがって，システムの入力や条件・環境
参照信号
が，設計時のそれから変化すると，システムは
適応アルゴリズム
望ましい応答を出力できない．そこで，そのよ
誤
うな変化に応じて，システムの特性を適応的に
差
図 10.1 適応信号処理システムのブロック図
変化できる機能が必要であり，そのような機能
を備えたシステムが適応システムである．
環境変化に応じてシステムの特性を変化させる適応システムでは，図 10.1 のように，参照
信号（理想応答とも言う）とシステムの応答との誤差を最小にするように，可変システムの
パラメータを変化させる．その方法を適応アルゴリズム(adaptive algorithm)と言う．
10.2 適応信号処理の例
1) 自動波形等化器
データ伝送系では，受信信号が送信信号と一致することが必要である．しかし，図 10.2
に示すように，伝送路の特性や伝送路に混入する雑音のため，受信信号は送信信号と一致し
ない．すなわち，受信信号 y(k)は
y (k )
a(n)w(k n) v(n) a(k )w(0)
n
a(n)w(k n) v(n)
(10.1)
n , n k
となる．右辺の第 2 項は符号間干渉で，第 3 項は雑音である．そのために，これらの歪みを
ゼロにするために，波形等化器3 (equalizer)が必要となる．
雑音 v(n)
a(n)
伝送路
w(n)
y (n)
＋
自動等化器
h(n)
判
aˆ (n)
定
図 10.2 データ伝送系と自動波形等化器
符号間干渉などをゼロにするように波形等化器を設計しても，経年変化や温度などのため
等化器 システムや構成成分の好ましくない振幅や位相の周波数応答を補正する目的で設
計された回路網．
3
信号処理 8 章-10 章.docx
- 78 -
大西
伝送路特性は変化する．したがって，それらの変動に
昇（名古屋大学）
SA
適応できる自動波形等化器が必要である．自動波形等
SA
化器とは，伝送系の変動に応じて，等化器特性を適応
的に調整し，符号間干渉を抑圧する器械である．
SB
SB E
2) 伝送系のエコーキャンセラー
4 線式と 2 線式回線を結合する伝送系では，結合部
でエコーを発生する．図 10.3(a)で，話者 A は，話者
B の声 SB と，自身の声 SA のエコーE とを聞く．このよ
うなエコーを除去するため，図(b)のように，エコー
キャンセラーを挿入する．エコーキャンセラーは，実
際のエコーをキャンセルするために，エコー経路のイ
ンパルス応答を推定する必要がある．その推定に適応
図 10.3 エコーキャンセラー
アルゴリズムを使用する．
3) システム同定
未知システムの特性を，その入出力から推定することを同定(identification)という．図 10.4
では，未知システムの特性を，FIR システムで近似する．そのために，実際のシステムの出
力 d(n)と近似システムのそれ y(n)との誤差 e(n)を最小にするように，
FIR のパラメータ h(n)
(n=0,1,...,M-1)を適応的に決める．この際に，適応アルゴリズムを使用する．
図 10.4 未知システムの同定
図 10.5 FIR 型適応フィルタ
10.3 適応フィルタ
1) FIR 型適応フィルタ
次式の FIR 型適応フィルタについて説明する(図 10.5 参照)．
M 1
y ( n)
h( k ) x ( n k )
(10.2)
k 0
理想応答 d(n)と適応フィルタの応答 y(n)との差 e(n) d (n) y (n) の２乗平均値
MSE E[e 2 (n)]
(10.3)
を最小にするよう，フィルタ係数 h(n) (n=0,1,...,M-1)を決める．上式に，(10.2)を代入し
て整理すると，次式となる．
信号処理 8 章-10 章.docx
- 79 -
大西
昇（名古屋大学）
MSE E[e 2 (n)] E[d 2 (n)] 2 E[d (n) y (n)] E[ y 2 (n)]
M 1
M 1 M 1
E[d 2 (n)] 2 h(k ) E[d (n) x(n k )] h(k )h(m) E[ x(n k ) x(n m)] (10.4)
k 0
k 0 m 0
M 1
M 1 M 1
k 0
k 0 m 0
Pd 2 h(k )rdx (k ) h(k )h(m)rxx (k m)
ここで，
Pd E[d 2 (n)]
rdx (k ) E[d (n) x(n k )]
(10.5)
rxx (k m) E[ x(n k ) x(n m)]
，rdx(k)は相互相関関数，rxx(k)は自己相関関数
で，Pd は理想応答 d(n)の 2 乗平均値（電力）
である．
2) 最適フィルタ係数
２乗平均値
M 1
M 1 M 1
k 0
k 0 m 0
MSE Pd 2 h(k )rdx (k ) h(k )h(m)rxx (k m)
(10.6)
は，h(k)の 2 次式であるので，次式を解いて h(k) (k=0,1,...,M-1)を求める．
M 1
MSE
2rdx ( k ) 2 h(m) rxx (k m) 0
h( k )
m 0
k 0,1, , M 1
(10.7)
これを書き下すと
rxx (1)
rxx ( M 1) h(0) rdx (0)
rxx (0)
r (1)
rxx (0)
rxx ( M 2) h(1) rdx (1)
xx
rxx (0) h( M 1) rdx ( M 1)
rxx ( M 1) rxx ( M 2)
(10.8)
となる．これは正規方程式と呼ばれる（線形予測係数についての連立方程式（9.10）との違
いは右辺である）． 式(10.8)の行列，ベクトルを
rxx (1)
rxx ( M 1)
rxx (0)
h(0)
r (1)
h(1)
(
0
)
(
2
)
r
r
M
xx
xx
, h
,
R xx
(
1
)
(
2
)
(
0
)
r
M
r
M
r
h
(
M
1
)
xx
xx
xx
rdx (0)
r (1)
p dx
r
(
M
1
)
dx
(10.9)
のように表すと，最適なパラメータ hopt は次式となる．
hopt R 1 p
(10.10)
つぎに，最適なパラメータ hopt での２乗平均値 MSEmin を求める．式(10.6)に hopt を代入し，
信号処理 8 章-10 章.docx
- 80 -
大西
昇（名古屋大学）
展開すると
M 1
M 1
M 1
k 0
m0
MSEmin Pd 2 hopt (k ) rdx ( k ) hopt (k ) hopt ( m) rxx (k m)
k 0
(10.11)
T
p
Pd 2 hopt ( k )rdx ( k ) hopt ( k )rdx ( k ) Pd hopt ( k )rdx (k ) Pd hopt
M 1
M 1
M 1
k 0
k 0
k 0
M 1
となる．上式の導出において，式(10.7)から得られる関係
h
m 0
opt
(m)rxx (k m) rdx (k ) を用
いている．
3) 行列を用いた式の導出
式(10.9)の行列とベクトルを用いると，式(10.6)の MSE は,
MSE Pd 2h T p h T Rh
となる．MSE の h による微分(付録３参照)は，R が対称行列であるので，
(10.12)
MSE
2 p 2 Rh
h
(10.13)
となる．したがって，最適な係数 hopt は
hopt R 1 p
となる．これを式(10.12)に代入すると，
T T
T T
MSEmin Pd 2hopt
p hopt Rhopt Pd 2hopt
p hopt RR 1 p
T T
T
T
p hopt p Pd hopt
p Pd hopt
p
Pd 2hopt
となり，これは式(10.11)と同じである．
10.4 適応アルゴリズム
FIR 型適応フィルタの最適係数 hopt は，自己相関行列 R と相互相関ベクトル p が既知であれ
ば，式(10.10)で求まる．しかし，逆行列の計算（O(M3)）が必要で，実時間処理に適さない．
また，自己相関行列 R と相互相関ベクトル p が既知でない場合もある．そこで，繰り返し計
算により， hopt を求める方法が必要となる．その方法として，最急降下法，最小２乗平均ア
ルゴリズムおよび逐次最小２乗法を以下で説明する．(詳細は，文献(25)を参照)
a) 最急降下法
最急降下法(steepest descent method)では，自己相関行列 R と相互相関ベクトル p が既知であ
るとする．そして，係数 hk を次式で更新する．
hk 1 hk k
, where k
MSE
h
h hk
2 Rhk 2 p
(10.14)
信号処理 8 章-10 章.docx
- 81 -
大西
昇（名古屋大学）
上式で，下付の文字は反復回数を表し，右辺第２項のμ（>0）はステップサイズである．
最急降下法が収束するために，μが満足すべき条件を求める．式(10.14)の両辺から，最
適値 hopt をひくと，次式となる．
hk 1 hopt hk hopt k hk hopt (2 Rhk 2 p )
p Rhopt
hk hopt 2 ( Rhk Rhopt )
(10.15)
上式で， k hk hopt とすると
k 1 k 2R k ( I 2R) k
(10.16)
となる．行列 R の固有値分解（付録４参照）を R DD とし，上式の両辺に，直交行列 DT を
T
左から掛ける．
D T k 1 D T k 2D T R k D T k 2D T DD T k D T k 2D T k
ここで， v k D T k とすると，上式は
v k 1 v k 2 v k ( I 2 )v k
(10.17)
(10.18)
上式の漸化式を解くと，次式となる．
v k ( I 2) k v 0
ベクトル v k の i 番目の要素を， v k (i ) とすると，上式は
(10.19)
v k (i) (1 2i ) k v 0 (i)
(10.20)
となる．これより， v k (i ) がゼロに収束するための条件は次式となる．
1 2 i 1 すなわち 0
1
(10.21)
i
したがって， v k の全ての要素がゼロに収束するための条件は
0
1
(10.22)
max
となる．ここで， max は行列 R の最大固有値である．なお， max が不明な場合には，μを十
分小さくすれば，最適解 hopt に収束する．
つぎに，ステップサイズμを固定でなく，反復に応じて変化させることを考える．その場
合，k 回目のステップゲインμk を，行列 R のｋ番目の固有値λk の逆数の半分，すなわち
k
1
2 k
(10.23)
とすれば，M 回の繰り返しで，最適解 hopt に収束する．（詳細は，文献(26)を参照）
b) 最小２乗平均アルゴリズム
実際には，自己相関行列 R と相互相関ベクトル p を正確に知ることはできず，勾配
k 2 Rhk 2 p が求まらないことが多い．そこで，最小２乗平均 (least mean squares, LMS)アル
信号処理 8 章-10 章.docx
- 82 -
大西
昇（名古屋大学）
ゴリズムでは，時刻 k での誤差の瞬時値 e 2 (k ) を使用して，次式のように勾配を求める．
T
2
d
k
h
x
(
)
e
k
(
)
ˆ
k k
k
hk
hk
2
(10.24)
T
こ こ で ， h k は 時 刻 k で の パ ラ メ ー タ ベ ク ト ル で hk (h(0), h(1), , h( M 1)) ，
x k ( x ( k ), x ( k 1), , x ( k M 1)) T である．この勾配を用いて，次式のように hk を更新す
る．
ˆ
hk 1 hk
k
(10.25)
ここで，μはステップサイズである．式(10.24)の勾配を計算すると，次式となる．
T
2
d
k
h
x )
(
(
)
e
k
e
k
(
)
(
)
ˆ
k k 2e(k ) x k
2e(k ) 2e(k )
k
hk
hk
hk
(10.26)
したがって，式(10.25)は
hk 1 hk 2 e(k ) x k
(10.27)
となり，これが最小２乗平均(LMS)の更新アルゴリズムである．LMS は時刻ｋにおける入力
xk ( x(k ), x(k 1), , x(k M 1))T と誤差 e(k)のみを使用するため，計算量は O(M)（M は次数）
と少ない．したがって，実時間処理に適したアルゴリズムである．
LMS アルゴリズムでの係数 hk が
lim
E[hk ] hopt
(10.28)
k
となるためのステップサイズμの条件は，最急降下法のそれ（式(10.22)）と同じである．
c) 逐次最小２乗アルゴリズム
自己相関行列 R と相互相関ベクトル p を，その時点までに入手した入力と理想応答から，
逐次推定しながら係数 hk を求める逐次最小２乗アルゴリズム(recursive least squares, RLS)を
説明する．
xi ( x(i ), x (i 1), , x(i M 1)) T として，時刻 k－1 までに得た入力信号行列 Xk-1, 理
想応答ベクトル d k 1 を次のように表す．
X k 1
x0T
T
x
1 ,
xT
k 1
d k 1
d (0)
d (1)
d (k 1)
(10.29)
時刻ゼロから k－1 までの，誤差の平均は
信号処理 8 章-10 章.docx
- 83 -
大西
q d (i) h x
昇（名古屋大学）
T
1
(10.30)
d k 1 X k 1 h Qk 1 d k 1 X k 1 h
k
i0
となる．ここで，qi は時刻 i での誤差の重み係数で， Qk 1 は q 0 , q1 , , q k 1 を対角要素にもつ
E
1
k
T
k 1
i
2
i
対角行列，すなわち Qk 1 diag (q 0 , q1 , , q k 1 ) である．上式を最小にする係数 hk 1 は
E
0 の解である．
h
T
T T
T
T T
E 1 d k 1Q k 1 d k 1 h X k 1Q k 1 d k 1 d k 1Q k 1 X k 1 h h X k 1Q k 1 X k 1 h
h k
h
1
2 X kT1Q k 1 d k 1 2 X kT1Q k 1 X k 1 h 0
k
(10.31)
であるので，時刻 k－1 までの誤差の平均を最小にする hk 1 は次式となる．
hk 1 X kT1Qk 1 X k 1
1
X kT1Qk 1 d k 1
(10.32)
式(10.32)の直接解法には，逆行列の計算（O(M3)）が必要で，計算量が大である．そこで，
係数 hk を逐次的に求めることを考える．
入力 x(k)と理想応答 d(k)を新たに入手できる．
すなわち，
時刻 k-1 から k と一時刻進むと，
x0T
X k T ,
xk 1
xT
k
d (0)
dk
d (k 1)
d (k )
(10.33)
時刻 k での最適な係数 hk は，式(10.32)で k-1 を k で置き換えると，
hk X kT Qk X k
1
X kT Qk d k
(10.34)
となる．ここで， Q k diag (q 0 , q1 , , q k 1 , q k ) である．
Pk X kT Qk X k
1
(10.35)
として，これに
X
Q
X k kT1 , d k d k 1 , Q k k 1
x
d
(
k
)
0
k
0
q k
(10.36)
を代入すると，次式となる．
0 X k 1
Q
Pk ( X kT 1 x k ) k 1
q k x kT
0
1
Pk11 q k x k x kT
1
X kT 1Qk 1 X k 1 q k x k x kT
1
(10.37)
上式に，次式の逆行列の補題（詳細は付録５）を適用する．
信号処理 8 章-10 章.docx
- 84 -
大西
T 1
A 1b c T A 1
1
( A bc ) A
1 c T A 1b
昇（名古屋大学）
(10.38)
式(10.38)で， A Pk11 , b c q k x k とおくと，式(10.37)は次式のように， Pk の漸化式とな
る．
q k Pk 1 x k x kT Pk 1
T
Pk Pk 1
T
( I g k x k ) Pk 1
1 q k x k Pk 1 x k
ここで，M 次元ゲインベクトル g k を
q k Pk 1 x k
gk
1 q k x kT Pk 1 x k
としている．
(10.39)
(10.40)
つぎに，係数 hk の更新式を求める．式(10.34)に式(10.39)を代入し，展開すると
hk X kT Q k X k
1
X kT Q k d k ( I g k x kT ) Pk 1 X kT1
( I g k x kT ) Pk 1 ( X kT1Q k 1 d k 1 q k x k d (k ))
( I g k x kT )hk 1 q k ( I g k x kT ) Pk 1 x k d (k )
Q
x k k 1
0
0 d k 1
q k d (k )
(10.41)
となる．式(10.40)より
q k ( I g k x kT ) Pk 1 x k g k
(10.42)
であるので，この関係を，式(10.41)の右辺第２項に適用すると，次式の更新式が得られる．
hk hk 1 g k x kT hk 1 g k d (k ) hk 1 g k (d (k ) x kT hk 1 ) hk 1 g k e(k ) (10.43)
T
ここで， e(k ) d (k ) x k hk 1 である．
以上をまとめると，RLS のアルゴリズムはつぎのようになる．
Step 1
（初期化）k=-1, h1 0, P1 I . εは小さな正の数．
Step 2
k=k+1.
Step 3
q k Pk 1 x k
1 q k x kT Pk 1 x k
hk hk 1 g k e(k )
gk
,
e(k ) d (k ) x kT hk 1
,
Pk Pk 1 g k ( x kT Pk 1 )
go to Step 2.
Step 3 の Pk の計算において，括弧内を先に計算すると，計算量は O(M2)になることに注意
（ g k x kT を先に計算すると，行列の乗算が必要となり，計算量は O(M３)）
．
信号処理 8 章-10 章.docx
- 85 -
大西
昇（名古屋大学）
11. ケプストラム解析
畳み込みの形で結合した２つの信号を分離する方法であるケプストラム解析を説明し，音声認
識やエコー検出の応用例を示す．
11.1 ケプストラム解析とは
2 つの信号のスペクトルに重なりがなければ，フィルタにより 2 つの信号を分離できる．しか
し，フィルタでは，畳み込みの形で結合した２つの信号を分離することはできない．
ケプストラム解析は，2 つの信号が畳み込みの形で結合した信号から，元の 2 つの信号を分離
し，個々の信号の解析を可能とする信号処理技術の一方法である．応用例には，音声解析，エコ
ー解析，地震波解析などがある．
11.2 ケプストラム
信号 s(t)が，インパルス応答が h(t)の伝達系を経て観察されるとする．すなわち，
x (t ) s (t ) * h(t )
(11.1)
地震波の場合，s(t)は震源で，h(t)は地殻特性となり，音声では，s(t)は声帯で発生した音で，
h(t)は声道という音声の生成系となる．
図 11.1 に，ケプストラム解析の概要を示す．式(11.1)をフーリエ変換すると，
X ( ) S ( ) H ( )
(11.2)
と，積の形になる．つぎに， X ( ) のパワーを求め，その対数をとると，次式のように，和の形
になる．
2
2
log X ( ) log S ( ) log H ( )
2
(11.3)
これを逆フーリエ変換すると，パワーケプストラム(cepstrum, spectrum からの造語)が得られる．
F 1 [log X ( ) ] F 1 [log S ( ) ] F 1 [log H ( ) ]
2
2
2
(11.4)
ここで， F 1 [ P ] は P の逆フーリエ変換を表す．パワーケプストラムは周波数領域からの逆変換
であるので，パワーケプストラムは時間領域での表現であるが，元の時間と区別するために，ケ
フレンシー(quefrency, frequency からの造語) 領域での表現とする．
X ( )
=s(t)*h(t)
S ( ) H ( )
2
2
2
log S ( ) log H ( )
2
2
フーリエ
対数
フーリエ
変換
変換
逆変換
パワー
ケプストラム
log X ( )
2
x(t)
ケフレンシー
リフター
図 11.1 ケプストラム解析の概要
式(11.4)により，信号 x(t)はケフレンシー領域に変換され，s(t)のパワーケプストラム
F 1 [log S ( ) ] と h(t) のパワーケプストラム F 1 [log H ( ) ] の和となる．
（x(t)=s(t)+n(t)と和の
2
2
講義信号処理 11 章付録文献.docx
- 86 -
大西
昇（名古屋大学）
形の場合，
周波数領域で S(ω)と N(ω)が分離できることがあったように，）
音声や地震波の場合，
これらのパワーケプストラムは，ケフレンシー領域で分離できる．ケフレンシー領域での分離に
はリフター(lifter, filter からの造語)を用い，これは周波数領域でのフィルタに対応する．リフタ
ーには，short-pass，long-pass, notch リフターがあり，それぞれ low-pass, high-pass, notch フィル
タに対応する．
11.3 音声のケプストラム解析
s(t)
x(t)
有声音は，図 11.2 のように，声帯で発生したイ
声道
ンパルス列の振動（音）が，声道（舌,唇,顎,歯な
h(t)
どの調音器官からなる音響的な管）という音声生
T
成系を経た結果である．
図 11.2 有声音の生成
生成系への入力 s(t)が周期 T のインパルス列
s (t )
N 1
(t kT )
(11.5)
k 0
とすると，出力である音声 x(t)は
x (t ) s (t ) * h(t )
(11.6)
となる．両辺をフーリエ変換すると
sin( NT / 2) j ( N 21)T
X ( )
e
H ( )
sin(
T
/
2
)
(11.7)
上式のパワーを求め対数を求めると，次のようになる．
2
log X ( )
2
sin( NT / 2)
2
log
log H ( )
sin(T / 2)
(11.8)
右辺の第１項は，ω=2π/T の高調波からなる線スペクトルで，ωについて急激に変化する（高
いケフルレンシーを持つ）
．一方，第２項はωについて緩やかに変化する（低いケフルレンシー
に限定される）
．そこで，式(11.8)を逆フーリエ変換し，ケフレンシー領域で，low-pass リフタ
ー（境界は T）を使用することで，生成系に相当する成分（ log H ( ) ）を分離できる．リフタ
2
ーで分離した成分をフーリエ変換すると log H ( ) を得て，これを指数化することで H ( ) を
2
2
抽出できる．声帯での振動 s(t)は,ケプストラム上に周期 T で現れるので，そのケフレンシーか
らピッチを得ることができる．詳細は文献(26)を参照のこと．
適用例
図 9.4 の音声波形にケプストラム解析を施した．8192 点をハミング窓により切り出し，
FFT（8192 点）を施した結果，図 11.3 のパワーケプストラムが得られた．ピークの間隔からピッチ
が 116.4Hz となり，これは線形予測法で求めた値と一致する．図 11.3 の最初のピーク（379 点目）
よりも低ケフレンシー領域のケプストラムのデータに FFT を施し，指数化することで，スペクトルの
包絡 H ( ) を求めた．図 11.4 の実線が H ( ) で，振動の激しい点線が X ( )
2
2
2
である．なめら
講義信号処理 11 章付録文献.docx
- 87 -
大西
昇（名古屋大学）
かなスペクトル包絡が得られていることが分かる．そして，約 800Hz，1.1kHz, 1.9kHｚにフォルマ
ントが見られる．
図 11.3 音声「あ」のパワーケプストラム
ピークがはっきり見えるよう，0.00025 以下の領域の強いパワーを０とした．
図 11.4 パワースペクトル
講義信号処理 11 章付録文献.docx
- 88 -
大西
実線が H ( ) で，振動の激しい点線が X ( )
2
昇（名古屋大学）
2
11.4 ケプストラム解析によるエコー検出
地震波では，直接波 s(t)とその反射（エコー，ここでは簡単のため１つとする） as (t ) の和
を観察する．すなわち，
x (t ) s (t ) as (t )
(11.9)
これをフーリエ変換すると
X ( ) S ( ){1 ae j }
となる（ H ( ) 1 ae j ）． X ( ) のパワーを求め，対数を取ると，
log X ( )
2
log S ( ) 1 ae j
2
2
log S ()
log 1 ae j
2
(11.10)
2
(11.11)
log S ( ) log{1 a 2 2a cos( )}
2
となる．式(11.11)の第２項は，次式のようにフーリエ展開できる（ a 1 とする）
．
(1) m 1 a m cos(m )
log{1 a 2a cos( )}
m
m 1
2
(11.8)
すなわち，式(11.11)の第２項は周期性があり，
無限個の高調波成分を含む (図 11.3 参照)．一方，
0.5
直接波成分 s(t)が不規則な信号であれば，式
(11.11)の第 1 項に周期性は見られない．そこで，
100
式(11.11)を逆フーリエ変換すると，図 11.4 のよ
-0.5
うに，周期性はケプストラム上の鋭いピークとな
-1
200
300
400
500
る．鋭いピークの存在により，エコーの存在を知
ることができ，ピークの位置から遅れ時間τ，ピ
ークの値から減衰係数 a を推定できる．
パワーケプストラム
F 1 log S ( )
τ
2
図 11.3 式(11.8)のグラフ(a=0.5,
τ=0.02) (横軸は周波数)
2
F 1 log 1 ae j
2τ
3τ
ケフレンシー
図 11.4 反射を含む信号のパワーケプストラム
講義信号処理 11 章付録文献.docx
- 89 -
大西
付録１
昇（名古屋大学）
窓関数（長さ N=2M+1）
方形窓(rectangular)
1
wr ( n )
0
バートレット窓(Bartlett)
n
1
wbr (n)
M
0
n M
n M
n M
n M
n
1 1
cos
whn (n) 2 2
M
0
ハニング窓(Hanning)
n M
n M
n
0.54 0.46 cos
whm (n)
M
0
ハミング窓(Hamming)
n M
2n
n
0.08 cos
0.42 0.5 cos
wbl (n)
M
M
0
ブラックマン窓(Blackmann)
I 1 n M
wk ( n ) 0
0
カイザー窓(Kaiser)
n M
2
n M
n M
n M
( x / 2) m
I 0 ( x) 1
m!
m 1
0 次第１種変形ベッセル関数
I 0 (0)
n M
2
付録２ ディジタルバタワース低域フィルタ
アナログバタワースフィルタ（カットオフ角周波数は１rad/s）の２乗振幅特性は次式である．
H a ( )
1
2
1
2N
この式から，伝達関数 H(S)の式が次のようになる．
H ( s) H ( s)
双一次変換 s k
1
1 s 2
N
1 z 1
を，上式に施すと
1 z 1
講義信号処理 11 章付録文献.docx
- 90 -
大西
H ( z ) H ( z 1 )
昇（名古屋大学）
1
1 z 1 2
1 k
1
1 z
N
z e j を上式に代入して，ディジタルフィルタの周波数特性を求めると
1
2
H ()
2
1 kj tan
2
N
上式のｋは，アナログフィルタのカットオフ角周波数は１rad/s が，ディジタルフィルタのそれ
Ωc に対応するので
k
1
tan c
2
である．したがって，ディジタルバタワース低域フィルタの２乗振幅特性は次式となる．
1
2
H ()
tan
2
1
tan c
2
2N
付録３ ベクトルでの微分
a , b をｎ次元列ベクトル， f (a ) はスカラー値をとる関数，M を正方行列(n x n)とする．
ベクトルでの微分は
f ( a ) f (a )
a
a1
T
f (a )
f (a )
a 2
a n
(A.1)
である．そして，次の等式が成立する．
(a T b ) (a T b )
a,
b,
a
b
(a T Mb )
(a T Mb )
Mb ,
MTa
a
b
(A.2)
(A.3)
(a T Ma )
( M M T )a
a
付録４
(A.4)
固有値分解（スペクトル分解）
正方行列 A (mxm)の固有値と固有ベクトルを j , e j ( j 1,2, , m) とする．そして，m 個の固有ベ
講義信号処理 11 章付録文献.docx
- 91 -
大西
昇（名古屋大学）
クトルが一次独立（行列が単純）である場合は，行列 A は
A e1
T
1
e1
E E T
em
m emT
(A.5)
となり，これを行列 A の固有値分解（スペクトル分解）という．
付録５
逆行列の補題
(nxn)行列 A, (nxm)行列 B, (nxm)行列 C に対して，A, A+BCT, I+CTA-1B は全て正則とすると，
( A BC T ) 1 A 1 A 1 B( I C T A 1 B) 1 C T A 1
特に，m=1 の時，B，C は n 次元ベクトル b , c となり，
T 1
A 1b c T A 1
1
( A bc ) A
1 c T A 1b
(A.6)
(A.7)
講義信号処理 11 章付録文献.docx
- 92 -
大西
昇（名古屋大学）
参考文献
2章
(1)佐々木公男，ディジタル信号処理
(2)宮川 洋，佐藤拓宋，茅
基礎理論と方法論，丸善，3.3.1，2001 年．
陽一，不規則信号論と動特性推定,コロナ社，p.10，1973 年．
(3)赤尾保男，堀井憲爾，過渡現象論，廣川書店，pp.48-51,1986.
3章
(4)小倉久直，確率過程入門，森北出版，1998．
(5)添田喬，中溝高好，大松繁，信号処理の基礎と応用，日新出版，1979．
(6)森下
厳，小畑秀文，信号処理，計測自動制御学会，3 章，1992 年．
4章
(7)佐々木公男，ディジタル信号処理
基礎理論と方法論，丸善，4.2，2001 年．
(8)辻井重男，鎌田一雄，ディジタル信号処理，昭晃堂，2.4.3，1996 年．
(9)武部，高橋，西川：ディジタル信号処理，６章，丸善，2004
(10)貴家仁志，マルチレート信号処理，昭晃堂，1999．
(11)相良岩男，AD/DA 変換回路入門，日刊工業新聞社，1999．
5章
(12)森下
厳，小畑秀文，信号処理，計測自動制御学会，p.28, 1992 年
(13)辻井重男，鎌田一雄，ディジタル信号処理，昭晃堂，pp.58-60, 1996 年．
6章
(14)貴家仁志，ディジタル信号処理，昭晃堂，3.4，1997．
7章
(15)佐々木公男，ディジタル信号処理
基礎理論と方法論，丸善，pp.120-123, 2001 年．
(16)樋口龍雄，川又政征，MATLAB 対応
ディジタル信号処理，昭晃堂，11 章，2000．
(17)酒井英昭
(18)森下
編著：信号処理，オーム社，5.1，1998．
厳：わかりやすいディジタル信号処理，昭晃堂，7.3，1996 年．
(19)N. K. Bose, Digital Filters –Theory and Applications, Elsevier Science Pub., 3.5, 1985.
8章
(20)佐々木公男，ディジタル信号処理
(21)森下
基礎理論と方法論，丸善，5.2.1，2001 年．
厳･小畑秀文，信号処理，計測自動制御学会，5.1, 1992 年
(22)磯部 孝編：相関関数およびスペクトル，東京大学出版会，pp.28-46,1975
講義信号処理 11 章付録文献.docx
- 93 -
大西
昇（名古屋大学）
9章
(23)森下
厳,わかりやすいディジタル信号処理,8 章,昭晃堂
(24)森下
厳,小畑秀文，信号処理，計測自動制御学会，6 章, 1992 年
10 章
(25)飯國洋二，適応信号処理アルゴリズム，培風館，2000．
(26)酒井英昭
編著：信号処理，オーム社， 6 章，1998．
11 章
(27)斎藤収三，中田和男：音声情報処理の基礎，オーム社，7.3，1981．
講義信号処理 11 章付録文献.docx
- 94 -
＜A-Z＞
自動波形等価器
79
バタフライ計算
65
AIC 基準
75
時不変システム
41
バタワースフィルタ
53
Levinson-Durbin の方法
74
遮断周波数
51
ハミング窓関数
56
集合平均
14
非再帰型システム
44
周波数サンプリング法
56
ビットリバーサル
65
循環推移定理
32
標準ｚ変換法
57
32
標本化
22
標本化周波数
22
＜あ＞
アップサンプリング
アナログ信号
26
2
アパーチャ効果
27
循環畳み込み定理
安定なシステム
42
信号
位相特性
46
信号確率集合
14
標本化定理
22
因果的なシステム
42
振幅特性
46
フィルタ
51
因果的な信号
29
ステップ関数
9
不規則信号
14
インパルス応答
41
スペクトル
6
複素フーリエ級数展開
5
ウィナー・ヒンチンの定理
20
正規化角周波数
46
フーリエ級数展開
4
エイリアシング
24
z 変換
33
フーリエ級数法
エコーキャンセラー
79
線形システム
41
フーリエ変換
エルゴード性
19
線形予測法
69
偏自己相関係数
全極モデル
69
全零モデル
69
窓関数
67
モーメント関数
15
＜か＞
ガウス性不規則信号
18
1
確定信号
2
双一次変換法
58
確率過程
14
相互相関関数
20
完全直交系
4
＜た＞
ギブス現象
6
ダウンサンプリング
基本角周波数
4
畳み込み積分
26
54
6
71
＜ま＞
＜や＞
予測係数
70
予測残差
70
7
＜ら＞
キュムラント関数
17
チェビシェフフィルタ
54
ラプラス変換
10
極
34
直線位相
52
離散時間システム
40
極零モデル
69
定係数差分方程式
43
離散時間フーリエ変換
29
群遅延
46
定常
19
離散フーリエ変換
31
ケフレンシー
86
ディジタル信号
リフター
87
ケプストラム解析
86
適応アルゴリズム
78
留数
13
高速フーリエ変換
63
適応信号処理
78
量子化
24
零点
34
＜さ＞
2
デルタ関数
8
再帰型システム
44
同定
79
最急降下法
81
特性関数
15
最小 2 乗平均
82
時間平均
15
ナイキスト周波数
22
自己相関関数
19
ノッチフィルタ
52
時系列
69
システム関数
44
＜な＞
＜は＞
パーセバルの等式
8
95
講義信号処理索引.doc
© Copyright 2024 ExpyDoc
