RIAAフォノイコライザを FIR デジタルフィルタで実装する ② 係数計算編

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

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


今回はフリーで使える数値解析用プログラミング言語Octaveを使って、アナログレコード再生用のRIAAイコライザのシミュレーションを行うために、まずはRIAAイコライザ回路をFIRフィルタで構成するためのフィルタ係数を求めます。

最終的な目標は、
・IIRフィルタでRIAAイコライザ回路（アナログ等価=非直線位相）
・FIRフィルタでRIAAイコライザ回路（アナログ等価=非直線位相）
・FIRフィルタでRIAAイコライザ回路（直線位相）
という3種類のシミュレーションと試聴を行い、それぞれの特性と実際の音の違いを検証することです。

ではさっそく、FIRフィルタでRIAAイコライザ回路を構成する作業を説明します。
デジタルフィルタとはなにか？などについては上に書いた過去記事やWikiなどを参照してください。

FIRフィルタは、フィルタ係数を対称形にすることで直線位相フィルタを構成できることが大きな特徴で、今回は直線位相と非直線位相（アナログに近似）の2種類のFIRフィルタ設計でRIAAイコライザ回路を構成します。

任意の周波数特性に対し、逆フーリエ変換を行うと、その周波数特性のインパルスを求めることができ、そのインパルスを係数とするFIRフィルタを構成すれば、もとの周波数特性をもつフィルタができます。

この作業を行う際に、元の周波数特性に複素共役対称性を持たせる（=負の周波数領域を追加して、対称性を持たせる）ことで、逆フーリエ変換の結果が対称性を持ち、直線位相のFIRフィルタを構成することができます。

逆に、元の周波数特性に複素共役対称性を適用せず、正の周波数領域のみで逆フーリエ変換を行った場合は、非直線位相（=アナログで構成するフィルタに近似）となります。

それでは実際の作業に入りましょう。作業環境はWindows11にGNU Octave 8.3.0です。

まずはFIRフィルタの特徴である、直線位相特性をもつRIAAイコライザのフィルタ係数生成です。つぎのコードをOctaveのエディタにコピペして実行すれば、RIAAイコライザを構成するFIRフィルタの係数ファイルが生成されます。

%%%%%ここから

pkg load signal;

% サンプリング周波数とフィルタ設定
fs = 48000; 　　　　　　　　　　　　　　% サンプリング周波数 (Hz)
N = 8192; 　　　　　　　　　　　　　　　% FIR フィルタのタップ数
f = linspace(0, fs/2, N/2+1); 　　　　　　% 周波数範囲（片側スペクトル）
current_dir=pwd;

% RIAA伝達関数の定義
T1 = 3.18e-3; % 時定数1
T2 = 318e-6; % 時定数2
T3 = 75e-6; % 時定数3
s = 2 * pi * 1i * f; % s = jω
Hs = (1 + s*T2) ./ ((1 + s*T1) .* (1 + s*T3));

% 時間領域のインパルス応答を逆フーリエ変換で取得
Hf_full = [Hs, conj(flip(Hs(2:end-1)))]; 　　　　　　% 直線位相特性を得るため複素共役対称性を持たせる
h = real(ifft(Hf_full,N)); 　　　　　　　　　　　　　　% インパルス応答を取得
h = circshift(h, -N/2); 　　　　　　　　　　　　　　　% ゼロ位相シフト
h_6db = 20*h; 　　　　　　　　　　　　　　　　　　　%+6dB@1kHzに調整

% FIRフィルタの周波数特性
[Hfir, w] = freqz(h, 1, N, fs);

% gainとphase
gain = 20*log10(20*abs(Hfir)); 　　　　% ゲイン特性 (dB)
phase = angle(Hfir) * (180/pi); 　　　　% 位相特性 (degree)

% FIRフィルタの係数を保存
save_path=[pwd,sprintf("/coeff_lin6db_%d.txt",fs)]
save("-ascii", save_path, "h_6db");
#save(coeff_file, 'h_fir', '-ascii');

% ゲイン特性と位相特性のプロット
figure;
subplot(2,1,1);
semilogx(w, gain);
grid on;
xlim([10 24000]);
ylim([-30 30]);
xlabel('Frequency (Hz)');
ylabel('Gain (dB)');
title(sprintf('Linear FIR RIAA Response %dHzSampling',fs));

subplot(2,1,2);
semilogx(w, phase);
grid on;
xlim([10 24000]);
xlabel('Frequency (Hz)');
ylabel('Phase (Deg)');

% FIRフィルタ係数をファイルに保存
coeff_file = 'current_dir/coeff.txt';
disp(['FIR filter coefficients saved to ', coeff_file]);

%%%%%ここまで

サンプリング周波数ｆｓとタップ数Ｎは任意の数値を入れられますが、今回はｆｓ=48000，N=8192としています。
Nは2^nを選びます。
FIRフィルタの低域特性は分解能ｆｓ／Ｎで決まり、この場合は48000/8192=5.86 (Hz)
となるので、周波数帯域は5.86Hz～24kHzとなります。
ゲインは1kHzで+6dBになるように調整しています。

このコードの中で注目すべき箇所は、
Hf_full = [Hs, conj(flip(Hs(2:end-1)))]; 　　　　　　% 直線位相特性を得るため複素共役対称性を持たせる
h = real(ifft(Hf_full,N)); 　　　　　　　　　　　　　　% インパルス応答を取得
の部分です。直線位相特性を得るために、求める周波数特性に複素共役対称性をもたせて、周波数特性を正負対称にしています。
そして、ifft関数を使って、この周波数特性のインパルスを得ています。

生成されたFIRフィルタの係数はcoeff_lin6db_48000.txtというファイルとして保存されます。これには8192個の係数が書かれています。
係数ファイルが生成されると同時に、この係数を使った場合のFIRフィルタの特性のグラフが表示されます（図1）。



図1．FIRフィルタによるRIAAイコライザ（直線位相）のGain-Phase特性


次に、同じFIRですがこんどはアナログ回路に近い非直線位相のFIRフィルタ係数を生成します。

%%%%%ここから

pkg load signal;

% サンプリング周波数とフィルタ設定
fs = 48000; % サンプリング周波数 (Hz)
N = 8192; % FIR フィルタのタップ数
f = linspace(0, fs/2, N/2+1); % 周波数範囲（片側スペクトル）
current_dir=pwd;

% RIAA伝達関数の定義
T1 = 3.18e-3; % 時定数1
T2 = 318e-6; % 時定数2
T3 = 75e-6; % 時定数3
s = 2 * pi * 1i * f; % s = jω
Hs = (1 + s*T2) ./ ((1 + s*T1) .* (1 + s*T3)); %周波数特性

% 非直線位相を持つフィルタ設計のため、複素共役対称性を適用しない
%Hf_full = Hs; % 複素共役対称性をもたせず、そのまま使用

% 時間領域のインパルス応答を逆フーリエ変換で取得
h = real(ifft(Hs, N)); % 時間領域の係数
h_6db = 36*h; %+6dB@1kHzに調整

% FIRフィルタの周波数特性
[H, f_resp] = freqz(h, 1, N, fs); % 周波数特性の計算

% gainとphase
gain = 20*log10(36*abs(H)); % ゲイン特性 (dB)
phase = angle(H) * (180/pi); % 位相特性 (degree)

% FIRフィルタの係数を保存
save_path=[pwd,sprintf("/coeff_nonlin6db_%d.txt",fs)]
save("-ascii", save_path, "h_6db");

% ゲイン特性と位相特性のプロット
figure;
subplot(2, 1, 1);
semilogx(f_resp, gain, 'b');
xlabel('Frequency (Hz)');
ylabel('Gain (dB)');
title(sprintf('non-Linear FIR RIAA Response %dHzSampling',fs));
xlim([10 24000]);
ylim([-30 30]);
grid on;

subplot(2, 1, 2);
semilogx(f_resp, phase, 'r');
xlabel('Frequency (Hz)');
ylabel('Phase (Deg)');
xlim([10 24000]);
ylim([-100 0]);
grid on;

%%%%%ここまで

使い方は上と同じですが、周波数特性に複素共役対称性を適用せず、正の周波数特性のみを使用し、非直線位相特性になっています。
得られる係数ファイルはcoeff_nonlin6db_48000.txt、特性図は図2のとおりです。



図2．FIRフィルタによるRIAAイコライザ（非直線位相）のGain-Phase特性


これでFIRの係数が計算できたので、次回はその係数をFIRフィルタに実装して、再度Gain-Phase特性を確認し、実際に音楽を再生します。