Design Signal Processing

◎信号処理
◆デジタル信号処理とサンプリング
※デジタル信号処理の基本
「入力信号がどのような信号の経路(信号パス)を通って出力されるか」
_◇デジタル信号処理回路への入力
数値系列
有限のビット数の時系列データ
振幅、時間ともに不連続なデータ量
※時間に対する不連続性
sampling
標本化
※振幅に対する不連続性
quantization
量子化
①アナログ信号x(t)をある時間ごとに測定したデータ系列⇒離散信号(時間についてのみ不連続)
⇒通常、サンプリング間隔Tを一定とする
⇒Tを明示せず、t = k*T となるようなサンプル点kをもってx[k]と表す
⇒サンプリング間隔 T[s]の逆数をサンプリング周波数という。
fs = 1/T [Hz]
⇒サンプリング周波数は1秒間に府j組まれる離散信号の総個数
②離散信号に対して量子化をおこなって有限のビット長のデジタル値であらわす
⇒ xq[k]
_◇シャノンのサンプリング定理
標本化定理
※サンプリング周波数Fs
⇒Fs/2以上の信号成分を低域通過フィルタ(アンチエリアシング・フィルタ)を通して除去した後、離散的な値に変換しても元の情報は失われない。
⇒サンプリング周波数をアナログ情報の最大周波数の2倍以上に設定すれば、その情報からアナログ情報を生成できる。
T≦1/(2*f0)
fs≧2*f0
T:サンプリング間隔
f0:信号に含まれる最高周波数f0[Hz]
fs:サンプリング周波数[Hz]
※ナイキスト周波数 (Nyquist)
サンプリング周波数の半分の周波数
CD: 可聴範囲上限20kHzの2倍の40kHzに1割ほどの余裕をもたせ、44.1kHz
_◇エイリアシング
aliasing
異なる周波数のアナログ信号からサンプリングして得られるデジタル信号がまったく同じ値になってしまう現象
※雑音になるので「折り返し雑音」と呼ぶ
※ナイキスト周波数以上の信号をサンプリングしてデジタル化すると、エイリアシングのために違う周波数の信号が見えてしまう
⇒折り返した周波数に現れる
※アンチ・エイリアシング・フィルタ
折り返し雑音となる高い周波数をカットして、サンプリング定理を満たすようにするローパス・フィルタ
※サンプル&ホールド(S&H)回路
サンプリングすると同時に信号値を保持して離散信号を作り出す。
_◇量子化処理
※量子化幅⊿Q (量子化ステップサイズ)
⊿Q=(表現可能な最大幅)/2^L
=((信号の最大幅)-(信号の最小値))/2^L
L:信号を量子化するビット数
※量子化雑音
量子化によって生じる丸め誤差 ±LSB/2
_◇デルタ関数
(ディラックの)
インパルス。パルスの幅が限りなくゼロに近く、かつパルスの面積は1になる、という関数。
δ(t)=0 … t≠0のとき
∫[-∞:∞]δ(t)dt = 1
※正規分布の密度関数から考える
f(t) = 1/(√(2π)*σ)*exp[-(t-μ)^2/2*σ^2]
δ(t) = lim[σ->0] f(t)
_◇サンプリングの数学的表現
サンプリングされた信号をX(t)とするとき
X(t)=∑[n=-∞:∞]X(t)*T*δ(t-n*T)
 =∑[n=-∞:∞]X(nT)*T*δ(t-n*T)
X(t)は連続関数なので、0区間も意味を持つ
_◇ゼロ次ホールド効果
サンプリング定理が正確に成り立つのは、純粋なδ関数を使って出力した場合であるので「階段状」にある幅をもって出力したものをスムースフィルタにとおしても歪みが発生する。(高周波特性が劣化する。)
|sin(π*f/fs)|/(π*f/fs)
fs/2において-3.92dB
幅Tの矩形波の代わりにデューティを少なくして、その分、振幅を増やし、残り区間にゼロを挿入することでサンプリング定理の定義に近づけると歪特性が改善される。
デューティ50%にすることで-0.91dBで収まる。
_◇オーバサンプリング
Over Sampling
デューティを変えることは難しいことが多いので、サンプリング周波数に対して信号帯域を相対的に低くして、歪を目立たなくする。オーバサンプリングによりゼロ次ホールド効果が対策されるだけでなく、スムージングフィルタの特性も簡単化される。
①カットオフf/2のローパスフィルタに通された信号を周波数fでサンプルしても2fでサンプルしてもサンプリング定理によれば信号は等価(ただし、f/2以上にあらわれるエリアシング部分が異なる)
②サンプルされた信号はインパルスとゼロの信号なのでインパルス間にゼロを入れても歪は発生しない。(信号を歪ませないためにはゼロでなければならない)
③そこでオーバサンプリングをするためにゼロを挿入。(例えば2倍であれば中間点にゼロを挿入)
④オーバサンプリングフィルタでエリアシングスペクトルを取り除けば、オーバーサンプルをしたことになる。
_◇2次元光学画像におけるサンプリング
※CCDで光学的画像をとりこむ場合
⇒フォトダイオードでサンプリング
⇒画素の間隔で決まる空間サンプリング周波数の半分以下の解像度の映像成分のみがCCD表面で結像する必要がある
⇒そうでないとモアレ(一種のエリアシング)が発生
⇒モアレ対策:光学的ローパスフィルタ(水晶フィルタ)の挿入
◆実信号の複素数的解釈
現実世界の実信号は、2つの複素共役信号の合成されたものと考える。(複素数の虚数部分は打ち消しあって観測できるのは、実数部分だけになる)
実信号X
X = (Z + Zc) / 2
ここで
Z = X + jY, Zc = X – jY
例)振幅1の正弦波Xは
X = (e^(i*2π*Ti/n)+e^(-i*2π*Ti/n))/2
複素平面状の半径1/2の円周上を互いに反対方向に回転するベクトルの合成となる。
※AD変換されたデジタル信号を、区間Tで複素フーリエ変換すると、かならず共役複素数ペアの結果が現れる。
※デジタル信号処理の内部では、信号を抽象的な複素信号として扱い、最後にDA変換時に実数に戻す。
※複素信号で信号処理することにより、複素平面上の点で符号化(コンスタレーション)するQPSK,QAMなどが可能になる。
※複素信号は実際に観測される信号ではないので、解析信号と呼ばれる。
※アナログミキサー(周波数変換)
実信号の信号では周波数の符号は変えられず、周波数ゼロ以下に変換されるべき信号はプラス方向に折り返し歪となる。
※複素ミキサー
マイナスの周波数を持つものとして広がるので歪は発生しない。
例)スーパヘテロダイン受信機とダイレクトコンバージョン
①アナログ処理では、折り返し成分が無いような比較的高い周波数に中間周波数IFを選ぶ必要がある。
②DSPで解析信号を使えば、折り返しは心配ないので、IF周波数をゼロにできる。
→ダイレクトコンバージョン
アンテナ→帯域フィルタ1→ミキサ1<=第1局発
(ミキサ1)→帯域フィルタ2→第1IFアンプ→
ミキサ2<=第2局発
(ミキサ2)→帯域フィルタ3→第2IFアンプ→
→(10.7MHz, 455kHzなど)→検波→アンプ→出力
_◇信号の偶関数成分と奇関数成分への分解
※偶関数と奇関数は相関ゼロであり直交している
⇒1次元の信号を2次元平面上のベクトルの動きとして書ける。
※基本周期T、サンプリングはT/n秒
信号はTの中心にたいして、偶関数成分と奇関数成分の合成として表すことができる。
_◇複素表示
奇関数部分を虚数軸、偶関数部分を実数軸として表す。
※始まりの位相π/4の正弦波の場合
sin(π/4 + 2πi/n) … i=0~n-1
=sin(π/4)cos(2πi/n) + cos(π/4)sin(2πi/n)
=(1/√2)[cos(2πi/n) + sin(2πi/n)]
ベクトルの回転速度が信号の周波数に比例する。
_◇負の周波数
※同じ回転速度でも反時計方向と時計方向の2つの場合が考えられるので、これを周波数の符号と考える。
◆z変換
デジタル信号
離散時間 k (整数 -∞<k<∞)
数値 x[k]
サンプル時間 T[s]
kが非整数の時刻では、信号の値は定義されない(0ではない)
{x[k]}[k=-∞:∞]
⇒数値集合、順序集合
サンプル時間Tの時間遅れを変数 Z^-1であらわす
k*T秒のおくれは
Z^-k
※べき乗の形式で時間原点0からの遅れ時間を表す
※Z^0=1
※定数cは c*Z^0なので、時刻t=0における信号値が3を意味する。
※-Z^6は,時間原点より6*T秒左にある、信号値が-1であることを意味する
k*T[s]遅れた時間における振幅x[k]は
x[k]*z^-k
※このようにべき乗表現することで、デジタル信号全体を変数zのべき乗表現(べき級数和)であらわすことができる
X(z) = ∑[k=-∞:∞]{ x[k]*z^-k }
⇒両側z変換
X(z)はデジタル信号{x[k]}[k=-∞:∞]のz変換と呼ばれる。
⇒z変換はデジタル信号系列{x[k]}[k=-∞:∞]のすべての信号値を表す
※起点とする時刻をk=0として、そこからの信号を考えれば
X(z) = ∑[k=0:∞]{ x[k]*z^-k }
⇒片側z変換
※逆z変換
X(z)から{x[k]}[k=-∞:∞]を求める操作
※z変換対
x[k]⇔X(z)
_◇インパルス波形のz変換
δ[k]= 1 k=0
0k≠0
インパルス波形のz変換は
∑[k=0:∞]δ[k]*z^-k = 1
インパルス波形を n*T 秒遅らせると
z^-n
※インパルス応答
単位インパルス波形の入力に対するデジタル・システムの応答出力。その信号波形はフィルタ(広義にはデジタルシステム)の別表現となる。
_◇ステップ波形のz変換
U[k]= 1 k≧0
0k<0
※等比級数の公式
1+λ+λ^2+…+λ^K = ∑[k=0:K]λ^k
=(1 – λ^(K+1))/(1-λ)
K→∞のとき|λ|<1ならば
∑[k=0:∞]λ^k=1/(1-λ)
ところで
ステップ波形{U[k]}[k=0:∞]のz変換対は
U[k]⇔∑[k=0:∞]{U[k]*z^-k}
=1+z^-1+z^-2+…+z^-k+…
λ=z^-1とすれば、等比級数の公式をつかい
U[k]⇔ 1/(1-z^-1)
ただし、
|λ|=|z^-1|=|1/z|=1/|z|<1 より |z|>1のとき
※収束領域
級数和が収束しz変換を計算可能
上記の例では|z|≦1では級数和は発散し、z変換を計算不能
例)z変換
X(z)=1/(1-0.5*z^-1)
を有する x(t)
分母がない無限級数の総和の形に変形する
1/3=0.333…=3*0.1+3*(0.1)^2+…
=0.3*{1+(0.1)+(0.1)^2+…}
1/(1-0.5*z^-1)の割り算を実行すると
=1+0.5*z^-1+0.25*z^-2+0.125*z^-3+…
….._◇指数関数波形のz変換
_◇指数関数波形のz変換
x[k]=α^k(ただし、k<0のときx[k]=0, αは複素数も含む)
⇒ステップ波形{U[k]}[k=0:∞]を乗ずる形で
{α^k*U[k]}[k=0:∞] =α^k k≧0
=0k<0
再び等比級数公式で
|λ|=|α*z^-1|=|α/z|=|α|/|z|<1
よって|z|>|α|
の範囲でz変換は
1/(1-α*z^-1)
_◇極 (pole)と零点
X(z)=∞
となるzの値を極と呼ぶ
H(z)=A(z)/B(z)
と記すならば、H(z)の極はB(z)=0とおいた方程式の根である。
⇒極周波数におけるゲインは無限大
X(z)=0
となるzの値を零点とよぶ
H(z)=A(z)/B(z)
と記すならば、H(z)の極はA(z)=0とおいた方程式の根である。
零点周波数におけるゲインは0
例)
H(z)=(1+a*z^-1)/(1+b*z^-1)
の場合
B(z)=1+b*z^-1=0より, zp=-b
A(z)=1+a*z^-1=0より、zq=-a
※極はZ変換の分母を0とするZの値
例)
1/(1-α*z^-1)の場合、 Zp=α
※最終的に発散しないためには極Zpとα(複素数)は
k→∞で、|x[k]|→0より
|Zp|=|α|<1
⇒複素平面上の単位円の中
※極の位置とインパルス応答の波形
複素平面上の単位円の中=減衰波形
複素平面上の単位円上=持続振動
単位円の外側=増大波形
※零点
安定性には関係ない
_◇正弦波のz変換
x(t) = sin(ω*t)*U(t) = sin(ω*t) t≧0
0t<0
T秒でサンプリングして得られるデジタル信号
{x(k*T)}[k=0:∞]
のz変換
x(k*T)=sin(k*ω*T)*U[k]
ここでオイラーの公式により
e^jθ =cosθ+j*sinθ
e^-jθ=cosθ-j*sinθ
より
sinθ=(1/2j)*(e^jθ-e^-jθ)
ここでθ=k*ω*Tとおけば
sin(k*ω*T)=(1/2j)*(e^jk*ω*T-e^-jk*ω*T)
e^jk*ω*T=(e^j*ω*T)^k
のように考えれば、指数関数のz変換の式を適用でき、
{e^jk*ω*T}*U[k]は 1/(1-(e^j*ω*T)*Z^-1) のように変換できる
2項を通分して整理し、オイラー公式を使うことで
sin(k*ω*T)*U[k] は
      sin(ω*T)*Z^-1
————————–
1 – 2*cos(ω*T)*Z^-1 + z^-2
cos(k*ω*T)*U[k] は
      1 – cos(ω*T)*Z^-1
————————–
1 – 2*cos(ω*T)*Z^-1 + z^-2
とz変換できる
例)周波数5Hz, ω=10π[rad/s]
_◇時間遅れとz変換
{x[k]}[k=-∞:∞]
ただし、x[k]=0, k<0のとき
を1サンプル時間T秒おくらせたものを g[k] であらわすと
g[k] = x[k-1]
x[k]のZ変換がX(z)であれば
g[k] = z^-1 * X(z)
※z^-1は信号全体を1サンプル時間 T だけ遅らせる働きをする
⇒mサンプル遅らせるには z^-mを乗ずればよい
_◇z変換とブロック図
時間領域のブロック図はz変換領域に変換できる。
①遅延回路Tは、z^-1の乗算に置き換えられる
⇒遅延ブロックを通ると入力X(z)は z^-1 * X(z)に変化する
②係数乗算、信号加算、信号分岐は変化しない
③入力信号 x[k]はX(z)
④出力信号 y[k]はY(z)
_◇伝達関数と周波数特性
※アナログフィルタ
周波数によって値(複素インピーダンス)が変化するコイル、コンデンサにより周波数特性が得られる。
j*2πf*L
1/j*2πf*C
※デジタルフィルタの場合
例)
差分方程式 y[k]=a*x[k]+x[k-1]
⇒複素正弦波x(t)=e^(j*2πf*t)を代入したときの出力信号を求めれば周波数特性が分る
x[k]=x(k*T)=e^(j*2πf*k*T)
x[k-1]=e^(j*2πf*(k-1)*T)
これより
y[k]=e^(j*2πf*k*T)*(a+e^(-j*2πf*T))
y[k]/x[k]=a+e^(-j*2πf*T)
⇒ディジタルフィルタの周波数特性は、遅延回路による
※z変換
Y(z)=a*X(z)+z^-1*X(z)
H(z)=Y(z)/X(z)=a+z^-1
※z^-1 = e^(-j*2πf*T) =e^(-j*ω*T)
⇒z=e^(j*2πf*T) =e^(j*ω*T)
⇒伝達関数H(z)のzに e^(j*ω*T) を代入すれば周波数特性が得られる
_◇サンプリング定理とエイリアシング
※デジタルフィルタの周波数特性は、サンプリング周波数 fs ごとに繰り返す。 fs=1/T
nが整数のとき、e^(j*2π*n)=1なので
e^(j*2πf*T)=e^(j*2πf*T)*e^(j*2π*n)
=e^{(j*2π*T)*(f+n/T)}
=e^{(j*2π)*(f+n*fs)*T}
※正規化周波数Ω
Ω=f/fs
0~0.5の範囲の値をとる(サンプリング定理による)
※f=0Hzと f=fs/2 に相当する z の値
z=e^(j*2π*T*0)=e^0=1
z=e^(j*2π*T*fs/2)=e^(jπ)=cos(π)+j*sin(π)=-1
⇒z=1を伝達関数に代入すればDC, Z=-1を代入すれば最高周波数の、ゲインと位相の値を把握できる。
◆デジタルフィルタ:差分方程式による表現
入力┌────┐出力
  │デジタル│
─→┤信号処理├─→
  └────┘
処理=四則演算
※設計概念:
入力から希望出力が得られる不変の計算式を求める
(推定する、近似する)
⇒不変な計算式を得ることがデジタル・フィルタの設計そのもの
⇒不変な計算式は差分方程式と呼ばれる
※平均値
⇒信号の平滑化処理
y[k]=0.5*{x[k]+x[k-1]}
※差分
⇒エッジ検出
y[k]=0.5*{x[k]-x[k-1]}
※RCフィルタ相当
y[k]=0.5*y[k-1]+x[k]
※デジタル・フィルタとは、デジタル化された信号系列
{x[k]}[k=-∞:∞]
を別の信号系列
{y[k]}[k=-∞:∞]
に変換する機能を有する演算回路
⇒この数式表現が差分方程式
◆デジタルフィルタ:ブロック図による表現
※矢印(エッジ)は常に値を持つ。値は時間とともに変化する
※分岐(黒丸)は同じ値を2箇所に出力する
※ループには必ず遅延回路が含まれなければならない
_◇定数倍回路(乗算回路)
      │\a
x[k]─→┤ >─→y[k]=a*x[k]
      │/
_◇加算回路
2入力1出力加算回路
x1[k]─┐ ┌─┐
      └→│+├─→y[k]=x1[k]+x2[k]
      ┌→│ │
x2[k]─┘ └─┘
_◇遅延回路
(シフトレジスタ)
y[k]=x[k-1]
1サンプリング間隔T秒遅らせる処理
mサンプル遅らせる処理
y[k]=x[k-m]
※遅延回路をm個直列に接続する
時間領域での表現
     ┌─┐
x[k]→┤T├─→y[k]=x[k-1]
     └─┘
z変換領域での表現
     ┌────┐
X(z)→┤z^-1├─→Y(z)=z^-1*X(z)
     └────┘
_◇例
① y[k]=2*x[k]+3*x[k-1]-4*x[k-3]
x[k]┌─┐  ┌─┐┌─┐
 →●─┤T├●┤T├┤T├┐
  │ └─┘│└─┘└─┘││\-4
  │    │      └┤ >┐
  │    │   │\3 │/ │┌─┐
  │    └───┤ >┐   └┤+├→y[k]
  │        │/ │┌─┐┌┤ │
  │     │\2   └┤+├┘└─┘
  └─────┤ >────┤ │
        │/     └─┘
② y[k]=0.5*y[k-1]+2*x[k]
x[k]       ┌─┐
───┐   ┌───┤T├┐
   │   │   └─┘│
   ▽2  ▽0.5   │
   │  ┌┴┐     │
   └──┤+├─────●→y[k]
      └─┘
③ w[k]=4*x[k]-3*x[k-1], y[k]=2*w[k]+3*w[k-1]
x[k]┌─┐
───┬┤T├┐
   │└─┘│
   ▽4  ▽-3
   │  ┌┴┐w[k]┌─┐
   └──┤+├──●─┤T├┐
      └─┘  │ └─┘│
           ▽2   ▽3
           │   ┌┴┐
           └───┤+├→y[k]
               └─┘
※直列型構成
④ v[k]=x[k]-x[k-1], w[k]=2*x[k]+3*x[k-1], y[k]=v[k]+w[k]
        ┌─┐
     ┌─●┤T├┐
     │ │└─┘│
     │ ▽1  ▽-1
     │ │  ┌┴┐
x[k] │ └──┤+├─┐v[k]
─────●    └─┘ │
     │  ┌─┐   │
     └─●┤T├┐ ┌┴┐
       │└─┘│ │ ├→y[k]
       ▽2  ▽3└┬┘
       │  ┌┴┐ │
       └──┤+├─┘w[k]
          └─┘
※並列型構成
◆デジタルフィルタ:伝達関数
※伝達関数(システム関数)
     出力信号の全体 出力信号のz変換
H(z)=-------=--------
     入力信号の全体 入力信号のz変換
_◇インパルス応答
インパルス波形
δ[k] =1 k=0
=0k≠0
入力波形のz変換
X(z)=∑[k=0:∞]x[k]*z^-k=1
なので
Y(z)=H(z)*X(z)=H(z)
インパルス波形を入力すると伝達関数そのものの情報が得られる
_◇差分方程式から伝達関数を求める
例1)
y[k]=2*x[k]+3*x[k-1]-4*x[k-2]
z変換すると
Y(z)=2*X(z)+3*z^-1*X(z)-4*z^-2*X(z)
=(2+3*z^-1-4*z^-2)*X(z)
よって
H(z)=2+3*z^-1-4*z^-2
例2)
y[k]=0.5*y[k-1]+y[k-2]-x[k-1]
z変換すると
Y(z)=0.5*z^-1*Y(z)+z^-2*Y(z)-z^-1*X(z)
移項して整理すると
(1-0.5*z^-1-z^-2)*Y(z)=-z^-1*X(z)
よって
H(z)=Y(z)/X(z)=(-z^-1)/(1-0.5*z^-1-z^-2)
※z^-1に関する分母多項式の定数項はかならず1になるように正規化する。
※z^-1に関する多項式の比として分母、分子から構成される関数を有理関数とよぶ
例3)
w[k]=4*x[k]-3*x[k-1], y[k]=2*w[k]+3*w[k-1]
⇒2つの式をそれぞれz変換してW(z)を消去する
W(z)=4*X(z)-3*z^-1*X(z)
Y(z)=2*W(z)+3*z^-1*W(z)
=(2+3*z^-1)*(4-3*z^-1)*X(z)
H(z)=Y(z)/X(z)=(2+3*z^-1)*(4-3*z^-1)
例4)
v[k]=x[k]-x[k-1], w[k]=2*x[k]+3*x[k-1],
y[k]=5*v[k]-2*w[k]
それぞれz変換してV(z),W(z)を消去する
H(z)=9*z^-1
_◇ブロック図と伝達関数
①直列接続
前段
W(z)=H1(z)*X(z)
後段
Y(z)=H2(z)*W(z)
全体の伝達関数は
H(z)=Y(z)/X(z)=H2(z)*H1(z)
※関数の積の形=基本ブロックの直列接続
※前段と後段の基本ブロックを入れ換えても同一
②並列接続
2つの基本ブロックの出力を重み付けして足し合わせたもの
     ┌─────┐V(z)
    ┌┤H1(z)├┐
    │└─────┘│a
X(z)│      ┌┴┐
→───●      │+├→
    │      └┬┘
    │┌─────┐│b
    └┤H2(z)├┘
     └─────┘W(z)
H(z)=a*H1(z)+b*H2(z)
③フィードバック接続
X(z)┌─┐W┌─────┐
→───┤+├─┤H1(z)├●→Y(z)
    └┬┘ └─────┘↓
     ↑  ┌─────┐│
     └──┤H2(z)├┘
       V└─────┘
Y(z)=H1(z)*W(z)
W(z)=X(z)+V(z)
V(z)=H2(z)*Y(z)
WとVを消去すると
H(z)=Y(z)/X(z)=H1(z)/{1-H1(z)*H2(z)}
※再帰する信号経路の伝達利得=ループ利得
⇒伝達関数の分母を規定する
※非再帰の信号経路の伝達利得
⇒伝達関数の分子を規定する
フィードバック接続全体の伝達関数は
     非再帰経路の伝達利得
H(z)=-----------
     1-再帰経路の伝達利得
※分母の定数項は1になるように正規化する
※ブロック図から伝達関数を求める
入力から出力までのすべての信号経路を考え、経路に含まれる枝の積の総和が伝達関数になる
乗算回路を含まない信号枝は重み1
遅延回路は重みz^-1とする
加算回路は枝の結節点となるだけで積には入らない
例1)伝達関数を求める
x[k]
┬────┬─────┬─────┬─────┐
│    ▽a3   ▽a2   ▽a1   │
│┌─┐┌┴┐┌─┐┌┴┐┌─┐┌┴┐┌─┐┌┴┐y[k]
└┤T├┤+├┤T├┤+├┤T├┤+├┤T├┤+├→
 └─┘└─┘└─┘└─┘└─┘└─┘└─┘└─┘
H(z)=1+a1*z^-1+a2*z^-2+a3*z^-3+z^-4
例2)
x[k]  ┌─┐    y[k]
─→─┬>─┤+├──┬─→
   │a0└┬┘  ↓
   │  ┌┴┐  │
   │  │T│  │
   │  └┬┘  │
   │  ┌┴┐b1│
   └>─┤+├─<┤
    a1└┬┘  │
      ┌┴┐b2│
      │T├─<┘
      └─┘
※非再帰ブロックの伝達関数A(z)は
A(z)=a0+a1*z^-1
※再帰部分の伝達関数B(z)は
B(z)=b1*z^-1+b2*z^-2
全体の伝達関数は
H(z)=A(z)/(1-B(z))=(a0+a1*z^-1)/(1-b1*z^-1-b2*z^-2)
例3)
x[k] ┌─┐q[k]
─────┤+├─────┐
     └┬┘     │
   0.5△┌────┐│
      ├┤T   ├┤
     3▽└────┘▽4
      │ ┌─┐  │
      └─┤+├──┘
        └┬┘
         ↓y[k]
※非再帰ブロック
A(z)=4+3*z^-1
※再帰ブロック
B(z)=0.5*z^-1
H(z)=(4+3*z^-1)/(1-0.5*z^-1)
◆デジタルフィルタ:数式表現
※デジタルフィルタ
入力、出力1個ずつ
※入出力関係を表す特性パラメータ
伝達関数 H(z)=Y(z)/X(z)
周波数スペクトル
ゲイン特性A(f)
フェーズ特性θ(f)
極と零点
零点H(z)=0となるzの値
極H(z)→∞となるzの値
インパルス応答
伝達関数の逆z変換
_◇複素表示とフェーザ表示
※複素表示(直交座標)
w=a+j*b なる2次元量で扱う
※フェーザ表示(極座標)
A*e^(j*2πft)
⇒cos波とsin波を同時にあらわす(オイラーの公式)
A*e^(j*2πft)=A*cos(2πft)+j*A*sin(2πft)
⇒e^(j*2πft)は、単位円上の周波数f[Hz]で反時計方向に回転する周期関数
※伝達特性
周波数fの複素正弦波入力にたいする周波数特性を表す
H(f)=Y(f)/X(f)
⇒H(f)=Re{H(f)} + j*Im{H(f)}
⇒H(f)=|H(f)|*e^(j∠H(f))
|H(f)|=|Y(f)|/|X(f)|=出力の振幅/入力の振幅
。。。ゲイン
∠H(f)=∠Y(f)-∠X(f)
。。。位相
入力信号の周波数f[Hz]、時間のずれを t0[Sec]
∠H(f)=2πf*t0
位相が+で進み位相、負で遅れ位相
t0=位相/角周波数ω=2πf
入力 x(t)=e^(j*2πf*t)
出力 y(t)=A*e^(j*2πf*(t-t0))
y(t)/x(t)=A*e^(-j*2πf*t0)=A*e^(-jω*t0)=H(f)
◆デジタルフィルタの表現例
_◇伝達関数から他の表現を求める
H(z)= (4+3*z^-1+2*z^-2) / (1-0.5*z^-1-0.06*z^-2)
H(z)= Y(z)/X(z)
であるから
(1-0.5*z^-1-0.06*z^-2)*Y(z)=(4+3*z^-1+2*z^-2)*X(z)
ここで
X(z)をx[k], z^-1*X(z)をx[k-1]のように置き換えれば
差分方程式
y[k]=0.5*y[k-1]+0.06*y[k-2]+4*x[k]+3*x[k-1]+2*x[k-2]
を得る。
この差分方程式から以下のブロック図が得られる
x[k]┌─┐ ┌─┐
→──┬┤T├┬┤T├┐   ┌──→y[k]
   │└─┘│└─┘│   │ ┌─┐ ┌─┐
  4▽  3▽  2▽   ├─┤T├┬┤T├┐
  ┌┴┐ ┌┴┐ ┌┴┐ ┌┴┐└─┘│└─┘│
  │+├─┤+├─┤+├→┤+│0.5▽   ▽0.06
  └─┘ └─┘ └─┘ └┬┘  ┌┴┐  │
               └─←─┤+├──┘
                   └─┘
※信号処理プログラム
ブロック図の遅延素子Tを左からw1,w2,u1,u2とする
function[youtput] = FILTER(xinput)
w1=0; w2=0; u1=0; u2=0;
for k = 1:1:length(xinput)
xin = xinput(k);
y1=4*xin+3*w1+2*w2;
yout = y1+0.5*u1+0.06*u2;
w2=w1; w1=xin; u2=u1; u1=yout;
youtput(k) = yout;
end;
endfunction
◆フィルタリング
filtering
※filter
周波数選択性
必要な周波数成分を通過
不要な周波数成分を除去
「フィルタとは入力信号の周波数によってゲインや位相が変化し、出力信号の大きさと入力信号に対する時間のずれを同時にコントロールする回路である」
_◇フィルタの種類
①ローパスフィルタ LPF
遮断周波数 fc以下は通過、fc以上は除去
②ハイパスフィルタ HPF
遮断周波数 fc以下は除去、fc以上は通過
③バンドパスフィルタ BPF
fl < fhなる周波数に対し、flとfhの間の周波数は通過、それ以外を除去する
⇒中心周波数fc, fl, fhを遮断周波数、fh-flを(通過)帯域幅という
fc=√(fl*fh), W=fh-fl
④バンドエリミネートフィルタ BEF
fl < fhなる周波数に対し、flとfhの間の周波数を除去、それ以外を通過させる
⇒中心周波数fc, fl, fhを遮断周波数、fh-flを(阻止)帯域幅という
⇒バンドリジェクトフィルタ、バンドストップフィルタ、ノッチフィルタとも言う
fc=√(fl*fh), W=fh-fl
⑤オールパスフィルタ APF
すべての周波数を同じ大きさで通過させる。ゲインは一定だが位相だけが周波数によって変化する。
⇒位相ひずみの補正などに使われる
※通過域
pass band
※阻止域
stop band
あるいは、減衰域
attenuation band
※遮断周波数 あるいは カットオフ周波数
通過域と阻止域の教会
_◇フィルタのゲインと周波数特性
※通過/除去=ゲインの違い
入力と出力の振幅比の違い
利得
※実際のフィルタは、通過域から阻止域へは滑らかな曲線で遷移
※阻止域と通過域の境目は信号の大きさが半分になる点
(電力1/2)
ゲイン=1/√2(0.707) 3dB減衰する点
⇒遮断周波数
※実際の設計では、遮断周波数前後のゲイン1とも0とも見なせない範囲は遷移域として扱う
※理想の周波数条件に近づける
①通過域のゲインは平坦部分を長く
②境界付近は急峻
③阻止域のゲインは0に近く
_◇フィルタの特性カーブ
①バタワース特性(振幅平坦特性)
通過域平坦
遮断周波数での傾きなだらか
阻止域平坦
②チェビシェフ特性
通過域波状(リップル)
遮断特性急峻
阻止域平坦
③逆チェビシェフ特性
通過域平坦
遮断特性急峻
阻止域波状
④連立チェビシェフ特性
通過域波状
遮断特性最も急峻
阻止域波状
⇒楕円フィルタ(elliptic filter).
楕円関数を利用
⑤ベッセル特性
遅延平坦特性
遮断特性が緩やか
_◇位相遅れ
※入力信号 x(t)の位相を基準として、出力信号 y(t)の位相がどれだけずれているのかを、フィルタの位相(phase)という
位相θを各周波数 2πf で割ると時間となる
⇒位相は時間のずれと角周波数の積
※基準波形
θ=0
※進み位相
θ>0 時間軸に対して波形は左に移動
※遅れ位相
θ<0 時間軸に対して波形は右に移動
※入力信号の異なる周波数成分 f1, f2
それぞれ異なる時間 τ1, τ2だけずれて出力にあらわれる場合
⇒出力では周波数成分の位相関係の変化により波形ひずみが発生する
⇒ひずみを生じさせないためにはτ1=τ2でなければならない
ここでτと位相φ[rad]は
τ=φ/2πf
なる関係がある。2πf=ω[rad/sec]なので
τ=φ/ω
ともかける。
※遅延時間 delay time
τにマイナス符号をつけたもの
τ=一定値 なら「一定遅延特性を有する」
⇒位相φと各周波数ωもしくは周波数fが正比例しなければならない
⇒線形位相 linear phase
⇒線形位相でないフィルタでは出力に波形ひずみが生じる
※位相特性(偏角)の値は主値(-π~π)の範囲で表す
⇒ジャンプしたように見える部分は、2π平行移動した部分にあると考えるとよい⇒位相が一貫して遅れる(進む)ように見える
_◇線形位相
※位相特性θ(e^jωT)[rad]を角周波数ω[rad/Sec]で微分した値
dθ(e^jωT)/dω
⇒フィルタを通った後の出力が入力にくらべてどれくらいズレているかを示す物理量
位相=角周波数x時間
より、位相[rad]を角周波数[rad/秒]で微分すれば時間となる
※一定遅延
すべての角周波数ωの正弦波cos(ωT)に対して時間t0だけずれること
cos(ωT)⇒遅延⇒cos(ω(t-t0))
⇒入力と出力の波形変わらず歪なし
cos(ωt-ωt0)なので
時間のずれを位相量でみるとωt0[rad]
⇒角周波数ωに比例する
※線形位相
位相が角周波数ωに比例して、波形ひずみが生じない性質
_◇デジタルフィルタの全体構成
①アナログ入力
②帯域制限
アナログ信号をローパス・フィルタ(アンチエイリアシングフィルタ)によって標本化定理を満たすようにする。
⇒サンプリング周波数の半分の周波数成分しか含まれないようにする
③サンプル&ホールド
アナログ信号から離散信号への変換
④符号化
離散信号からデジタル信号への変換。ADCによる量子化と符号化。
⑤デジタル・フィルタ
ハードウエアでは乗算、加算、遅延回路の3要素。
⑥復号化
DACによりデジタル信号を離散信号に変換する。逆量子化。
⑦補間
離散信号からアナログ信号への変換。ローパス・フィルタ(補間フィルタ)を通す
_◇デジタル・フィルタの処理形態
①バッチ処理
いったんバッファメモリに蓄えて、蓄えた信号に計算処理を施して出力信号を得る。FFT処理が代表。全てのデータが出揃ってから計算する。複雑で、効率のよいアルゴリズムが使用できるが、リアルタイム的な処理には向かない
②リアルタイム処理
入力信号を実時間で逐次処理しながら出力信号を得る。現時点から見て過去の入出力データから現時刻の出力値を計算する。
_◇不規則雑音(ランダム雑音)の除去
※移動平均処理
平均値を計算する区間をずらしながら計算する
⇒ローパスフィルタでもある
⇒信号は比較的緩やかに時間変化し、雑音は不規則で時間的に速い(周波数が高い)
※通過域と阻止域
※ローパス処理後の雑音のみの取り出し
(処理後の雑音)
=(出力)-(入力を遅延時間遅らせたもの*利得)
利得(ゲイン):対象となる信号周波数の
_◇ドリフト雑音の除去
(測定データ)-(ドリフト雑音の推定値)
※ドリフト雑音の推定
移動平均
⇒ドリフト成分は低い周波数成分のため
※測定データから平均値を差し引く
⇒ハイパスフィルタ
_◇不規則雑音+ドリフト雑音の除去
※バンドパスフィルタ
⇒中間周波数に信号成分があるので、そこを取り出す
※ローパスフィルタをかけた波形にさらにハイパスフィルタを通す
_◇波形等化フィルタ
波形を復元するフィルタ
ゴーストを対象とする=ゴーストフィルタ
エコーを対象=エコー・キャンセラ
※ゴースト画像
主信号(直接波)G[k]と、ゴースト信号(間接波)g[k]が重なって起こる
⇒伝播伝搬ではマルチパス
⇒音声ではエコー
送信信号 x[k]
受信信号 y[k] = G[k] + g[k]
直接伝播が3サンプル遅れとすると
G[k] = x[k-3]
ゴーストはさらに3サンプルおくれ、かつ|α|<1だけ振幅がちいさくなるとすると
g[k] = α*x[k-6]
y[k]=x[k-3]+α*x[k-6]
とかける。
これに対して受信した信号y[k]を遅らせて-α倍した信号
-α*y[k-3]
を足す(ゴースト打ち消すが、ゴーストのα倍のゴースト残る)
さらにこれを遅らせて-α倍した信号を足す
α^2*y[k-6]
さらにまたこれを遅らせて-α倍した信号を足す
-α^3*y[k-6]
と、微弱なα^4倍の残留ゴーストがのこるだけで、主信号を取り出すことができる
_◇システム関数
例)
時間領域のブロック図
x[k] ┌─┐ ┌─┐
──┬┤T├┬┤T├┐
  │└─┘│└─┘│
  ▽1  ▽0.5▽0.25
  │  ┌┴┐ ┌┴┐
  └──┤+├─┤+├→y[k]
     └─┘ └─┘
差分方程式
y[k]=x[k]+0.5*x[k-1]+0.25*x[k-2]
⇒z変換すると
Y(z)=X(z)+0.5*z^-1*X(z)+0.25*z^-2*X(z)
=(1+0.5*z^-1+0.25*z^-1)*X(z)
⇒カッコ内をH(z)とおく
入力信号例
x[k]=0 k<0, k>2
=10≦k≦2
の方形波。z領域では
X(z)=1+z^-1+z^-2
Y(z) = H(z) * X(z)
がなりたつ
※H(z)を伝達関数とよぶ
出力信号の全体=伝達関数×入力信号の全体
_◇巡回形と非巡回形
※デジタルフィルタの入出力の一般形
y[k]=∑[m=0:M]{am*x[k-m]}+∑[n=1:N]{bn*y[k-n]}
※非巡回形(非再帰)
現時刻の出力 y[k] の計算で過去の出力 y[k-n] を用いない
右や帰還がない
y[k]=∑[m=0:M]{am*x[k-m]}
Y(z)=∑[m=0:M]am*z^-m*X(z)
H(z)=Y(z)/X(z)=∑[m=0:M]am*z^-m
x[k] ┌─┐ ┌─┐
──┬┤T├┬┤T├┐
  │└─┘│└─┘│
  ▽a0 ▽a1 ▽a2
  │  ┌┴┐ ┌┴┐
  └──┤+├─┤+├→y[k]
     └─┘ └─┘
※巡回形(再帰)
y[k]=∑[m=0:M]{am*x[k-m]}+∑[n=1:N]{bn*y[k-n]}
Y(z)=∑[m=0:M]{am*z^-m*X(z)}+∑[n=1:N]{bn*z^-n*Y(z)}
H(z)=Y(z)/X(z)=∑[m=0:M]am*z^-m/(1-∑[n=1:N]bn*z^-n)
x[k] ┌─┐ ┌─┐
──┬┤T├┬┤T├┐
  │└─┘│└─┘│
  ▽a0 ▽a1 ▽a2
  │  ┌┴┐ ┌┴┐ ┌─┐
  └──┤+├─┤+├→│+├────┬→y[k]
     └─┘ └─┘ └┬┘   ┌┴┐
              ↑    │T│
             ┌┴┐b1 └┬┘
             │+├─<├─┤
             └┬┘   ┌┴┐
              │    │T│
              │ b2 └┬┘
              └──<├─┘
_◇インパルス応答の継続時間による分類
※非巡回形デジタルフィルタ
インパルス波形を入力
⇒フィルタの係数に一致する応答出力が得られ、インパルス応答の継続時間は有限
⇒FIR (Finite Impulse Response)デジタル・フィルタ
※巡回形デジタルフィルタ
⇒インパルス応答の継続時間は無限
H(z)=∑[k=0:∞]{y[k]*z^-k}
⇒IIR (Infinite Impulse Resonse)デジタル・フィルタ
※例外的に巡回型でもFIRとなるものもある(特殊なもの)
_◇時間軸の伸張、圧縮
zをz^Lで置き換える
⇒時間軸をL倍に引き伸ばす
(応答出力の間に(L-1)個の信号値0を埋め込む)
のと等価
※時間軸をL倍に伸ばす⇒周波数軸は1/L倍に縮む
◆FIRフィルタ
※遅延回路の値を1サンプルずつずらしたものに係数を乗算し、それらの乗算結果を加算することで出力y[k]が得られる。
※インパルス応答の継続時間が有限⇒安定
y[k]=∑[m=0:M]{a_m*x[k-m]}
⇒積和演算
a_m: タップ係数、フィルタ係数、乗算係数
⇒フィルタの特性を定める
M: フィルタの次数(伝達関数の最高次数)
M+1:タップ数。フィルタ係数a_mの総個数
_◇伝達関数
Y(z)=∑[k=0:∞]{y[k]*z^-k}
=∑[m=0:M]{a_m*{∑[k=0:∞]x[k-m]*z^-k}}
H(z)=Y(z)/X(z)=∑[m=0:M]{a_m*z^-m}
_◇ブロック図
以下
┌────┐
┤Z^-1├
└────┘
┌─┐
┤Z├
└─┘
と略記する
※M次のFIRフィルタは通常、
M個の遅延回路
M個の加算回路
(M+1)個の乗算回路
で実現できる
①直接形
X(z)┌─┐ ┌─┐ ┌─┐
  ─┬┤Z├┬┤Z├┬┤Z├┐
   │└─┘│└─┘│└─┘│
   ▽a0 ▽a1 ▽a2 ▽aM
   │  ┌┴┐ ┌┴┐ ┌┴┐
   └──┤+├─┤+├─┤+├→Y(z)
      └─┘ └─┘ └─┘
※伝達関数の係数が直接乗算回路のタップ係数となる。
②転置形
X(z)
─┬────┬─────┬─────┐
 │    │     │     │
 ▽aM  ▽a2   ▽a1   ▽a0
 │┌─┐┌┴┐┌─┐┌┴┐┌─┐┌┴┐
 └┤Z├┤+├┤Z├┤+├┤Z├┤+├→Y(z)
  └─┘└─┘└─┘└─┘└─┘└─┘
③縦続形
伝達関数H(z)の因数分解
※次数Mが偶数
H(z)=H0*Π[k=1:M/2](1+a1(k)*z^-1+a2(k)*z^-2)
※次数Mが奇数
H(z)=H0*(1+a1(0)*z^-1)*Π[k=1:M/2](1+a1(k)*z^-1+a2(k)*z^-2)
ここでH0は利得定数
⇒各因数をそれぞれフィルタとして構成し、直列に接続する
⇒各因数を並べる順番:オーダリング自由
⇒1,2次だけでなく、より高次の多項式でも因数を作れる
X(z)┌─┐ ┌─┐
  ─┬┤Z├┬┤Z├┐
   │└─┘│└─┘│
   │   ▽a1 ▽a2
   │  ┌┴┐ ┌┴┐
   └──┤+├─┤+├┐
      └─┘ └─┘│
 ┌───────────┘
 │  ┌─┐ ┌─┐
 └─┬┤Z├┬┤Z├┐
   │└─┘│└─┘│
   │   ▽a1M▽a2M
   │  ┌┴┐ ┌┴┐  Y(z)
   └──┤+├─┤+├┤>→
      └─┘ └─┘H0
※例
H(z)=0.01+0.14*z^-1+0.49*z^-2+0.36*z^-3
①直接形
a0=0.01
a1=0.14
a2=0.49
aM=0.36
X(z)┌─┐ ┌─┐ ┌─┐
  ─┬┤Z├┬┤Z├┬┤Z├┐
   │└─┘│└─┘│└─┘│
   ▽a0 ▽a1 ▽a2 ▽aM
   │  ┌┴┐ ┌┴┐ ┌┴┐
   └──┤+├─┤+├─┤+├→Y(z)
      └─┘ └─┘ └─┘
②転置形
a0=0.01
a1=0.14
a2=0.49
aM=0.36
X(z)
─┬────┬─────┬─────┐
 │    │     │     │
 ▽aM  ▽a2   ▽a1   ▽a0
 │┌─┐┌┴┐┌─┐┌┴┐┌─┐┌┴┐
 └┤Z├┤+├┤Z├┤+├┤Z├┤+├→Y(z)
  └─┘└─┘└─┘└─┘└─┘└─┘
③縦続形
H(z)=0.01+0.14*z^-1+0.49*z^-2+0.36*z^-3
を因数分解する
H(z)=0.01*(1+z^-1)*(1+4*z^-1)*(1+9*z^-1)
それぞれ
1+z^-1, 1+4*z^-1, 1+9*z^-1を表す各段を直列に接続し、
最後に利得定数0.01を乗ずる
_◇周波数特性
※FIRフィルタの周波数特性
⇒伝達関数H(z)に
z=e^(j2πfT)=e^(jωT)
を代入すればよい
H(jωT)=∑[m=0:M]{a_m*e^-(jmωT)}
※これを極座標形式で表現すれば
H(jωT) = |H(jωT)| * e^(j∠H(ejωT))
⇒絶対値 |H(jωT)|がゲイン特性
⇒偏角∠H(ejωT)が位相特性
に相当する
H(jωT)=∑[m=0:M]{a_m*e^-(jmωT)}
にオイラーの公式を適用し
=∑[m=0:M]a_m*{cos(mωT)-j*sin(mωT)}
=∑[m=0:M]{a_m*cos(mωT)}-j*{-∑[m=0:M]a_m*sin(mωT)}
実部をHr(e^jωT), 虚部をHi(e^jωT)と書けば
(以下Hr, Hiと略す)
|H(jωT)|=√(Hr^2 + Hi^2)
∠H(ejωT)=arctan(Hr, Hi)
※周波数特性の計算例
移動平均
y[k]=(x[k]+x[k-1])/2
両辺をz変換すれば
Y(z)=(X(z)+z^-1*X(z))/2
伝達関数は
H(z)=(1+z^-1)/2
周波数特性は、zにe^jωTを代入すれば求められるので
H(e^jωT)=(1+e^-jωT)/2
=(((e^jωT/2)+(e^-jωT/2))*e^-jωT)/2
オイラーの公式を適用して整理すると
=cos(ωT/2)*e^-jωT/2
この式がサンプリング定理をみたすためには
0≦ωT<π
より
0≦ωT/2<π/2
なので cos(ωT/2)≧0
ゲイン特性は
cos(ωT/2)
ωが0付近ではゲインは1
ωが2π*fs/2=π/T付近では0
位相特性は
-ωT/2
負数なので遅れ位相
◇ローパス型FIRフィルタ
移動平均処理はローパスフィルタに相当する
_◇ハイパス型FIRフィルタ
移動差分処理はハイパスフィルタに相当する
例)
差分方程式
y[k]=(x[k]-x[k-1])/2
で記述されるFIRフィルタの周波数特性。
サンプリング間隔Tは0.01Sとする
差分方程式の両辺をz変換
Y(z)=(X(z)-z^-1*X(z))/2
よって
H(z)=(1-z^-1)/2
周波数特性はzにe^jωTを代入し
H(e^jωT)=(1-e^-jωT)/2=((e^jωT/2 – e^-jωT/2)/2)*e^-jωT/2
オイラーの公式を適用すると
H(e^jωT)=j*sin(ωT/2)*e^-jωT/2
さらに j=e^jπ/2 であることを利用して整理すると
H(e^jωT)=sin(ωT/2)*e^-j(π/2-ωT/2)
よって
ゲイン特性 sin(ωT/2)
位相特性 π/2-ωT/2
(位相は進み位相となる)
直流付近ではゲイン0、fs/2付近の高い周波数(π/T)付近ではゲイン1
_◇バンドパス型FIRフィルタ
ローパスとハイパスを組み合わせるとバントパス特性が実現できる
例)
p[k]=(x[k]+x[k-1])/2
y[k]=p[k]-p[k-1]
サンプリング間隔Tは1秒
差分方程式をz変換して伝達関数を求める
P(z)=(1+z^-1)*X(z)/2
Y(z)=(1-z^-1)*P(z)
H(z)=Y(z)/X(z)=(1-z^-2)/2
周波数特性は
H(e^jωT)=sin(ωT)*e^j(π/2-ωT)
よって
ゲイン特性 sin(ωT)
位相特性 π/2-ωT
直流付近と、0.5Hz付近ではゲイン0.0.25Hz付近ではゲイン1のバンドパス特性
※サンプリング周波数を1Hzとすると正規化周波数とおなじ
※接続順序を変えても伝達関数は同じ。
_◇バンドエリミネート型FIRフィルタ
ローパスフィルタの時間軸を操作してバンドエリミネート特性を実現することもできる
例)
1次のローパスフィルタ
H(z)=(1+z^-1)/2
H(z^2)の特性
⇒時間軸が1/2に縮む⇒周期性によりバンドエリミネート特性を示す
_◇線形位相 FIRフィルタ
※FIRフィルタは完全な直線位相特性の実現が可能
※無歪伝送条件
x[k]=x[k-k0]
ここでk0は定数
z変換すれば
Y(z)=z^-k0*X(z)
H(z)=Y(z)/X(z)=z^-k0
周波数特性をもとめるため z=e^jωTを代入し
H(z)=e^-j*k0*ω*T
ゲイン特性 |e^-j*k0*ω*T|=1
オールパス特性(全周波数で1)
位相特性 ∠e^-j*k0*ω*T=-k0*ω*T
ωに比例するので線形位相
※遅延時間をτ(e^jωT)で表し、
τ(e^jωT)=d∠H(e^jωT)/dω
と定義し
τ(e^jωT)=-k0*ω*T
※線形位相フィルタ=タップ係数が左右対称
⇒インパルス応答波形が時間の対称中心に対して線対称
⇒2個の乗算を1個に集約可能
⇒次数が偶数の場合
乗算回路 (N/2) + 1
例) H(z)=a0+a1*z^-1+a2*z^-2+a1*z^-3+a0*z^-4
⇒次数が奇数の場合
乗算回路 (N+1)/2
例) H(z)=a0+a1*z^-1+a2*z^-2+a2*z^-3+a1*z^-4+a0*z^-5
X(z)┌─┐    ┌─┐
─┬─→┤Z├─┬─→┤Z├─┐
 ↓  └─┘ ↓  └─┘ │
┌┴┐    ┌┴┐     │
│+├┐   │+├┐    │
└┬┘│   └┬┘│    │
 ↑ │┌─┐ ↑ │┌─┐ │
 └←┼┤Z├←┴─┼┤Z├←┤
   ↓└─┘   ↓└─┘ ↓
   ▽a0    ▽a1  ▽a2
   │      │    │
   │     ┌┴┐  ┌┴┐
   └─────┤+├──┤+├→Y(z)
         └─┘  └─┘
◆IIRフィルタ
※一般に急峻な遮断特性を実現するにはIIRフィルタがFIRフィルタより有理
⇒FIRフィルタよりはるかに少ない次数で実現できるが、FIRフィルタでは問題とならなかった安定性が問題となる
※インパルス応答の継続時間が無限
⇒出力が入力に戻されるフィードバック構成
例)1次のIIRフィルタ
⇒ステップ応答類似
x[K]
 │
 ▽2
┌┴┐
│+├←──┐
└┬┘   △0.5
 │ ┌─┐│
 ├─┤T├┘
 │ └─┘
 ↓y[K]
※安定なフィルタの条件
インパルス応答の出力y[k]に対して
lim[k→∞]|y[k]|→0
※一般的なIIRフィルタ
y[k]
= ∑[m=0:M]a_m * x[k-m] + ∑[m=1:N]b_m * y[k-n]
_◇伝達関数
y[k]=∑[m=0:M]a_m * x[k-m] + ∑[m=1:N]b_m * y[k-n]
をz変換して
Y(z)=Y(z)*∑[m=1:N]b_m*z^-n + X(z)*∑[m=0:M]a_m*z^-m
H(z)=Y(z)/X(z)=A(z)/B(z)
A(z)=∑[m=0:M]a_m*z^-m
B(z)=1-∑[m=1:N]b_m*z^-n
a_m, b_m タップ係数(フィルタ係数)
M 非再帰部分の次数
N 再帰部分の次数
⇒分母多項式B(z)の定数項は正規化し、必ず1とすること
⇒b_mが全て0であるときFIRフィルタに一致する
※分母多項式B(z)があるので、B(z)=0を満たすzの極の存在により安定性が問題となる
_◇安定性
※安定なフィルタの条件
インパルス応答の出力y[k]に対して
lim[k→∞]|y[k]|→0
⇒高次のIIRフィルタの分母多項式B(z)は1次と2次の式を因数として因数分解される
⇒1時と2次のフィルタについて安定性が分れば分る
※1次のIIRフィルタ
y[k]=b1*y[k-1]+x[k]
当然b1は0ではない
⇒x[k]=δ[k], y[-1]=0とすると
y[k]=(b1)^k
|b1|>1 出力発散する
|b1|=0 減衰することなく発振
|b1|<0 0に収束する
H(z)=1/(1-b1*z^-1)
極をzpとすると
1-b1*zp^-1=0 より zp=b1
⇒安定条件は
|zp|<1
と書き換えられる。
⇒極の絶対値が1より小さい
※2次のIIRフィルタ
差分方程式
y[k]=b1*y[k-1]+b2*y[k-2]+x[k]
H(z)=1/(1-b1*z^-1)-b2*z^-2)
⇒安定条件は
|zp1|<1
|zp2|<1
⇒極の絶対値が1より小さい
とやはり書ける
B(z)=(1-b1*z^-1)-b2*z^-2)
と2次式なのでB(z)=0となる極は、実数でも複素数でもありえる
判別式 D=b1^2+4*b2
なので
1) D≧0の実根の場合は、
zp1=α、zp2=βとおいて
H(z)=h1/(1-α*z^-1)+h2/(1-β*z^-1)
ただし、
h1=1/(1-α^-1*β)
h2=1/(1-α*β^-1)
と部分分数に分解できる
また、2次方程式の解と係数の関係から
α+β=b1, α*β=-b2
インパルス応答は
y[k]=h1*α^k+h2*β^k
と求まる
2) D<0の複素根の場合は
共役な複素根を re^jθ re^-jθとして
やはりH(z)を部分分数で表す
H(z)=h1/(1-re^jθ*z^-1)+h2/(1-re^-jθ*z^-1)
ただし、
h1=1/(1-re^j2θ)
h2=1/(1-re^j2θ)
ここでオイラーの公式より
re^jθ+re^-jθ=r(cosθ+jsinθ)+r(cosθ-jsinθ)
=2*r*cosθ
2次方程式の解と係数の関係から
2*r*cosθ=b1
また
re^jθ * re^-jθ=r^2=-b2 つまり -b2≧0
インパルス応答
y[k]=r^k(sin((k+1)*θ))/sinθ)
インパルス応答y[k]がk→∞で0に収束するような条件を求めると
実根の場合
D≧0より b2≧-(b1^2)/4
b2<-b1+1
b2<b1+1
複素根の場合
r^kが収束すればよいので
D<0よりb2<-(b1^2)/4
r^2=-b2<1 上の-b2≧0とあわせ
0≦-b2<1
0≧b2>-1
※横軸にb1, 縦軸にb2をとると
b2=-1, b2=-b1+1, b2=b+1の三本の直線に囲まれた二等辺三角形の中が安定領域となる
※Z平面では
B(z)=0を満たす根が
単位円の内部は安定。単位円の円周上は発振、外部は発散
_◇安定性の判別
①yk=0.5*y_k-1 + 0.24*y_k-2 + x_k
差分方程式
y[k]=b1*y[k-1]+b2*y[k-2]+x[k]
とくらべ
b1=0.5
b2=0.24
これは、b1-b2平面で安定領域にある
②yk=1.4*y_k-1 – 0.8*y_k-2 + x_k
差分方程式
y[k]=b1*y[k-1]+b2*y[k-2]+x[k]
とくらべ
b1=1.4
b2=-0.8
これは、b1-b2平面で安定領域にある
③yk=1.8*y_k-1 – y_k-2 + x_k
差分方程式
y[k]=b1*y[k-1]+b2*y[k-2]+x[k]
とくらべ
b1=1.8
b2=-1
これは、b1-b2平面で安定領域にない
④yk=-0.8*y_k-1 – 2*y_k-2 + x_k
差分方程式
y[k]=b1*y[k-1]+b2*y[k-2]+x[k]
とくらべ
b1=-0.8
b2=-2
これは、b1-b2平面で安定領域にない
_◇ブロック図
①直接型構成
伝達関数
H(z)=H1(z)*H2(z)=A(z) * 1/B(z)
と考え、A(z)部分と B(z)部分を直列に接続したもの
②標準型構成
直接型をB(z)が先、後からA(z)と順序を交換すると
遅延回路を共有できるようになる。
⇒遅延回路の個数を最小化したもの
⇒最小実現のシステム構成
 ↑y[k]
┌┴┐ ┌─┐ ┌─┐ ┌─┐
│+├←┤+├─┤+├─┤+├──┐
└┬┘ └┬┘ └┬┘ └┬┘  │
 △a0 △a1 △a2 △a3 △a4
 │┌─┐│┌─┐│┌─┐│┌─┐│
 ├┤T├┼┤T├┼┤T├┼┤T├┤
 │└─┘│└─┘│└─┘│└─┘│
 ↑   ▽b1 ▽b2 ▽b3 ▽b4
┌┴┐ ┌┴┐ ┌┴┐ ┌┴┐ ┌┴┐
│+├←┤+├─┤+├─┤+├─┤+│
└┬┘ └─┘ └─┘ └─┘ └─┘
 ↑x[k]
③並列型構成
伝達関数H(z)を部分分数に展開し、
その結果にもとづき、1次と2次のIIRフィルタを並列接続する。
④継続型構成
H(z)を分子、分母でそれぞれ因数分解
     A1(z)A2(z)
H(z)=----------=H1(z)H2(z)
     B1(z)B2(z)
⇒1次と2次のIIRフィルタを直列接続する形にできる
※ペアリング
各因数A1,A2,B1,B2の組み合わせ方
⇒極と零点の組み合わせ
※オーダリング
直列に接続する順序
⑤転置型構成
標準型の入出力を逆転したもの
例)
伝達関数
H(z)=A(z)/B(z)
A(z)=0.043716-0.072699*z^-1+0.104426*z^-2-0.072699*z^-3+0.043716*z^-4
B(z)=1-3.0159*z^-1+3.7965*z^-2-2.2880*z^-3+0.55931*z^-4
_◇1次のIIRフィルタの極と零点
伝達関数
H(z)=H0*(1+a1*z^-1)/(1-b1*z^-1)
|b1|<1
H0:利得定数
①LPF
a1=1
H(z)={(1-b1)/2} * {(1+z^-1)/1-b1*z^-1)}
※ω=0 (z=e^jωT=1)における振幅を1として正規化)
⇒零点 ガウス平面の -1
⇒極 b1
②HPF
a1=-1
H(z)={(1+b1)/2} * {(1-z^-1)/1-b1*z^-1)}
※ω=π/T (z=e^jωT=-1)における振幅を1として正規化)
⇒零点 ガウス平面の 1
⇒極 b1
③APF
H(z)=(-b1+z^-1)/(1-b1*z^-1)
_◇2次のIIRフィルタの極と零点
伝達関数
H(z)=H0*(1+a1*z^-1+a2*z^-2)/(1-b1*z^-1-b2*z^-2)
b1,b2は、横軸にb1, 縦軸にb2をとったとき
b2=-1, b2=-b1+1, b2=b+1の三本の直線に囲まれた二等辺三角形の中が安定領域
H0:利得定数
①LPF2
a1=2, a2=1
H0=(1-b1-b2)/4
⇒ガウス平面の-1が二重零点となる
⇒極は単位円の内部に2つ(複素共役)
※ω=0の振幅を1に正規化
②HPF2
a1=-2, a2=1
H0=(1+b1-b2)/4
⇒ガウス平面の1が二重零点となる
⇒極は単位円の内部に2つ(複素共役)
※ω=π/Tの振幅を1に正規化
③BPF2
a1=0, a2=1
H0=(1+b2)/2
⇒ガウス平面の1と-1が零点となる
⇒極は単位円の内部に2つ(複素共役)
※振幅ピークの振幅を1に正規化
④BEF2
|a1|<2, a2=1
H0=(1-b2)/2
a1=-2*b1/(1-b2)
⇒単位円上に零点(複素共役)
⇒極は単位円の内部に2つ(複素共役)
⑤APF2
H(z)=(-b2-b1*z^-1+z^-2)/(1-b1*z^-1-b2*z^-2)
⇒単位円の外に零点(複素共役)
極までの距離をrとすると零点は1/rの位置にある
r=√(-b2)
⇒極は単位円の内部に2つ(複素共役)
_◇IIRフィルタにより正弦波を作る例
周波数f0[Hz], ω0=2*π*f0 [rad/s]の音y(t)
y(t)=sin(ω0*t)
①デジタル化する
サンプリング間隔T[s], t=k*T (kは整数)とすると
y[k] = y(k*T) = sin(k*ω0*T)
⇒このsin部分を算術で作り出す
②k=(k-2)+2と分解し、加法定理を適用する
sin(α+β)=sinαcosβ+cosαsinβ
y[k] = sin(k*ω0*T)
=sin{(k-2)*ω0*T+2*ω0*T}
=sin{(k-2)*ω0*T}cos{2*ω0*T}+cos{(k-2)*ω0*T}sin{2*ω0*T}
③2*ω0*Tなる項に2倍角の公式を適用する
sin2θ=2sinθcosθ
cos2θ=2*cos^2θ-1
y[k]=
sin{(k-2)*ω0*T}*{2*cos^2(ω0*T)-1}
+cos{(k-2)*ω0*T}*2*sin(ω0*T)*cos(ω0*T)
=2*cos(ω0*T)*[sin{(k-2)*ω0*T}*cos(ω0*T)
+cos{(k-2)*ω0*T}*2*sin(ω0*T)]-sin{(k-2)*ω0*T}
④上記式の[]内は、加法定理を逆に用いて
sin{(k-2)*ω0*T+ω0*T}=sin{(k-1)*ω0*T}
と変形できる
⑤差分方程式とするために
sin{(k-1)*ω0*T}はy[k-1]
sin{(k-2)*ω0*T}はy[k-2]
で置き換えれば
y[k]=2*cos(ω0*T)*y[k-1]-y[k-2]
を得る。
⇒2*cos(ω0*T)は定数なので、単純なフィルタとなる
⑥k=0, k=1の初期値として
y[0]=sin(0)=0
y[1]=sin(ω0*T)
を与え、以下k≧2について差分方程式を計算すれば正弦波波形が得られる。
    ┌─┐      y[k]
┌───┤+├────────┬→
│   └┬┘        │
│    │         │
│    △2cos(ω0T)│
│┌─┐ ↑ ┌─┐     │
└┤T├←┴←┤T├─────┘
 └─┘   └─┘
初期値      初期値
sin(ω0T)  0
◆周波数スペクトル分析
_◇FFT
Fast Fourier Transform
◆過渡応答
※信号処理でもでだしのところには過渡応答が生じる
◆雑音除去
◆信号発生
◆信号検出
◆特徴抽出
◆数値微分
※測定データ x1..xN ⇒ y1..yN に対して
(y_n+1 – y_n)/(x_n+1 – x_n)
を計算すると、多くの場合、ノイズが強調され所望の結果が得られない。
⇒2点でなく、もう少し広い範囲のデータを見る
⇒広い範囲を多項式で近似する
⇒その多項式の微分を、もとの関数の微分とする
※変数変換により、範囲の中心を0に移動することで
⇒近似式の係数が微分値に等しくなる
⇒測定値 y_n+1に b_liという重みをつけて足しこむ
⇒多項式近似の係数 ai を求めることができる
⇒b_liは測定間隔んと近似に使うデータ数 m が固定なら不変
◆増幅
◆誤り訂正
_◇伝送路で発生する誤り
①ノイズ、フェージングによるビット誤り
②バーストエラー(ブロック状にドロップアウトする。例えばビル陰に入った)
_◇有限体
四則演算が有限個の数の集合内で規定された体系。
※1ビットのデータの枠組みで計算が成り立つ世界をGF(2)とする。
G: ガロア
F: Field (体)
ガロア体
※GF(2)で加算を定義するとXORとなるが、減算と区別できない。
_◇ガロア拡大体
X^2+X+1=0
はGF(2)では根を持たないが、仮に上の式を満たすαがあると仮定し、2個のGF(2)の変数を使ってガロア拡大体GF(2^2)を考える。
Z=x+αyx,y:GF(2)
この場合、
X^2+X+1=0
は、αとα^2の2つが根となる。
※1つのガロア拡大体状の代数方程式を解いて得られる複数の解も、複素数と同じく共役根と呼ばれる。
GF(2^2)は
0,1,α,1+α
の4つの元から成り立っているが
α^2+α+1=0
なので
α^2=α+1
となり、
0,α^0,α^1,α^2
このとき、α^3以上は巡回し、独立した元ではない。
_◇原始多項式
N個のGF(2)の変数を使ってガロア拡大体GF(2^N)を作る場合、N次多項式が必要になる。これはαのベキ乗で表現できる独立な2^N-1個のGF(2^N)の元が存在するものでなければならない。
※αのことをGF(2^N)の原始元とよぶ
次数 GF原始多項式
2 X^2+X+1
3 X^3+X+1
4 X^4+X+1
5 X^5+X^2+1
6 X^6+X+1
7 X^7+X+1
8 X^8+X^4+X^3+X^2+1
9 X^9+X^4+1
10 X^10+X^3+1
11 X^11+X^2+1
12 X^12+X^6+X^4+X+1
13 X^13+X^4+X^3+X+1
14 X^14+X^10+X^6+X+1
15 X^15+X+1
16 X^16+X^12+X^3+X+1
_◇ガロア拡大体の四則演算
GF(2^2)の2つの変数x,yの加算
GF(2^2)の元となる
減算、加算に置き換えられる
乗算
x,yは0またはαのべき乗で表されるので乗算は指数の足し算となる。ガロア拡大体の元はαのべき乗で巡回するので、掛け算の結果もGF(2^2)に含まれる
除算
やはりべき乗の差となるので、GF(2^2)に含まれる。
※複素数式にX+Y*αと記述せず、ベクトル(X,Y)として表現することも可能。加算はベクトル和。掛け算はベクトルマトリクス演算となる。
※用途に応じて使い分ける
関数表現
ベキ表現
ベクトル表現
※ガロア拡大体GF(2^4)での4則演算例
16個の元があり、原始多項式はX^4+X+1
(1010)と(0110)の掛け算
(1010)=α^3+α^1=(加算表参照)=α^9
(0110)=α^2+α^1=(加算表参照)=α^5
よって
(1010)*(0110)=α^(9+5)=α^14=(加算表参照)=(1001)
◆画像処理
_◇ノイズ
①ブロックノイズ
MPEG,JPEG等でブロックの境目がわかるようなブロック状のノイズ
⇒ビットレートが十分でない場合に起こる。
②モスキートノイズ
⇒蚊の大群が飛び交っているようなもやもやしたノイズ
⇒激しく動く物体にまとわりつくように発生する。
③エッジノイズ
⇒テロップなど境界線がハッキリしたものの周囲に現れる。
_◇マッハ・バンド
Mach band
※グラデーションなど階調をなだらかに変化させる表現
⇒階調の境界線が目だって見える現象が発生
⇒マッハの帯
※人間の視覚
階調の境界を強調する特性を備える
例)黒画面→グレースケールの画面(段階的に明るくする)→白画面を接する形で配置
⇒白画面の領域でグレースケール画面と接する場所が,帯状により白が強調されて見える。
※錯視が原因
_◇ウィーナフィルタ(最小2乗フィルタ)
※復元フィルタで得られた画像における復元誤差を測る評価基準
E∥ε∥^2=E[{f(x,y)-f(x,y)}^2]
※原画像と復元画像との間の平均2乗誤差を考える。
※平均2乗誤差を最小にする画像を得るフィルタ
⇒ウィーナフィルタという。
_◇輪郭抽出
微分処理
信号差を求める
⇒信号の変化店を抽出する
◆オーバーサンプリング、サンプリングレート変更
サンプリングによって発生する量子化ノイズを広帯域幅に均等に分散させ、信号帯域のノイズパワーを減少させることができる。
_◇デジメーション
デジタル信号のサンプリングレートを下げる処理
※単純なデジメーションではオーバサンプリングの効果を損なう
⇒デジメーションフィルタをデジメーションの前段に設ける
デジメータ decimator
※ポリフェーズフィルタ
⇒シリアルパラレル変換の利用
_◇補間処理
サンプリングレートを高める処理
◆波形生成
_◇正弦波生成
①テーブル参照法
②再帰フィルタ法
単位遅延の漸化式から求める
y_n = a1 * y_n-1 + a2 * y_n-2
a1 = 2 * cos(ω*t)
a2 = -1
③テイラー展開法
◆音声符号化
_◇符号化方式
①波形符号化
通常の時間対電圧波形などを、可逆もしくは非可逆に圧縮
②スペクトル符号化
音源+音響管モデルで符号化する
③ハイブリッド符号化
スペクトル符号化で表現しきれない部分を波形符号化する
⇒CELP (Code Excited Linear Prediction)
A-b-S法(Analysis-by-Synthesis)
_◇LPC
Linear Prediction
線形予測分析
◆画像符号化
_◇MH符号化
_◇JPEG
Joint Photographic Experts Group
◆動画符号化
_◇MPEG
Moving Picture Experts Group
◆プログラム例
_◇量子化誤差
quantizedError … 正弦半波の8ビット量子化と誤差
gnuplotによるプロット結果…quant.png