ディジタルフィルタは、難しい理論から、当り前の技術にならなければならない。コンデンサと抵抗と OP アンプしかない時代においては、信号を数値として扱うディジタルフィルタは、高価な夢の技術であった。 ところが今、ディジタルコンピュータの普及と高度化によってその環境はあり得ないほど一変してしまった。
オーディオへの応用で、私は 1999 年に、PC の能力で 44.1 kHz のステレオのサンプルに対してリアル タイムに 400 個の係数の積和処理が可能であることを確認したことがある。そのとき思ったことは、 数100個の FIR が使えるときに私の知識は、まだ数タップのフィルタを扱うことしかできていない。 大量のタップ数を扱う技術は、古典的なフィルタの設計技術から別の領域の別ものに変るだろう、 ということである。
44.1 kHz ステレオサンプルには、1 サンプル 1 ch あたり 11.34 usec の時間が許される。11マイクロ秒に できる積和数は、巨大である。Web で FIR フィルタを PC の MMX 命令を使って 1 タップあたり 0.625 サイクルで行うソースが公開されていた。http://www.bdti.com/ これを使えば 1GHzクロックあたり、 18,140個のタップを処理できる。それは 0.4 秒のインパルス応答の長さに相当し、 ほとんどの残響処理がすでに可能であることを示している。
2016/4/6 gnuplot を使用する等ソース、filter.tgz を置きました。
1. 入力信号とフィルタと出力信号
2. インパルス応答
3. フィルタの縦続接続
4. 離散フーリエ変換 (DFT)
5. Z 変換
6. FIR フィルタと IIR フィルタ
7. 逆フィルタ
8. Convolution 定理
9. 相互相関関数と自己相関関数
10. FIR フィルタは全ゼロ型、IIR は全極型
11. 多項式展開
12. 発散をしない IIR フィルタへの修正
13. DFT と多項式展開の簡単なプログラム
14. FIR フィルタの例
15. 振動現象と複素回転
16. 逆フィルタと画像復元
17. 超解像
18. FIR と IIR の定義について
線形系とは、入力の定数(スカラー)倍と和とが出力の定数倍と和をもたらすことである。 x(i) が ax(i) になると出力は、y(i) から ay(i) となる。和 x(i)= x1(i) + x2(i)を入力すると、 出力はそれぞれの出力の和 y(i)= y1(i) + y2(i) になる。
現実のシステムに、普通の入力値の 1000 倍もの入力をいれるとシステムは全く違う動作をするものである。 大きすぎる入力では波形はクリップされる。逆に小さすぎると雑音に埋もれるだけでなく、非線型性によって 消えるのが普通である。ある限定された範囲で線形とみなしてよい領域があり、線形性は、そこだけで成立する。
デバイスは、信号を入力し信号を出力する。入出力が時間に離散的でないなら信号は時間関数であり、 デバイスは、関数の関数、汎関数であるが、ディジタル信号処理では入出力は数列、デバイスは数列の関数である。 数列をベクトルとみると、デバイスはベクトルを入出力する行列である。
デバイスの定常性を表すのに入出力の添字を一定値シフトしても成立するという、"シフト不変"という言葉が 使われる。y(i)= f(x(i)) のとき y(i+j)= f(x(i+j))。線形シフト不変な系は、インパルス応答だけで記述できる。
離散的な時間を使う前、システム理論は、どのような数学的道具立てを使ったか。交流理論では交流に対する 抵抗がコイル(インダクタンス)では周波数に比例する 2πfL、コンデンサでは周波数に反比例した 1/(2πfC) という抵抗と考えて回路を扱うことができた。過渡特性を解析するためにラプラス変換が使われ、周波数領域 を扱うのにフーリエ級数展開、フーリエ変換が使われてきた。
ラプラス変換では、最初から信号の変換形式を代数式で表現した。1 は、インパルス、1/s は、ステップ関数 である。デバイスと出力とも式で表示する。s は単なる複素数だった。積分要素は、1/s 、微分要素は s を、 1次遅れ要素(*)は 1/(s+a) で表せた。ラプラス変換は、変換後の式の変形だけで、実用的な問題を扱えるよう にした。インパルス応答、ステップ応答の代数式でシステムが記述でき、その対応する時間波形も容易に類推できた。
ラプラス変換は、 f(t) exp(-st) の t=0〜∞の定積分である。
∞ F(s)= ∫ f(t) exp(-st) dt t=0f(t) が面積1で幅が0極限のインパルス(クロネッカーのデルタ関数)ならその積分は、t=0 の exp(-st)の値 1 になる。f(t) がステップ関数 f(t)=1 (t>=0) なら、[-1/s exp (-st)]_0^∞= 1/s になる。f(t)= exp(-t/T) なら exp(-t/T -st) の 0〜∞の積分は、[-1/(1/T +s) exp(-t/T-st)]_0^∞ = 1/(1/T +s) = 1/(s+a), (a= 1/T) となる。
フーリエ級数展開は、連続関数である波形 f(t) を基本周波数 f0= 1/T0 の整数倍の周波数をもつ三角関数 cos と sin の系列の和、級数で表わした。
T0 T0
A(n)= ∫ f(t) cos(nwt) dt, B(n)=∫ f(t) sin(nwt) dt,
t=0 t=0
∞
f(t)= Σ{ A(n)cos(nwt) + B(n)sin(nwt) }
n= 0
ここで -B(n) を B(n) と書き直して、
T0 T0
A(n)= ∫ f(t) cos(nwt) dt, B(n)= -∫ f(t) sin(nwt) dt,
t=0 t=0
∞
f(t)= Σ{ A(n)cos(nwt) - B(n)sin(nwt) }
n= 0
この A(n), B(n) のペアは、同じ周波数の振幅 H(n)= √(A^2 + B^2) と位相 arg(n)= arctan(B/A) を考えること
ができ、複素数 F(n)= A(n) + j B(n) (j=√-1) をなす。オイラーの法則 exp(jwt)= cos(wt) + j sin(wt) から、
T0
F(n)= ∫ f(t) exp(-jnwt) dt
t=0
∞
f(t)= Σ Re{F(n) exp(njwt)}
n= 0
なぜなら、F(n) exp(njwt) = { A(n)cos(nwt) - B(n)sin(nwt) } + j { B(n)cos(nwt) + A(n)sin(nwt) }
フーリエ級数展開した周波数軸上の離散数列は、複素数列 F(n) であり、もとの実数波形 f(t)は、なにかの実部
を見ているといってよい。逆変換の式の中に複素数の実部をとる式 Re{} は、そういう意味にとれる。それなら、
原波形も実数に制限する理由はなく、複素数に一般化してよいとして、逆変換の実部をとるのをやめると次の
フーリエ変換になる。順変換と逆変換がほとんど同じ式で、乗算する exp() が共役であるだけの違いになる。
フーリエ変換は、連続から連続への変換であり、ラプラス変換とおなじ積分(複素数 s を周波数に比例する純虚数 s=jw として) を半無限0〜∞の積分から、t= -∞〜∞の積分に変えている。これはフーリエ級数展開を複素表示し て周波数を連続にしたものでもある。フーリエ変換は、次の定積分を使って、複素数も許す時間関数 f(t) を 周波数の関数 F(jw) へ変換する。
∞
F(jw)=∫f(t)exp(-jwt)dt
t=-∞
∞
f(t)=∫F(jw)exp(jwt)dt
w=-∞
ラプラス変換が半無限0〜∞の積分であること、 f(t) の t＜0 の部分を扱わない理由は、現在と過去にだけ依存する
物理的に実現可能な系に関係している。現実の系では未来の値を使用するフィルタは存在できないという物理的な
因果律からの制約がある。システムを時間関数のインパルス応答として表すとき、t＞＝0 だけを考え、積分範囲を
0〜∞ だけにするのである。つまり、ラプラス変換は、時間関数は任意の時間波形ではなく、システムのインパルス
応答として存在するものだけを扱うのである。次章のインパルス応答 h(k) において k= 0〜∞とするフィルタを
因果的 (causal) という。
離散的フーリエ変換 (DFT) は、時間、周波数ともに離散的で、複素数列から複素数列への変換である。フーリエ変換 のほとんどの性質を受け継いでいるし、有限数値の有限列の積和であるから結果も有限であり、収束の問題がない。
離散フーリエ変換した周波数領域で眺めると、フィルタを通した出力信号 y の周波数特性 Y(f) は、入力信号 x の 周波数特性 X(f) とフィルタの周波数特性 H(f) との乗算になる。(Convolution 定理) これは、ラプラス変換、 フーリエ変換、Z 変換にも共通した性質である。信号にフィルタを掛けるという言葉がぴったりしているように、 その出力信号の周波数特性は、入力信号の周波数特性にフィルタの周波数特性を乗算したものである。
(*) 抵抗 R[Ω] とコンデンサ C [F]との直列回路に電圧信号を入力し、コンデンサの両端電圧を出力とする系。 出力は、電流の時間積分 y = 1/C ∫idt、入力は、出力と抵抗の両端電圧(電流と抵抗 R の積)の和 x= y+Ri である。 x = y + CR y' ここで y'/y を s と置けば、この系は、1/(1+sCR) である。時定数 CR は時間の単位 [ΩF=sec] をもつ。 たとえば、1kΩと0.1μFは、0.1msec の時定数をもつ。s に j2πf をいれ f= 1/(6.28CR)となる 1.6 kHzで折れ曲がり、 それ以下では平坦、それ以上で -6dB/octave の周波数特性をもつ。
コンデンサは、電流の時間積分/Cが電圧、インダクタンスでは、電流の微分に L を乗算したものが電圧だった。 定常な交流では、インピーダンスは、インダクタンスには 2πfL、コンデンサには 1/2πfC ですむ。
入力 x(i)がインパルス {1,0,0,0,...} であるなら、フィルタの出力 y(i) は、フィルタのインパルス応答 h(i) である。 なぜならフィルタの動作を表すコンボルーション、
N-1
y(i)= Σ h(k)x(i-k) (1)
k=0
において x(i)= 1 (i==0), 0 (else) とすると、Σは、i==k のときの h(k) だけを値にもつから y(i)= h(i) となる。
周波数領域で考えると、入力信号のインパルスの周波数特性は、f によらずに X(f)= 1 であり、
フィルタのインパルス応答 h(k) の周波数特性 H(f) がそのまま出力 Y(f) になる。
なお、フィルタをベクトルθ= (h(0),h(1),...,h(N))とし、Χ(i)= (x(i),x(i-1),x(i-2),...x(i-N))^T として、 積和を内積で記述することもある。 y(i)= θ・Χ(i)
そのため縦続接続のインパルス応答 h12(k) は、各インパルス応答 h1(k), h2(k)のコンボリューションである。 縦続接続の周波数特性は、各フィルタの周波数特性の乗算である。縦続接続の順は交換可能である。
H12(f)= H1(f)H2(f) h21(i)= h12(i)= Σ h2(k)h1(i-k) k 短いインパルス応答の FIR フィルタの縦続接続のインパルス応答は、簡単な計算でできる。 (例1) 11型 2段の縦続接続は、121型である。 1,1 + 1,1 ------- =1,2,1 (例2) 1,-1型 2段の縦続接続は、1,-2,1 型。 1,-1 -(1,-1) --------- 1,-2, 1 (例3) 121 と 11 の2段の縦続接続は、1,3,3,1 型。 (例4) 11 と 1,-1 の2段の縦続接続は、1,0,-1型。 (例5) 121 と 1,-1 の 2 段の縦続接続は、1,1-1,-1 型。
N-1
X(k)= 1/√N Σ x(i) w^{-ik}
i=0
N-1
x(i)= 1/√N Σ X(k) w^{ik}
k=0
である。変換は、直交(orthogonal)であり、逆変換は、原信号 x(i) を再現する。
1/√N によって正規化された直交変換は、ユークリッド距離を変えない座標変換である。
ΣX(k)＾2 = Σ x(i)^2
原信号 x(i) が実数のとき、 X(k)= X(-k)^* である。そのため、その実部は偶関数 (X_r(N/2)は値をもつ。)、虚部は奇関数である(X_i(N/2)= 0)。
X_r(k) = X_r(-k) = X_r(N-k)
X_i(k)= -X_i(-k) = -X_i(N-k)
そのため、原信号 x(i) N個の実数を実部 N/2+1 個と虚部 N/2-1 個の N 個で表す。 X(k)の実部と虚部は、それぞれ、
X_r(k)= Σ x(i) cos(-2πik/N)
X_i(k)= Σ x(i) sin(-2πik/N)
である。
DFTは、Z変換のような解析的な設計手法ではなく、数値計算による周波数特性の分析手段と考えられた。 1965 年に Cooley と Turkey によって高速演算 FFT (Fast Fourie Transform)が開発され、コンピュータと DSPの処理能力の向上に伴って、リアルタイムの処理にも利用されてきた。
例えば、オーディオの符号化において AC3 では、512 点窓の(256点) MDCT のために 128 点複素 FFT が使用され、MPEG-2 AAC (NBC) では、窓点数 1024点 (又は 128点) MDCT のために窓点数の1/4の複素 FFT が使用される。
オーディオの符号化ではサブバンドフィルタやDFTによって周波数軸に変換して粗く量子化し、符号化すること が基本的な方法である。MPEG Audio Layer II や mp3 (MPEG Audio Layer III) では、32のサブバンドに 分解するポリフェーズフィルタバンクが使われ、mp3 では、サブバンド標本をさらに 36 点 DCT変換する。 それ以降、さらに変換時間幅を大きくして、周波数分解能を上げる方が圧縮率を上げることができることが 明らかになってきた。DFT 自体が直接にオーディオの符号化の変換として使われたのではなく、MDCT という 変換に DFTが使われる。
N-1
Z(z)=Σ h(i) z^{-i}
i=0
= h(0) + h(1)z^-1 + h(2)z^-2 + h(3)z^-3 + ... + h(n)z^-n
z の複素平面上の単位円 z= exp(jω) 上が周波数特性を表す。単位円の右端、(1,0) は、ω= 0 (DC)、
左端 (-1,0) は ω= 1 (fs/2) に対応する。また、DFT は、Z 変換の z に z= exp(2πk/N) という単位円の
N 分割点の値の計算である。
Z(z)= 0 をみたす z は、 n 次方程式 h(0)z^n+ h(1)z^(n-1) + ... + h(n)= 0 の解である。 一般に n 次方程式の解は、n 個の複素解をもつ。n 次の FIR フィルタは一般に n 個のゼロ点をもつ。
例として 11型 FIR フィルタは、z から 1+z^-1 に変換する写像であり、z=0 でない限りゼロ点は、 z+1= 0 の解、z= -1 である。121 型の解も重解で z= -1 になる。
N
y(i)= Σ h(k)x(i-k) (1)
k=0
で表されるフィルタである。i 時点の出力 y(i) が h(k) と k (k＞=0)だけ過去の入力 x(i-k)
との積和で作られるフィルタである。 FIRフィルタは、つねに安定である。
フィルタが不安定とは、 (1)出力値がどこまでも大きくなる(発散)、(2)入力がある時点か
らずっと 0 であるのに、いつまでも出力がある現象(発振)をさす。
FIR では y(i) がともに有限値である入力とインパルス応答との積和であるため有限値をとり、 発散せず、内部状態をもたないから発振をする可能性もない安定なフィルタである。 設計法には、窓による方法、周波数サンプリング、最適化(ミニマックス誤差)法があるという。
N
y(i)= h'(0)x(i) + Σ h'(k)y(i-k) (2)
k=1
現在の入力 x(i) と過去の出力値 y(i-k) (k＞0) とを使った定数の積和で出力 y(i) をつくるフィルタ。
過去の出力値を再利用するとは、それを内部状態としてもつことであり、無限の過去の影響を残す可能性
があり、安定保証はない。
FIR フィルタよりも一般に短い係数値で鋭いフィルタが可能といわれ、係数の個数の点で有利とされる。 しかし、係数の数値精度から当然それには繰り返し限界があり、IIR に可能なことは、FIR でもある程度 できるとされる。
IIRフィルタには、バターワース(最大平坦振幅特性)、ベッセル(最大平坦群遅延)、チェビシェフ(通過帯域 または阻止帯域の一方で等リップル)、楕円フィルタ(通過帯域と阻止帯域の両方で等リップル)などの過去の アナログフィルタの設計方法が使え、鋭いしゃ断特性のフィルタを作るには 1 桁以上も効率がよいといわれる。 つまり、昔風の優秀な特性のフィルタは、すべて IIR フィルタなのである。
ここで(1) と (2) の定義式に疑問の方は、 18. FIR と IIR の定義について を参照。
FIR の逆フィルタは、IIR である。また IIR の逆フィルタは、FIR である。FIR フィルタ (1)の x と y を入れ換えた式、
N
x(i)= Σ h(k)y(i-k)
k=0
は、 y(i) から x(i) を積和で作る式だが、これは逆フィルタの入出力関係式でしかない。これを因果的な系、
y(i)= (現在と過去の x(i)と 過去の y(i)だけを含む式) にできる。y(i) の項をΣから外に取り出し、
N
x(i)= h(0)y(i) + Σ h(k)y(i-k)
k=1
最初の 0 でない h(k) を h(0) としても一般性は変らないから、 h(0)!= 0 と仮定して、
N
y(i)= {x(i) -Σ h(k)y(i-k)}/h(0) (3)
k=1
これは、IIR フィルタである。IIR の式 (2)と比較して、
h'(0)= 1/h(0), h'(k)= -h(k)/h(0) (4)とする IIR フィルタが元の FIR フィルタの逆フィルタである。しかし、このように逆フィルタは、容易に 実現できると考えてはならない。例えば、h(0)= 0 である遅延を伴ったインパルス応答は、いくらでも存在 するが、その逆フィルタは、存在しない。時間遅れを引き戻すフィルタはないからである。h(0) が 0 でなくて IIR による逆フィルタ (3) が存在しても、それが発散、発振する可能性がある。ある周波数を完全にカット するフィルタの逆フィルタは、その周波数成分をわずかでも含んだ信号を入力するだけで発振すると考えられる。
FIRフィルタは、時間遅れを z^-1 としてその係数を z の多項式の係数として書き、=0 と置くと、 Σ h(k) z^-k = 0 となる。これは、通常の N 次方程式 h(0)z^N + h(1)z^(N-1) + h(2)z^(N-2) +....+ h(N)= 0 と解が等しく、N 次方程式は、多重解も別に数えるとつねに N 個存在する。解は FIR フィルタのゼロ点である。 (5)の y の項を左辺に、xの項を右辺に移項して Z 変換すると、Y(z)(1+Σh(k)/h(0)z^-k)= X(z)/h(0)、
N Y(z)/X(z)= 1/Σh(k)z^-k k= 0であり、FIR の正確に逆数の系であることが分かる。ゼロ点は、逆フィルタ IIR ではそのまま極になる。 一般に単位円の外にひとつでも極があれば発振する。
例として、11型のLPFは、z+1= 0 から、z= -1 にゼロをもつ。これは、値が毎回反転するサンプリング周波数 (fs)の半分 fs/2 の周波数を完全にカットするフィルタである。y(i)= x(i)+x(i-1) の逆は、y(i)= x(i)-y(i-1) である。インパルス応答は、y(i)= {1,-1,1,-1,...} という発振になる。
121型では、z^2+2z+1= 0 から z= -1 に重根がある。逆フィルタは、一般の波形には少しはその成分が含まれる fs/2 成分で発振する。y(i)= x(i)+2x(i-1)+x(i-2) の逆は、y(i)= x(i)-2y(i-1)-y(i-2)、インパルス応答は y(i)= {1,-2,3,-4,5,...} という線形拡大発振になる。
1-1型 HPF は、z-1=0 つまり z= 1 (DC)にゼロがあり、逆フィルタは、DC 成分が継続するだけで発散する。 これは逆フィルタが値を累積する和分になるからである。y(i)= x(i)-x(i-1) の逆は、y(i)= x(i)+y(i-1)である。 インパルス応答は、y(i)= {1,1,1,1,...}になる。
また、逆フィルタが安定であっても、逆フィルタは、その周波数のノイズに対して、敏感にノイズを拡大する フィルタとなって実用的には使えないことがある。そのような、基本的な難しさがあるために逆フィルタの様々な 実現方法が提案されてきたのである。Wiener Filter も逆フィルタの一つの解であり、Mourjopoulos の方法、 最小位相成分のみのスペクトル平坦化、(ゲイン等化) も、逆フィルタ実現のためのアルゴリズムである。
IIRが発散するかどうかは、IIRフィルタを実際に動作させてみればよい。しかしこれには、フィルタの長さ N と比較してはるかに長い動作をさせないと分からないことが問題である。
それ以外に判定する方法として、その極を求めて、全ての極のうち中心から最大に外れた極の振幅を調べれば よいが、N 次方程式の解を求めるのは計算量の点で有効とは思えない。
IIR の安定判別には、すでにいくつか解答がある。(1) Kharitonov の定理を使って4つの多項式のテストをする方法 (2)IIRフィルタを梯子型に変換し、反射係数の絶対値が 1 を超えないことを判定する方法がすでに使われている。
N-1
y(n) =Σ h(i) x(n-i)
i=0
N-1 N-1
Y(k) =Σ Σ h(i) x(n-i) w^{kn}
n=0 i=0
N-1 N-1
= Σ h(i) w^{ki} Σ x(n-i) w^{k(n-i)}
i=0 n=0
j= n-i とおいて、
i+N-1
= H(k) Σ x(j) w^{kj}
j=i
w が 1 の N 乗根だから、 w^{k(j mod N)}= w^{kj}
もし、x の系列 x_0〜x_{N-1}が N の周期関数と仮定するなら、
x(n mod N) = x(n)、上のΣは、j=0 から N-1 までの和と等しい。
N-1
= H(k) Σ x(n) w^{kn} = H(k) X(k)
n=0
N-1
φ_xy (k) = Σ x(i) y(i+k)
k=0
であり、そのフーリエ変換(周波数領域)は、X(f)Y(f)^* である。
信号 x 自身の相関を自己相関関数 (Auto Correlation) といい、 信号の統計的性質を抽出する重要な手段である。
N-1
φ_xx (k)= Σ x(i) x(i+k)
k=0
Σの上の N-1 は、周期関数でない場合、N-k-1とすることが多い。
自己相関関数は、
(1) k について偶関数であり、φ(-k) = φ(k)
(2) φ(0)は正又は0、φ(0)＞＝0、
(3) φ(0)が最大、φ(0)>= φ(k)、
(4) 最大の絶対値をもつ。φ(0)＞＝|φ(k)|
などの性質をもつ。
自己相関関数φ_xx (k)のフーリエ変換は、パワースペクトルΠ(f)である。 パワースペクトルは、 X(f)X(f)^*= |X(f)|^2 である。
Σ h(k) z^{-k} = 0 (6)
で記述され、 FIRフィルタのゼロ点に対応する。IIRフィルタでは、これは極点である。
その周波数の入力がなくても出力が存在する(発振)の可能性がある。
安定の判別は、多項式 =0 の根が、z = C exp(-j2πk/N) (Cは実数)なら、|C| >= 1 なら 発振し、|C| < 0 なら安定である。C は、1遅延時間における振幅の拡大、縮小率である。
あるフィルタの FIR 型のインパルス応答 h(i) を求めて、これの特性を修正し多項式 Σ h(k) z^-k = 0 の根を安定に変更する方法を求める必要がある。
線形位相フィルタ
h(n)= +- h(M-n)
H(w)= A(e^(jw))exp(-jwM/2)
n=M/2 で折り返す偶関数のインパルス応答は、 A(e^(jw)) が実数、奇関数のインパルス応答は A(e^(jw)) が純虚数。
対称的なインパルス応答をもつフィルタは、線形位相と言われる。k の中央を 0 にする対称なインパルス応答 を考えることが数学的単純化によく使われる。インパルス応答が対称的なフィルタは、全ての周波数応答の位相が 0 になる。周波数振幅特性は cos 関数の合成で作られる。入力を先に加算すれば積和量も半分ですむ。 位相特性をさしおいて周波数軸上の振幅特性を合わせることは、むしろ容易だからよく使われる。
一般の h(k)は、対称 h(k) と非対称 h(k) を含み、これを対称 h(k) に置き換えることはできない。 現実のインパルス応答が対称であることはほとんどあり得ないことを考えると、このフィルタがとくに重視される 正当な理由はない。逆フィルタの形成において、最小位相フィルタという方法があり、その最終的な出力は、 対称インパルス応答である。"最小位相 FIR だけが逆フィルタが安定な IIR になる。"というのは、誤った俗説 である。反例は、(0.5, 0.5) の逆数は発散するが、(0.51,0.5) の逆数が発散しないことである。
y(i)= r(0)x(i) + Σ b(k)y(i-k) Σk=1〜Nのy(i-1)項を分離、この式を適用。
= r(0)x(i)+ b(1){r(0)x(i-1)+Σ b(k-1)y(i-k)} + Σ b(k)y(i-k) 両Σk=2〜Nをまとめ、
= r(0)x(i)+ b(1)r(0)x(i-1) +Σ(b(1)b(k-1)+b(k))y(i-k) Σk=2〜Nのy(i-2)項を分離、この式を適用。
= r(0)x(i)+ b(1)r(0)x(i-1) +(b(1)^2+b(2)){r(0)x(i-2)+Σb(k-2)y(i-k)} + Σ(b(1)b(k-1)+b(k)) y(i-k) 両Σk=3〜をまとめ、
= r(0)x(i)+ b(1)r(0)x(i-1) +(b(1)^2+b(2))r(0)x(i-2) +Σ((b(1)^2+b(2))b(k-2)+b(1)b(k-1)+b(k)) y(i-k) Σのy(i-3)分離適用
= r(0)x(i)+ b(1)r(0)x(i-1) +(b(1)^2+b(2))r(0)x(i-2) +(b(1)^2+b(2))b(1)+b(1)b(2)+b(3)) y(i-3) + Σ k=4〜
= r(0)x(i)+ b(1)r(0)x(i-1) +(b(1)^2+b(2))r(0)x(i-2) +(b(1)^2+b(2))b(1)+b(1)b(2)+b(3)) r(0)x(i-3) + Σ k=4〜
r(1) 以降の最初の 3 係数を示すと、
r(1)= b(1)r(0) r(2)= (b(1)^2+b(2))r(0) r(3)= (b(1)^2+b(2))b(1)+b(1)b(2)+b(3))r(0)逆フィルタでは、r(0)= 1/h(0), b(k)= -h(k)/h(0) であるから、
r(0)= 1/h(0) r(1)= - h(1)/h(0)^2 r(2)= - h(2)/h(0)^2 - h(1)^2/h(0)^3 r(3)= - h(3)/h(0)^2 + h(1)h(2)/h(0)^3 - h(1)^3/h(0)^4 + h(2)h(1)/h(0)^5つまり、逆フィルタは、y(i)= 1/h(0) x(i) - h(1)/h(0)^2 x(i-1) - (h(0)h(2)-h(1)^2)/h(0)^3 x(i-2) +... と書ける。
こうして、逆フィルタを FIR 実現するための逆数のインパルス応答 r(k) を求めることができるが、このような 繁雑な展開は実は不要である。r(k)は、IIRのインパルス応答であるから、次のように記述できる。
r(1)= b(1)r(0)
r(2)= b(2)r(0) + b(1)r(1)
r(3)= b(3)r(0) + b(2)r(1) + b(1)r(2)
r(4)= b(4)r(0) + b(3)r(1) + b(2)r(2) + b(1)r(3)
i
r(i)=Σ b(k) r(i-k) = -Σh(k)/h(0) r(i-k)
k=1
これは、IIRの動作自体である。IIR 逆フィルタの式に戻れば、y(i) は、kだけ過去の出力 y(i+k) からの影響
h(k)y(i+k)を x(i) から引算するものだった。
N h(0)y(i)= x(i)- Σ h(k)y(i-k) (5) k=1これは、次のようにも記述できる。i を i+k に置き換えΣ中の処理を分散させ y(i) を作ったとき k (0＜k＜N) だけ未来の x(i+k) から y(i) の影響 h(k)y(i) を引算する処理となる。
h(0) y(i+k)= x(i+k) - h(k)y(i) (for 0 < k < N) (5'')この因果的でない (5'') は、(5)と同じ結果をもつ。C 言語記述を次に示す。x[] はインパルス(1,0,0,0,..) に初期化する。処理後、配列 x[] が FIR の逆フィルタの FIR 記述となる。
float x[N]= {1};
int i,k;
for(i=0;i<N;i++){
y[i]= x[i]/h[0];
for(k=1;k<N;k++) x[i+k]-= y[i]*h[k];
}
多項式展開によって逆インパルス応答 r(k) を求める処理は、IIRフィルタのインパルス応答を求めることと合同 であり、これは、単に形式的に FIR フィルタにしているだけである。 形式的に FIR にすれば安定な系になった と考えることはできない。逆フィルタが発散する場合、それを FIR 型にするための逆インパルス応答係数 r(k) が発散するのである。これは、信号とシステムとが対等であることを示している。
室内音響のインパルス応答を使っての室内残響の改良、ぼやけた画像の鮮明化などの逆フィルタの実現には、 順特性のインパルス応答を測定し、h(k) とすると、
(1) それを使って、式 (3) で IIRフィルタ動作をさせることである。基本的に h'(0)= 1/h(0), h'(k)= -h(k)/h(0) という IIR フィルタ動作が必要であり、中心インパルス値は 1、周辺がぼやけ インパルス応答の負のインパルス応答 -h(k)を与える。全体的大きさには、1/h(0) が掛かる。 しかしそれを発散しない IIR 動作させる必要がある。
(2) そこにおいてインパルス応答の修正を必要な場合行うことである。IIRフィルタとしての安定性を 確保する限りにおいて、逆フィルタは動作しうる。発散をしないだけでなく、発散を止める限界から 多少離れたほうが、多項式展開で出て来る値の絶対値も巨大にならなくて済む。これは、発散ぎりぎり の場合にリンギング現象がでるが、そのような歪み、アーティファクトを防止するためである。
発散するIIRを発散しないものに修正するには、いくつかの経験的な方法があって、実用的にそれら は改良されていくものである。インパルス応答の修正は微小量ですむものである必要がある。 インパルス応答の修正が微小であれば、逆フィルタとしての特性劣化も微小であるからである。 それらの代表的な修正方法には、2 種類あり、両者を併用してもよい。
(1) h(0) の振幅をわずかに(0.数%)増加させる。
1 次方程式 a x + b = 0 は、a の増大によって、解 x= -b/a をわずかに 0 に近付ける。
2 次方程式 ax^2 +bx + c= 0 では、a の増大は、対称軸 -b/2a を 0 に近付け、
判別式 D をわずかに c の符号によって移動させる。x= (-b +- √(b^2- 4ac))/ 2a
しかし、そのような修正が高次において一般的に IIR を安定化する方向であるかどうかについては、
初等的に明確とはいえない。h(0) を大きくするのは、順方向の FIR フィルタを入力との線形加算であり、
薄められた FIR の逆フィルタの効果になる。そのとき応答周波数の変化があるかどうかも明確でない。
(2) 相似変換の利用 (h(k) に exp(-k*0.01)を掛ける。)
例えば、h(k) に exp(-k*0.01)を掛けることによって z の複素平面上の関数を相似的に縮小する。
Σ h(k) z^-k = 0 を複素数 z が満足するとき、h'(k) = h(k) exp(-ak) とする、Σ h'(k) z^-k= 0 を満足する z は、 z'= z exp(-a) である。なぜなら、Σ h'(k) z'^-k にh'(k)とz'とを代入すれば、 =Σ h(k)exp(-ak)(zexp(-a))^-k = Σh(k)z^-k= 0。
つまり、k>=0 について h(k) を exp(-ak) 倍に指数縮小することは、すべてのゼロ点を原点に向かって 比例的に exp(-a) 倍に縮小する相似変換である。相似変換のよい所は、ゼロと極の周波数を保存して、 そのディッピング、ピーキングの鋭さを減少するだけであることである。
ゼロ点が単位円上とすこし外にある場合、この方法によって単位円の内側に移動することができる。 大きく単位円から外れたゼロ点は、h(k) に exp(-ak) を掛ける方法では、インパルス応答を短縮する 特性変化が大きく、それがこの方法の欠点となる。現実に測定した室内音響のインパルス応答からは、 単位円から数%以上外にでるゼロ点は、極くまれであるから、この方法だけで実用可能かもしれない。
複素数 z をゼロ点にもつインパルス応答は、その逆フィルタの入力に与えなくても z 出力をもつ、 z が自然発生するということである。z が単位円上つまり、三角関数波形であるなら、定常な発振の 発生であり、z が単位円の外側なら、拡大波形の自然発生である。波形拡大はいつまでも続かないか ら問題ないかというと、そうではなく大振幅時の非線形性によって定常波形に変る。z が単位円の少し 内側なら、それは三角関数の減衰波形であり、自然発生しても発振ではない。
例えば、11型の FIR (0.5,0.5) の多項式展開は、(+2,-2) の繰り返しであるが、(0.51,0.5) に修正 すると、時定数 50 の減衰振動になる。(1.96, -1.922,..)
また、121 型の FIR (0.25,0.5,0.25)の多項式展開は、(4., -8., 12., -16., 20. ...) という振幅が 直線的に拡大していく振動であるが、h(0)を 0.26 に修正すると、やはり時定数 50の減衰が掛かった ビート振動波形になる。11 型と 121 型の発散は、0次、1次の発散であり、(1) の h(0) のわずかな 増加で安定化できる。
1331型の逆 IIR は 2次の発散をするが、相似変換で安定化できる。指数的収縮は、高次の多項式的発散 をも押えることができると考えられる。
#include <stdio.h>
#include <math.h>
#define N 256
static float x[N]= { .5, .5, 0, 0, };
static float y[N]; // inversed impulse response
static float z[N];
static float c[N]= {1.,0,0};
static float Xr[N];
static float Xi[N];
static void
dft(void)
{
int i,j;
for(j=0;j<N/2;j++){
Xr[j]= Xi[j]= 0;
for(i=0;i<N;i++){
double t= -2*M_PI*i*j/N;
Xr[j]+= x[i]*cos(t);
Xi[j]+= x[i]*sin(t);
}
z[j]= sqrt(Xr[j]*Xr[j]+Xi[j]*Xi[j]);
}
for(j=0;j<=N/2;j++) printf("%f\n",z[j]);
}
static void
pol_inv(void)
{
int i,j;
double r;
// polynominal inversion
for(j=0;j<N;j++){
y[j]= r= c[j]/x[0];
for(i=j;i<N;i++) c[i]-= r*x[i-j];
}
// filter to test XY= 1
for(j=0;j<<N;j++) {
double t= 0;
for(i=0;i<10;i++) t+= x[i]*y[j-i];
z[j]= t;
}
for(j=0;j<N;j++) printf("%f\n",z[j]);
}
int
main()
{
pol_inv();
//dft();
}
0 を中心にする対称型 FIR フィルタ h(i) の周波数応答は、実部のみだから H(k)= |X_r(k)|= Σh(i)cos(2πik/N)。周波数応答は、緩やかなコサイン関数数個でできている。ここで DFTの N は、 k/N の分解能を与えているだけで、フィルタの特性ではない。
11 型は、0 を中心に i= +-1/2 で値 1/2 をもつから H(f)= cos(πf/N)。1/2 に 1 サンプル離れて1/2 があるからとすると、H(f)= 1/2 √((1+cos(2πf/N))^2+sin(2πf/N)^2) = √((1+cos(2πf/N))/2) = cos(πf/N)
フィルタの縦続接続は、フィルタ係数の簡単な計算でできるが、フィルタの並列接続は、 さらに簡単である。121 型のLPFを 1 から引くと、121の中央の山を 1 に合わせて、1 から引くと、-1/4,1/2,-1/4 となり、例 5 の符号反転になる。 周波数応答特性は、線形であり、フィルタの係数の加減算も周波数応答上の加減算になる。
DCゲインを 1 にする(a, 1-2a, a)では、上の b'= 0.5 - a であるから、a < 0.5 では軸は負(LPF)、 a>0.5 で軸は正(HPF)。D'= b'^2 - a^2= (0.5-a)^2-a^2= 0.25 -a >0で左右、<0 で上下である。
線形性がその理由かもしれない。sin も cos も同じように扱う必要があるとき、それらをひとつに表現できるのは、 複素数である。cos(x) = 1/2 (exp(jx) + exp(-jx)), sin(x)= j/2 (exp(jx) - exp(-jx)) 線形性から、sin も cos も同じ複素数の定数倍と和で表現できる。それなら、それを構成するもっとも単純な要素で表現する、それが e の 虚数乗である。exp(jw) を使わずに、sin と cos の 2 成分で表現すればいいのではないか。その通りだ。それが オイラーの exp(jw) の成分表示 exp(jw)= cos(w) + j sin(w) である。sin と cos とが同様に系を通ることが exp(jw) の意味である。しかし、1 成分の信号を表すのに、どうして複素数に拡張する必要があったのだろうか。
実関数をフーリエ級数展開するとき sin と cos のふたつの展開を使う必要があった。任意の実関数は、偶関数と 奇関数の和であり、それぞれ cos(nwt) と sin(nwt) の級数で表現される。同じ周波数 n の sin と cos の対は、 振幅と位相をもつ複素数である。実信号の波形が周波数軸上の複素数表現を要求しただけである、と言えるか。
フーリエ変換では原信号も複素信号である。その線形性から虚数jを掛けることもスカラー倍とほぼ同じことだし、 入力の和は出力の和となるわけだから、信号を複素にするのも何も問題ではない。しかしそれは、"オッカムの剃刀" から無駄な2次元への拡張かと疑われる。複素信号は、ふたつの独立した実信号を扱うのと同じく冗長である。 独立 2 系統信号を扱う系、例えばステレオ音響信号は、複素信号として扱うことが適当だが、実信号系に複素信号 を使うと処理を倍化する。逆に、複素 FFT は、独立二実信号の FFT 処理を実際にすることができるのである。
同時に二信号を扱うことが物事をシンプルにするはずがないという疑問に対して、いや複素表現が物事をシンプル にするのだ、信号の本質は複素数だということかもしれない。信号自身は実信号でも、その信号の意味は違うかも しれない。動力学(ダイナミックス)から来る、位相空間概念 (位置と速度の 2次元表示)には、複素数が最適である。
時間的な変位が周波数軸上で exp(-iwt) という複素回転であることには、何か物理的意味があるだろうか。
時間の関数から周波数の関数へという変換において sin と cos が取られ、それらがペアになって扱われることに、 深い理由があるのであろう。両者が別の任意の関数でなく、両者だからこそ exp(iw) を使うことも意味がある。 しかし、もとの信号を偶関数に制限すれば、cos(nw) だけですむわけで、sin (nw) と cos(nw) とのように互いに 独立であるなら、別の信号系列に変換した表現もあり得るのである。
そこでの考え方は、"画像は、すべて劣化している。復元するためのヒントは、画像自身にある。" であろう。 点光源の例として代表的な星の像は、画像処理のすべての過程を通した広がり関数を表している。広がり関数は FIR のインパルス応答そのものであり、それを使って多少なりとも尖鋭な画像を得ようとする努力は、当然の作業 である。画像を得るまでの労力の大きかった宇宙望遠鏡画像などは、そのような画像尖鋭化の処理を経るものであろう。
画像は、音声、音響とはかなり違う。音声を 1 次元というのは、時間軸上の音圧又は速度が波形になっているのであり、 画像が 2 次元というのと整合しない。その同じ言い方なら 2 次元画像が、時間軸上で変化していくのであるから 画像は 3 次元である。というのは、画像において、IIRフィルタというのは、何であろうかという疑問が湧くから である。時間軸上の再現を再利用するのは、画像符号化の"予測"であり、2 次元画像内 IIR は、フレーム内予測 であろうか。FIRフィルタの逆フィルタを直接行うのに FIR の逆が IIR であるというのは、画像においてなにか。 それは、かなり視点が変わるが、神経生理学の"側抑制"(Lateral Inhibition)ではないだろうか。神経結合回路の 用語が実は、画像復元の本当のフィルタである。FIR フィルタの逆を行うのは IIRフィルタである、それは正しい。 それは、画像において、近傍の画素出力を抑制入力するものだろう。y(n)= {x(n) - Σ h(k)y(n-k)}/h(0) という式に相当する形は、近傍画素(入力でなく)出力から側抑制することなのである。
y(m,n)= {x(m,n) - Σh(j,k)y(m-j,n-k)}/h(0,0)
この式が一意に値を与え、安定であるという保証はない。これは、一種の無理難題かもしれない。 視覚神経系において中心には、正の反応画素があり、それを取り囲む円形の負の反応をする抑制領域をもつことは、 以前から有名なことである。ここで注意すべきは、考え方の容易な FIR で FIR フィルタの逆フィルタを行おうと しがちなことである。 FIR を FIR で打ち消すことは、容易な手法ではあるが基本的に膨大化を招くものである。 側抑制層は、その出力を互いに近傍の細胞に抑制入力化している。神経回路網は、そのことを教え、自然はそのこと をすでに熟知し駆使しているのかと思われる。画像自身の情報による画像の尖鋭化を恐らく仕組として生物は、行っ ているのである。そこで我々の技術がまだ採り入れていない、必要な機構は、自動的なホワイト化、白化である。 つねに、必要な情報を得るためには、自動的な逆フィルタ処理がなされる必要があるからである。
(1)まず画像の解像度を形式的に上げる。例えば画素をLPFを通して分解能を2倍(画素数を4倍)にする。 このとき LPF は、もとの画素分解能のもっている周波数である fs/4 までをできるだけ通し、それ以上の周波数を 遮断する LPF であろう。元の画素の存在する位置の画素値はそのままにして、内挿画素だけ偶数タップの全域通過 フィルタ簡単には (-1,5,5,-1) のようなフィルタを使って作成する。画像拡大に使われるフィルタを考えればよい。
(2)画像のある局所領域、例えばブロックが他の複数の画像のどこに近似するかを調べる。これは、動き推定と同様 な処理である。MPEG の半画素までのサーチと話が類似するが、先に画素を増やしているので、これは整数画素まで のサーチである。この動き推定は、符号化用でないので、ブロックサイズも任意であり、並行移動だけでなく、 画像を傾けるとか拡大縮小を伴うアフィン変換も許されるし、その変換パラメータは、捨て去るものであるから、 複雑でもよい。しかし、どこまでも任意というわけではない。自然物体の動き、変形、及び画像化に伴う変形を 除去する限りのものであろう。
(3)複数の画像の類似するブロックを平均する。このとき、平均化が解像度を失わないように画像の動き推定過程を 正しく行っている必要がある。(2)の段階の誤りはすべて解像度劣化につながるし、画像のもとの物体の変形など による画像の変化も大きいとこの処理自体が画像劣化となり得る。この平均処理も多くの工夫の入る余地がある。 通常の非線形フィルタのように、平均でなく、メディアン(中央値)を取ることも考えられる。
同じ物体の同じ位置に対応する画素値を取り出し平均することは、雑音を除去し微小な画像の特徴を明確にするだけ でなく、もとの解像度以上の解像度を引き出すことができるのである。超解像は、それ自身、理にかなっていて、 動画像の複数の画像が1枚の画像よりも情報量をもっている理由の一つの側面を説明するものである。動画像が画像の 時間的変動のためだけに使われているのだろうか。パニングは画角の拡大のためであり、もう一つの超解像に対応する ものであろう。そうでない同じ視野角をとった動画像は、ほとんどが無駄となるが、そういうことはあるはずがない。 それは、分解能を上げるために存在するのではないか。神経生理的な確認はないと思うが、ヒトの目は恐らく超解像 もやっているのではないか。
すでに超解像の工学的な実験は沢山行われて、効果が確認、報告されている。よく分解された動画像中の画像4枚 を使えば恐らく、分解能は2倍程度まで上げることができる。画素の位置のゆれが画像毎の画素位置を揺らし、 固定したカメラ以上の情報を与える。つまり、カメラは固定しないほうがよい。それ自身ゆっくり動く物体の画像は、 固定したカメラでも超解像できる、という逆説的な説明になる。 動物体領域の超解像は、それ自身難しいとされているのであるが。
(3) の平均化で得られる画像のブロックはサーチのもとにすると、もとの入力画像のブロックよりもサーチにより 適したものといえる。すでにノイズ除去されているからである。
M y(n)= Σ b(r) x(n-r) r= 0それに対して、IIR システムは、FIR を包含した混合型で示す。
N M y(n) - Σa(n)y(n-k) = Σb(n)x(n-r) k=1 r= 0こう書いて便利なのは、両辺を Z 変換して、
Y(z) (1 - Σa(k)z^-k)= X(z) Σb(r)z^-rその伝達関数 H(z)= Y(z)/X(z) の分子が FIR の式、分母が 1 - Σの形になる。
M N H(z)= Y(z)/X(z)= Σb(r)z^-r/(1 - Σa(k)z^-k ) r=0 k=1もし M＜N なら、H(z) は部分分数の和で表せ、
N
H(z)= Σ{ A_k/(1- d_k z^-1) }
k= 1
(係数は、A_m= H(z) (1- d_m z^-1)|z=d_m (m= 1,2,..,N) から求まる。)
因果的システムでは、
N h(n)= Σ A_k(d_k)^n u(n) k= 1(これには、A_k(d_k)^n u(n) の z 変換が A_k/(1 - d_k z^-1) (|z|＞|d_k|)であることを使う。) ゆえに、 h(n) は無限に長い持続時間をもつ。と、説明する。
しかし、この IIR の定義が不便であるのは、FIR と IIR の逆フィルタの明白な関係が式から見えないことである。 定義式の分子の FIR 部分を b(r)= 1,(r==0),0 (r>0) として、分子の Z 変換を 1 にしても、分子と同様のΣの 形の式を分母に作ることを困難に思わせる。もちろん、本当はそうではなく、b(0)= 1/h(0)、b(r)= 0 (r>0) として a(k)= -h(k)/c(0) とすると次式になる。
N Y(z)/X(z)= 1 / Σh(k)z^-k k=1FIR の逆フィルタの IIR への直接対応をいうと、私に IIR フィルタの定義の式をお教え下さる方々がおられた。 その方々が指導的な立場であるとき、そのことがどのようなことであるかを御存知ないのである。 このような不用意な御意見に対して常に耐えることが私には必要なことであった。IIR の混合型定義は、 理解に苦しむものであり、それによる損失は、大であると思う。