RIAAフォノイコライザを FIR デジタルフィルタで実装する ① 予告編

この記事は次の参考記事の続編です。
今回はRIAAフォノイコライザをFIRフィルタで実装するシミュレーションを行いました。

【参考記事】
RIAAフォノイコライザをIIRデジタルフィルタで実装する　①準備編
RIAAフォノイコライザをIIRデジタルフィルタで実装する　②係数計算編
RIAAフォノイコライザをIIRデジタルフィルタで実装する　③シミュレーション編
RIAAフォノイコライザをIIRデジタルフィルタで実装する　④実装・評価編
RIAAフォノイコライザをIIRデジタルフィルタで実装する　⑤Octaveによる位相検証

デジタルRIAAフォノイコライザ・サブ基板（基板頒布あり）
SSDAC128 and デジタルRIAAイコライザ デュアル基板（基板頒布あり）

上の参考記事で書きましたが、アナログレコードはRIAAカーブという周波数特性で録音されているので、正しい音で再生するにはRIAAイコライザを通す必要があります。従来はこれをアナログ回路で行っていましたが、IIRデジタルフィルタで構成することも可能で、IIRフィルタではゲイン特性、位相特性ともにアナログフィルタと等価な特性を持たせることができます。

ところで、デジタルフィルタにはIIRの他にFIRというタイプのものがあって、特徴は直線位相特性をもつことです。つまり全帯域において一律にサンプリング時間の遅れのみ発生し、アナログフィルタのように周波数による位相回りが起こらないという特徴があります。

※すみません、最初に記事を書いた時点で間違った思い込みがありました。
正しくは、「係数が対称なFIRフィルタは直線位相になる」であって、係数が対称でなければ直線位相にはなりません。
また直線位相特性が実現できることはFIRフィルタの重要な性質で、IIRでは実現できません。

過去記事ですでに述べたとおりRIAA規格によるレコード録音再生では、ゲイン特性に伴い位相特性も再現する必要があるため、当然ながら位相特性も含めてアナログ回路によるRIAAカーブと等価な特性を再現できるIIRフィルタで構成するのが妥当と考えられます。
ただ、正しい再生ができるかどうかとは別に、FIRフィルタでRIAAカーブを再現するとどのようになるのか、というのは興味があります。

FIRフィルタの係数はどのように求められるかというと、周波数特性を逆フーリエ変換したものがそのインパルス関数となるため、それをFIRフィルタの係数に使うことで実装できます。
ところが、検証に使っていたプログラミング言語”Octave”の逆フーリエ変換"ifft"の使い方がわからず、保留にしていました。

今回別の課題があって、Linux上でsoxやらbrutefirやらで遊んでいたときに、ふと、
「ひょっとしてChatGPT先生に訊いたら、FIRでRIAAイコライザを設計してくれるかも……！？」
と思い、実際にChatGPT先生に訊いてみたら、教えてもらえました＼(^o^)／
なんと、Pythonのnumpyに逆フーリエ変換があって、それが使えるのです！

そういうわけで、ChatGPTに教えてもらったPythonのコードを手直ししたものを以下に示します。


############　以下　コード

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import freqz

# サンプリング周波数
fs = 44100
nyquist = fs / 2

# RIAAカーブの時間定数
t1 = 3180e-6 # 3180 µs
t2 = 318e-6 # 318 µs
t3 = 75e-6 # 75 µs

# 周波数特性の計算
frequencies = np.linspace(10, nyquist, 4096) # 10Hzから始める
s = 2j * np.pi * frequencies # sはフーリエ演算子、jは虚数
riaa = (1 + t2 * s) / ((1 + t1 * s) * (1 + t3 * s)) #RIAAの伝達特性

# 振幅特性（対数スケールでのプロット用）
gain = 20 * np.log10(np.abs(riaa))
gain -= gain[np.argmin(np.abs(frequencies - 1000))] # 1kHzで0dBに正規化

# 位相特性（ラジアンを度に変換）
phase = np.angle(riaa, deg=True)

# FIRフィルタ設計（逆FFT）
N = 8192 # フィルタのタップ数
impulse_response = np.fft.ifft(riaa, n=N) #逆FFTを求める
fir_coeffs = np.real(impulse_response[:N])

# FIR係数を保存
np.savetxt("riaa_coefficients.dat", fir_coeffs)

# FIRフィルタの周波数応答を確認
w, h = freqz(fir_coeffs, worN=4096, fs=fs)
fir_gain = 20 * np.log10(np.abs(h))
fir_gain -= fir_gain[np.argmin(np.abs(w - 1000))] # 1kHzで0dBに正規化
fir_phase = np.angle(h, deg=True)

# グラフ描画
fig, ax1 = plt.subplots(figsize=(12, 8))

# 振幅特性（ゲイン）のプロット
ax1.set_title("RIAA EQ Frequency and Phase Response")
ax1.semilogx(frequencies, gain, label="Ideal RIAA Gain", linewidth=1.5)
ax1.semilogx(w, fir_gain, label="FIR Filter Gain", linewidth=1.5)
ax1.set_xlabel("Frequency (Hz) [Log Scale]")
ax1.set_ylabel("Gain (dB)")
ax1.set_xlim(10, 30000) # 横軸を10Hzから30kHzに設定
ax1.set_ylim(-30, 30) # 縦軸を30dBから-30dBに設定
ax1.grid(which="both", linestyle="--", linewidth=0.5)
ax1.axvline(1000, color="gray", linestyle="--", linewidth=0.8, label="1 kHz Reference")
ax1.legend(loc="upper left")

# 位相特性のプロット
ax2 = ax1.twinx() # 第二軸（位相特性用）
ax2.semilogx(frequencies, phase, label="Ideal RIAA Phase", color="tab:orange", linestyle="--", linewidth=1.5)
ax2.semilogx(w, fir_phase, label="FIR Filter Phase", color="tab:red", linestyle="--", linewidth=1.5)
ax2.set_ylabel("Phase (degrees)")
ax2.set_ylim(-180, 180) # 位相特性の範囲を-180度から180度に設定
ax2.legend(loc="lower right")

plt.show()

#############ここまで


この中でキモとなるのは、

s = 2j * np.pi * frequencies#フーリエ演算子ｓ
riaa = (1 + t2 * s) / ((1 + t1 * s) * (1 + t3 * s))#RIAAの伝達関数

の、RIAA特性の伝達関数の部分と、

# FIRフィルタ設計（逆FFT）
N = 16384 # フィルタのタップ数
impulse_response = np.fft.ifft(riaa, n=N) #逆FFTを求める
fir_coeffs = np.real(impulse_response[:N])

# FIR係数を保存
np.savetxt("riaa_coefficients.dat", fir_coeffs)

ここの、逆FFTを行ってインパルス応答を求めて、FIR係数として保存している部分です。
このFIR係数を使えば実際のFIRフィルタに実装できます。

このプログラムを実行した結果は次のとおり。

図1．FIRで実装したRIAAイコライザの特性

図1は、青いプロットがRIAA特性の理論値、橙のプロットがFIRフィルタのシミュレーション結果です。
位相特性も含めてなかなかきれいに出ていますね。
低域で誤差が大きく出ているのは、サンプリング周波数（Ｆｓ）とタップ数（Ｎ）によって決まる分解能によるものです。この例では分解能Δｆは

Δｆ＝Ｆｓ／Ｎ＝96000／8192＝１１．７２（Ｈｚ）

となります。

そういうわけで、音が出せる状態になったらまた報告します(^-^)