🐙

FIRシステムにおけるリサンプリングの基礎づけ

2020/11/20に公開

Python 3

はじめに

画像処理でも音声処理でもリサンプリングはたまに必要になります。
本記事では、リサンプリングで何が起きているのかをFIRシステムの観点から基礎づけします。
あくまでFIRシステムにおける基礎づけなので万能では無いですし見えている世界は狭いです。
でも、理論が何も無いよりはマシなので、リサンプリングがよくわかって無いという人には本記事を入門とする事をお勧めします。
なお、本記事はデジタル信号に限っています。

FIRシステムって何？

有限インパルス応答となる線形時不変システムの事です。場合によっては因果的である事を要請する場合もあります。

システムとは、入力に対して何か出力するものです。
線型性とは掛け算と足し算の分配法則や結合法則が成り立つ事です。
時不変とは、時間（空間）が移り変わっても仕組みは変わらない（不変）という事です。
因果的とは、現在の出力値は現在の入力値および過去の入出力値によってのみ決まるという制約です。

これらの条件（制約）を満たすシステムを線型時不変システムと呼びますが、FIRシステムはその内のインパルス信号を入力すると出力信号（インパルス応答）が有限長で打ち切られるようなシステムの事です。

一言で言うと入力値列に対する畳み込み和（Convolution）で実現できる処理全てがFIRシステムです。
$y$ を出力値、 $x_i$ を入力値列、 $a_i$ を係数列として式で書くと、FIRシステムは以下の通りです。

y=\sum_{i=0}^{N-1}a_{i}x_{i}

係数列と入力値列をそれぞれベクトルと見なせば内積を計算しているだけです。入力値列に対してスライドしながら内積を計算していくのが直線畳み込み和です。

巡回畳み込み和と周波数スペクトル

直線畳み込み和だとちょっと困るので、係数列を巡回化して巡回畳み込み和と考えます。
巡回行列化とは、例えば $N=4$ で係数列 $a_0,a_1,a_2,a_3$ とした場合、以下のような行列を構成する事です。

C = \begin{pmatrix}\begin{array}{cccc} a_0 & a_1 & a_2 & a_3 \\\\ a_3 & a_0 & a_1 & a_2 \\\\ a_2 & a_3 & a_0 & a_1 \\\\ a_1 & a_2 & a_3 & a_0 \end{array} \end{pmatrix}

そうすると、巡回畳み込み和の計算は巡回行列による行列積と見なせます。
係数列を巡回化して構成した巡回行列を $C$ とすると、入力値列ベクトル $\vec{x}$ と出力値列ベクトル $\vec{y}$ の関係は以下の通りに表せます。

\vec{y}=C\vec{x}

$C$ をかける事はそのまま巡回畳み込み和を表していますが、その性質がどのようなものかを調べるために、固有値分解を行います。
正方行列で表現可能な正則（ここでは逆行列を持つ）行列は、その変換を回転、係数倍、逆回転の行列積に分解できます。この時の係数倍が固有値で、回転と逆回転が固有行列とその逆行列になります。

巡回行列 $C$ の固有行列を $P,P^{-1}$ 、固有値が対角成分に並んだ固有値行列を $\Lambda$ とすると以下の関係になります。

C = P \Lambda P^{-1}

両辺の左から $P^{-1}$ 、右から $P$ をかけると、巡回行列 $C$ を対角行列 $\Lambda$ 化できます。この操作を行列の対角化と呼びます。 $PP^{-1}=P^{-1}P=I$ （ $I$ は単位行列で対角要素に $1$ が並びその他の要素は $0$ ）です。

P^{-1}CP = \Lambda

巡回行列の固有行列を求めると、離散フーリエ変換（Discrete Fourier Transform:DFT）行列を導出できます。
よって、巡回行列はDFTで対角化でき、対角化で現れる成分がスペクトルです。

$P^{-1}=M_{DFT}$ 、 $P=M_{IDFT}$ として、 $\vec{y}=C\vec{x}$ の式の巡回行列 $C$ を固有値分解した形で表現してみます。

\begin{array}{ccl} \vec{y} & = & C\vec{x} \\\\ \vec{y} & = & M_{IDFT} \Lambda M_{DFT} \vec{x} \end{array}

そして両辺にDFTを施します。単に行列 $M_{DFT}$ をかけるだけです。

\begin{array}{ccl} M_{DFT} \vec{y} & = & M_{DFT} M_{IDFT} C M_{DFT} \vec{x} \\\\ M_{DFT} \vec{y} & = & \Lambda M_{DFT} \vec{x} \end{array}

$\Lambda$ は対角行列なので、ベクトル $M_{DFT}\vec{x}$ との行列積の結果は、各対角要素をベクトルの各要素と掛け算する結果になります。つまり、 $N=4$ の時、固有値を $\lambda_0,\lambda_1,\lambda_2,\lambda_3$ とし、 $\vec{X} = M_{DFT} \vec{x}$ 、 $\vec{Y} = M_{DFT} \vec{y}$ とすると、以下の通りの計算です。

\begin{array}{ccl} \begin{pmatrix} Y_0\\\\ Y_1\\\\ Y_2\\\\ Y_3 \end{pmatrix} &=& \begin{pmatrix} \begin{array}{cccc} \lambda_0 & 0 & 0 & 0 \\\\ 0 & \lambda_1 & 0 & 0 \\\\ 0 & 0 & \lambda_2 & 0 \\\\ 0 & 0 & 0 & \lambda_3 \end{array} \end{pmatrix} \begin{pmatrix} X_0\\\\ X_1\\\\ X_2\\\\ X_3 \end{pmatrix}\\\\ \begin{pmatrix} Y_0\\\\ Y_1\\\\ Y_2\\\\ Y_3 \end{pmatrix} &=& \begin{pmatrix} \lambda_0 X_0\\\\ \lambda_1 X_1\\\\ \lambda_2 X_2\\\\ \lambda_3 X_3 \end{pmatrix} \end{array}

この通り、巡回畳み込み和をDFTでスペクトル領域に変換すると、入力値ベクトルのスペクトルを、巡回畳み込み和のスペクトルと入力値ベクトルのスペクトルとの要素積（アダマール積）で表現できます。
これがDFTにおける畳み込みの定理です。
ちなみに、逆に元の空間でのベクトル同士の要素積は、DFTのスペクトル領域ではスペクトル同士の巡回畳み込み和になります。

さて、DFT行列は1のN乗根（複素数）で構成されており、N乗すると必ず1に戻ってきます。これは $\sin$ 関数と $\cos$ 関数で表現でき、周期性のある波と見なせるので、DFTによるスペクトルを周波数スペクトルと呼びます。

窓関数と周波数スペクトル

窓関数を単にDFTで周波数スペクトルに変換すれば、周波数スペクトル領域でどのようなスペクトルを巡回畳み込みする事になるかわかります。

import numpy
from scipy.fftpack import fft, ifft
import scipy.signal as sig
import sys
import matplotlib.pyplot as plt

T = 256
wf = sig.hann(T)
amp = numpy.roll(numpy.abs(fft(wf) / T), T//2)
amp_db = 20 * numpy.log10(amp + sys.float_info.min)
plt.plot(amp_db)
plt.ylim(-150, 0)
plt.show()
plt.close()

コードでは振幅値を見るために絶対値を計算していますが、実際に畳み込まれるのは複素スペクトルのままである事に注意してください。

リサンプリング

ダウンサンプリング

ダウンサンプリングは間引き処理ですが、その時に何が起きているでしょうか。
ダウンサンプリングの代わりに、 $1$ と $0$ を交互に繰り返すくし状の窓関数をかけると考えて、その周波数スペクトルを見てみます。

import numpy
from scipy.fftpack import fft, ifft
import scipy.signal as sig
import sys
import matplotlib.pyplot as plt

T = 256
wf = numpy.zeros(T)
wf[::2] = 1
amp = numpy.roll(numpy.abs(fft(wf) / T), T//2)
amp_db = 20 * numpy.log10(amp + sys.float_info.min)
plt.plot(amp_db)
plt.ylim(-150, 0)
plt.show()
plt.close()

くし状窓関数のスペクトルもくし状になりました。単純に半分に間引く場合はスペクトル領域では最大整数周波数の半分の間隔でくしが立ちます。 $\frac{1}{\Delta T}$ に間引く場合はスペクトル領域では $\frac{T}/{\Delta T}$ 間隔でくしが立ちます。

ダウンサンプルするということは、このようにスペクトル領域でくし状スペクトルの巡回畳み込み和を行うことになります。そして、 $\frac{1}{2}$ のダウンサンプルの際は、対称成分含めてくし状の巡回畳み込み和が計算されることになるので、最大正規化整数周波数を $F=\frac{T}{2}$ とすると、 $\frac{F}{2}$ 以上の成分を入力値スペクトルが持っていた場合、成分がオーバーラップして畳み込まれます。

以下はオーバーラップが発生しない範囲までしか原信号が周波数成分を持たない場合の例です。

import numpy
from scipy.fftpack import fft, ifft
import scipy.signal as sig
import sys
import matplotlib.pyplot as plt

def __make_cos_wave(f, t):
    return numpy.cos(2 * numpy.pi * f * t)

def __to_amp_spectrum_db(values):
    amp = numpy.abs(fft(values) / len(values))
    return 20 * numpy.log10(amp + sys.float_info.min)

T = 256
wf = numpy.zeros(T)
wf[::2] = 1

signal = numpy.zeros(T)
t = numpy.linspace(0, 1, T, endpoint=False)
for f in range(T//4):
    # 振幅値を減衰しながらミックス
    signal += __make_cos_wave(f, t) / (10 * (f+1))

signal_amp_db = __to_amp_spectrum_db(signal)
downsampled_amp_db = __to_amp_spectrum_db(signal * wf)
plt.plot(signal_amp_db, label="plain signal spectrum")
plt.plot(downsampled_amp_db, label="downsampled spectrum", linestyle="--")
plt.ylim(-150, 0)
plt.legend()
plt.show()
plt.close()

以下はオーバーラップが発生する例です。生成する $\cos$ 波形の周波数範囲を $\frac{F}{2}$ 以上に広げています。

import numpy
from scipy.fftpack import fft, ifft
import scipy.signal as sig
import sys
import matplotlib.pyplot as plt

def __make_cos_wave(f, t):
    return numpy.cos(2 * numpy.pi * f * t)

def __to_amp_spectrum_db(values):
    amp = numpy.abs(fft(values) / len(values))
    return 20 * numpy.log10(amp + sys.float_info.min)

T = 256
wf = numpy.zeros(T)
wf[::2] = 1

signal = numpy.zeros(T)
t = numpy.linspace(0, 1, T, endpoint=False)
for f in range((T*3)//8): # 範囲を変更
    # 振幅値を減衰しながらミックス
    signal += __make_cos_wave(f, t) / (10 * (f+1))

signal_amp_db = __to_amp_spectrum_db(signal)
downsampled_amp_db = __to_amp_spectrum_db(signal * wf)
plt.plot(signal_amp_db, label="plain signal spectrum")
plt.plot(downsampled_amp_db, label="downsampled spectrum", linestyle="--")
plt.ylim(-150, 0)
plt.legend()
plt.show()
plt.close()

周波数 $\frac{F}{2}$ をはみ出した成分がオーバーラップして畳み込まれています。これがダウンサンプルにおける高域成分のオーバーラップによるエイリアシングです。

ちなみに、入力信号の周波数成分を周波数 $\frac{F}{2}$ 以下に制限できていれば、くし状スペクトルによる巡回畳み込み和で生成される繰り返し成分は、元になる周波数成分と完全に一致します。この状態で間引くとサンプリングの定理によりエイリアシングが発生しますが、繰り返し成分の外形が完全に元の周波数成分と一致するので、アップサンプリングにより完全に復元可能です。

要はダウンサンプリング時にローパスフィルタで高域成分をカットしておかないと、繰り返し成分のオーバーラップにより思わぬアーティフィシャルノイズ（エイリアシングノイズ）が発生してしまうということです。

アップサンプル

アップサンプルの場合はダウンサンプルの逆を考えます。つまり、アップサンプルしたい倍率になるように入力信号に対してゼロ値をくし状に挿入します。
ゼロ値挿入した入力信号は、ダウンサンプル時のくし状窓関数をかけた状態と同じ状態になります。つまり、アップサンプル前の入力信号の最大正規化整数周波数 $F$ までの成分が繰り返し出現する信号になります。
アップサンプル後の最大正規化整数周波数は $2F$ になるので、 $F$ 以上の成分をカットするローパスフィルタをほどこせば、綺麗にアップサンプリングできます。

ここまでの説明でわかる通り、FIRシステムによって基礎付けしたリサンプリングにおいて重要なのは、ローパスフィルタの性能です。理想的には通過域を完全フラットな状態で通過し、阻止域は完全にゼロに落とすようなローパスフィルタですが、残念ながらそのようなフィルタは有限長のフィルタカーネルでは実現できません。本記事では出てきませんが、IIRフィルタでも実現できません。

また、入力信号が元々ダウンサンプリングによって生成された信号であり、それをアップサンプリングで復元しようと考えた場合、失われた高域成分をどのように補償するかが超解像の研究テーマです。ここまでのFIRシステムに基礎付けられたリサンプリングでは、必ず高域成分を捨てる必要があり、それを復元する手立てはありません。

非整数倍のリサンプリング

例えば画像を $1.5$ 倍にアップサンプリングしたい場合にどうするかというと、これは有理数比なので、 $3$ 倍してから $\frac{1}{2}$ にすると考えます。そうすると、整数倍アップサンプルと整数倍ダウンサンプルの組み合わせで実現できます。

計算機上で扱うリサンプリングは必ず有理数比になるはずなので、整数比ダウンサンプリングと整数比アップサンプリングを実装し、その組み合わせで必ず実現できます。（単純な実装だとメモリを大量に消費する可能性があります。）

おわりに

ひとまず理論体系を簡単に、かつ知らなければならない事を極力排除して説明しました。グラフはスクリプトをPython対話モードでコピペ実行すれば各自の環境で閲覧可能なので載せていません。

本記事はあくまでFIRシステムによる基礎付けです。おそらく音声処理の場合は本記事に書いてある事が大分役に立つと思いますが、画像の場合はそもそも非定常非周期な信号であることがほとんどなので、周波数という概念に限らず色々と工夫が必要になると思います。ただし、モアレなんかは本記事の内容である程度説明できたりします。画像処理系の方は本記事の内容を保守的な手法として覚えておくと良いかもしれません。

そもそもなんでこんな古典的な内容の記事を今更書いたかというと、技術書典7で頒布したデジタル信号処理の本に対するフィードバックが全然無いからです。
何か疑問点やもう少し掘り下げて説明してほしい箇所があれば、私の能力の範囲で対応します。

以上です。

はじめに

FIRシステムって何？

巡回畳み込み和と周波数スペクトル

窓関数と周波数スペクトル

リサンプリング

ダウンサンプリング

アップサンプル

非整数倍のリサンプリング

おわりに

Discussion