インターネットに解説があんまり無くて永遠に理解できてなかったけどやっとわかってきた (本を探せばあったのかな……)　 $z$ 変換をちゃんと理解してなかった

参考にしたページ
https://www.rle.mit.edu/dspg/documents/Proc_Frequency_1999.pdf

$z$ 変換

インパルス応答とかの離散信号をフーリエ変換すると正規化角周波数 $\omega$ (単位が rad / sample のやつ) の式 $F(\omega)$ が得られる。 $\omega$ は $-\infty$ から $\infty$ までの数直線上の値 ( $2\pi$ ごとにループするけど) である。数直線上から自由に角周波数 $\omega$ を選んで $F(\omega)$ を計算すると、元の信号がインパルス応答ならその角周波数でのフィルタの振幅特性や位相特性がわかる。

これを、何を思ったか突然 $z=\cos\omega + j\sin\omega$ で変数変換してみたものが $z$ 変換らしい。このとき $\omega$ と $z$ の対応は下のようになる。

つまり、 $\omega$ が数直線上を動くのは、 $z$ が単位円上をぐるぐるするのと同じ。 $\omega$ は数直線上しか動かないので $z$ は必ず単位円上の値だと考えていい。(本当はだめかも？とりあえず今は大丈夫) ただの変数変換なので、 $z$ 領域の式 $X(z)$ で、 $z$ に $j$ を代入することで正規化角周波数 $\omega=\frac{\pi}{2}$ のときの周波数特性が求められたりする。

周波数ワーピング

$\tilde z^{-1}=\frac{z^{-1}-\alpha}{1-\alpha z^{-1}}(|a|\lt1)$ によって、またしても変数変換をする。(この式はオールパスフィルタそのもの。) $z^{-1}$ について解けば、 $z^{-1}=\frac{\tilde z^{-1}+\alpha}{1+\alpha \tilde z^{-1}}$ となる。

変数変換なので、対応を確認することにする。

$\frac{-2\alpha + j(1-\alpha^{2})}{1+\alpha^{2}}$ ってなんだよという感じだが、この値は絶対値を取ると $\alpha$ に関わらず 1 になる。実は、 $z$ が単位円上にあれば、 $\tilde z$ も必ず単位円上にある。これは以下のように確認できる。(この式はオールパスフィルタの振幅特性が平坦なことを示すものと同じ。)

式でも確認できるけど、 $\tilde z^{-1}=\frac{z^{-1}+\alpha}{1+\alpha z^{-1}}$ の分母の絶対値と分子の絶対値が等しいことを言えばいいだけので、図を書いたほうがわかりやすい。図は $\alpha \lt 0$ なので注意。

f:id:nagiss:20220201171836p:plain

$z=\cos \omega + j\sin\omega$ のとき、 $\tilde z$ の絶対値は $1$ であることがわかったので、新しい変数 $\tilde\omega$ を導入して $\tilde z = \cos \tilde\omega + j\sin\tilde\omega$ と表すことができる。 $z=\cos \omega + j\sin\omega$ のとき $\tilde\omega$ はどんな値になるかも確認したい。(この計算はオールパスフィルタの位相特性の計算と同じ。) さっきと同じように図で考えると、 $\omega$ と $\tilde\omega$ は下図の角度になる。

f:id:nagiss:20220201175536p:plain

図から、高校数学を駆使して $\omega$ と $\tilde\omega$ の関係を計算する。筋のいい方法がなかなか分からなくて大変だった。ネタバレをすると、図の二等辺三角形の各辺の長さを求めて余弦定理を使うと良い。結果は、
$\displaystyle \cos\tilde\omega = \frac{-2\alpha + (1+\alpha^{2})\cos\omega}{(1+\alpha^{2})-2\alpha\cos\omega}$
となる。 $-\pi\lt\omega\lt\pi$ の範囲では、 $\omega$ と $\tilde\omega$ の符号が等しいことを使えば、 $\omega$ から $\tilde\omega$ が一意に定まる。左辺を $\tan\tilde\omega$ にした場合は、
$\displaystyle \tan\tilde\omega = \frac{(1-\alpha^{2})\sin\omega}{(1+\alpha^{2})\cos\omega-2\alpha}$
となる。

$\omega$ と $\tilde\omega$ の対応を描画してみるとこんな感じ。

f:id:nagiss:20220201221057p:plain

描画コード

import numpy as np
import matplotlib.pyplot as plt

plt.figure(figsize=(4, 4), dpi=200)
for alpha in [0.8, 0.4, 0.0, -0.4, -0.8]:
    omega = np.linspace(0, np.pi, 1000)
    omega_tilde = np.arccos(
        (-2 * alpha + (1 + alpha * alpha) * np.cos(omega)) / ((1 + alpha * alpha) - 2 * alpha * np.cos(omega))
    )
    plt.plot(omega, omega_tilde, label=fr"$\alpha={alpha:.1f}$")
plt.legend()
plt.xlabel(r"$\omega$")
plt.ylabel(r"$\widetilde\omega$")
plt.xlim(0, np.pi)
plt.ylim(0, np.pi)
plt.xticks(np.linspace(0, np.pi, 5), [r"$0$", r"$\pi/4$", r"$\pi/2$", r"$3\pi/4$", r"$\pi$"])
plt.yticks(np.linspace(0, np.pi, 5), [r"$0$", r"$\pi/4$", r"$\pi/2$", r"$3\pi/4$", r"$\pi$"])
plt.grid()
plt.show()

結局、ここまで長い道のりを辿ってやったのは、 $\omega$ で表されたフーリエ変換の式 $F(\omega)$ を、変数変換によって $\tilde\omega$ で表したことである。 $\omega$ と $\tilde\omega$ の関係はさっき計算した通りなので、 $\tilde\omega$ の値から周波数特性を求められるようになった。

オールパスフィルタ

適当なフィルタ $H(z)$ を考える。イメージしやすいようにとりあえず $H(z)=1-0.2z^{-1}$ のローパスフィルタとしておく。このフィルタの周波数特性 $F_{H}(\omega)$ は、
$F_{H}(\omega)=H(\cos\omega+j\sin\omega)=1-\frac{0.2}{\cos\omega+j\sin\omega}$
と求められる。

ここで、唐突にフィルタ内の遅延要素 $z^{-1}$ をオールパスフィルタ $\frac{z^{-1}-\alpha}{1-\alpha z^{-1}}$ で置き換えたフィルタ $G(z)$ を考えると、
$G(z) = 1-0.2\frac{z^{-1}-\alpha}{1-\alpha z^{-1}}$
となる。

この式にさっきの変数変換を適用すれば、
$G(z)=1-0.2\tilde z^{-1} = H(\tilde z)$
と表せる。

$H(\tilde z)$ の周波数特性を考える。置き換える前のフィルタ $H(z)$ の周波数特性は $F_{H}(\omega)=1-\frac{0.2}{\cos\omega+j\sin\omega}$ なので、置き換えたフィルタの周波数特性は同様の計算で $F_{G}(\omega)=1-\frac{0.2}{\cos\tilde\omega+j\sin\tilde\omega}$ となる。

これにより、遅延要素をオールパスフィルタで置き換えたフィルタの特性は、元のフィルタの周波数特性の、周波数の軸を伸縮させたものになっているとわかる。

メル目盛の近似

サンプリング周波数が 16kHz のとき、 $\alpha=0.42$ くらいにするとメル目盛をよく近似する周波数目盛になる……とよく言われるけど 0.42 の根拠をあんまり見たことが無かった気がするので、描画して確認してみる。

f:id:nagiss:20220201223757p:plain

描画コード

plt.figure(figsize=(4, 4), dpi=200)

alpha = 0.42
sample_rate = 16000

omega = np.linspace(0, np.pi, 1000)
omega_tilde = np.arccos(
    (-2 * alpha + (1 + alpha * alpha) * np.cos(omega)) / ((1 + alpha * alpha) - 2 * alpha * np.cos(omega))
)
plt.plot(omega, omega_tilde, label=fr"$\alpha={alpha:.2f}$")

freq = omega * sample_rate / (2.0 * np.pi)
mel = 1000 * np.log2(freq / 1000.0 + 1.0)
max_mel = mel[-1]
target_omega = mel * (np.pi / max_mel)
plt.plot(omega, target_omega, color="black", label="mel")

plt.legend()
plt.xlabel(r"$\omega$")
plt.ylabel(r"$\widetilde\omega$")
plt.xlim(0, np.pi)
plt.ylim(0, np.pi)
plt.xticks(np.linspace(0, np.pi, 5), [r"$0$", r"$\pi/4$", r"$\pi/2$", r"$3\pi/4$", r"$\pi$"])
plt.yticks(np.linspace(0, np.pi, 5), [r"$0$", r"$\pi/4$", r"$\pi/2$", r"$3\pi/4$", r"$\pi$"])
plt.grid()
plt.show()

よく近似できてるっぽい。

おわりに

SPTK の freqt がどういう仕組みになっているのかわかりません　だれかおしえてください