« RIAAフォノイコライザを FIR デジタルフィルタで実装する ② 係数計算編 | トップページ | RIAAフォノイコライザを FIR デジタルフィルタで実装する ④ 実践&聴き比べ編 »

2024年12月16日 (月)

RIAAフォノイコライザを FIR デジタルフィルタで実装する ③ シミュレーション編

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

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

この記事はアナログレコード再生の際に使用するRIAAフォノイコライザをFIRフィルタで実装し、シミュレーションと実際の音出しを行って評価を行うことを目的とするものです。


前回までで、今回の目標である2種類のFIRフィルタの係数が計算できました。
2種類のFIRフィルタとはすなわち、直線位相タイプと非直線位相タイプを指します。
RIAA特性はアナログフィルタの時定数で定義されており、当然ながら位相も含めて録音時の逆特性で再生することが求められるため、RIAAイコライザをデジタルフィルタで再現する場合もアナログRIAAイコライザと等価な特性を目指すべきです。

しかしながらFIRデジタルフィルタでは、アナログでは実現できない直線位相フィルタが構成できるため、FIRフィルタで直線位相のRIAAイコライザを構成したらどのような音になるのか?アナログの特性に近似の、非直線位相のFIRフィルタや、IIRフィルタで構成した場合とで聴感上の違いがあるのか?という検証を行うことが今回の目的です。

それではさっそく始めましょう。

前回、直線位相FIRと非直線位相FIRの2種類のRIAAの係数を計算し、それぞれ8192タップの係数のファイルを得ました。
今回はこれらの係数ファイルを使って、特性のグラフ表示と、実際に音楽信号(wavファイル)を処理した出力(wavファイル)を得ることを目標とします。
今回もOctaveを使って進めます。

まずはFIRフィルタの係数ファイルと入力のwavファイルを与えると、結果をwavファイルで出力するプログラムです。

%%%%%ここから

pkg load signal;

% FIRフィルタの係数を読み込み
coeff = load('coeff_nonlin6db_48000.txt');                                 %  ①使用する係数ファイル

% 音声ファイルの読み込み
inputFile = "alice4816_direct_15s.wav";                                    %  ②処理するwavファイル
[~,name,ext] = fileparts(inputFile);
inputFileName = [name,ext];
[input_signal, fs] = audioread(inputFile);                                  % 入力信号とサンプリング周波数を取得

% フィルタリング
filtered_signal = filter(coeff, 1, input_signal);                           % フィルタ適用

% 出力ファイルを保存
outputFileName = sprintf('FIR_RIAA_%s',inputFileName);
audiowrite(outputFileName, filtered_signal, fs);

% 周波数特性を計算
[h, w] = freqz(coeff, 1, 4096, fs);                                          % 周波数特性を計算

% ゲイン(dB)と位相(degree)を計算
gain = 20 * log10(abs(h));                                                   % ゲイン(dB)
phase = angle(h) * (180 / pi);                                              % 位相(degree)

% プロット
figure;

% ゲイン特性
subplot(2, 1, 1);
semilogx(w,gain);
grid;
title(sprintf('FIR RIAA EQ %dHz Sampling',fs));
xlabel('Frequency (Hz)');
ylabel('Gain (dB)');
xlim([10, 30000]);
ylim([-20, 30]);                  % 必要に応じて調整

% 位相特性
subplot(2, 1, 2);
semilogx(w, phase);
grid;
xlabel('Frequency (Hz)');
ylabel('Phase (Deg)');
xlim([10, 30000]);

%%%%%ここまで

使い方は簡単で、まず上のコードをOctaveのエディタにコピペして、適当な名前でOctaveのワークスペースのディレクトリに保存します。
次に、使用する係数ファイルを①のところに記述します。前回のファイルはcoeff_lin6db_48000.txtとcoeff_nonlin6db_48000.txtの2種類でしたので、都度書き換えて使います。
次に、実際に処理したいwavファイルを②のところに指定します。
これらのファイルはすべてOctaveのワークスペースに配置します。
処理したいwavファイルが他の場所にある場合は、フルパスで指定します。形式は"c:/AAA/bbb/CCC.wav"のような書式で、バックスラッシュではなくスラッシュなので注意が必要です。フルパスを打ち込むのが面倒な場合は、Octaveのコマンドウィンドウを表示して、そこに目的のファイルをドラッグ&ドロップすればフルパスが表示されるので、それをコピペすると便利です。

あとは実行ボタンを押せば、特性が表示され、"FIR_RIAA_(入力ファイル名)"でRIAA処理されたファイルが出力されます。
今回は48kHzサンプリングでレコードからダイレクト録音したソースを使って検証しており、処理時間はだいたいwavファイルの再生時間と同じくらいかかります。つまり5分の音楽ファイルなら、処理に5分ほどかかりました。特性は図1のように表示されます。
実行ボタンを押したらすぐにコマンドウィンドウでエラーが出ていないか確認してください。エラーが出ているといくら待っても何も起こりません(^-^;

Fir_riaa_nonlinear_towav

図1.表示された非直線位相FIRの特性


直線位相タイプの係数も入れてみましょう。
リストの①のところにcoeff_lin6db_48000.txtと入力して実行します。図2に表示された特性を示します。


Fir_riaa_linear_towav
図2.表示された直線位相FIRの特性

ここまでで、処理されたwavファイルも出力されたと思います。出力wavファイルは同名の場合は上書きされてしまいますので、都度どこかに待避させてください。


次に、IIRフィルタで同じ処理を行ってみましょう。プログラムは次のとおりです。

%%%%%ここから

pkg load signal;

% フィルタ係数を定義
%{
%96kHz プリワーピングなし +6dB
#ゲインはb0のみで調整
b0 = 0.26;
b1 = [1 0];
b2 = [1 -0.96777105 0];
a1 = [1 0];
a2 = [1 -1.866859545 0.867284263];
%}

%
%48kHz プリワーピングなし +6dB
#ゲインはb0のみで調整
b0 = 0.5;
b1 = [1 0];
b2 = [1 -0.936564324 0];
a1 = [1 0];
a2 = [1 -1.749567588 0.751160265];
%}

% IIRフィルタ係数の計算
b = b0 * conv(b1, b2);
a = conv(a1, a2);

% 音声ファイルの読み込み
inputFile = "alice4816_direct_15s.wav";                  %③
[~,name,ext] = fileparts(inputFile);
inputFileName = [name,ext];
inputFileName
[input_signal, fs] = audioread(inputFile);

% フィルタリング
filtered_signal = filter(b, a, input_signal);

% 出力ファイルを保存
outputFileName = sprintf('IIR_RIAA_%s',inputFileName);
audiowrite(outputFileName, filtered_signal, fs);

% 周波数特性を計算
[h, w] = freqz(b, a, 4096, 'whole', fs);

% ゲイン(dB)と位相(degree)を計算
gain = 20 * log10(abs(h)); % ゲイン (dB)
phase = angle(h) * (180 / pi); % 位相 (degree)

% プロット
figure;

%ゲイン特性
subplot(2,1,1);
semilogx(w,gain);
grid;
title(sprintf('IIR RIAA EQ %dHz Sampling',fs));
xlabel('Frequency (Hz)');
ylabel('Gain (dB)');
xlim([10 30000]);

%位相特性
subplot(2,1,2);
semilogx(w,phase);
grid;
xlabel('Frequency (Hz)');
ylabel('Phase (Deg)');
xlim([10 30000]);
ylim([-100, 0]); % 必要に応じて調整

%%%%%ここまで

このコードでは、96kHzサンプリングと48kHzサンプリングの2種類のフィルタ係数が書かれていて、96kHzはコメントアウトされ、48kHzサンプリングが使えるようになっています。IIRフィルタでの係数計算については
RIAAフォノイコライザをIIRデジタルフィルタで実装する ②係数計算編
RIAAフォノイコライザをIIRデジタルフィルタで実装する ③シミュレーション編
を参照してください。

処理したいwavファイルを③のところに指定して実行すれば、同じように"IIR_RIAA_(入力ファイル名)"でRIAA処理されたファイルが出力されます。
このIIRフィルタの処理は時間がほとんどかからず、即座に結果が出力されます。

図3に表示された特性を示します。


Iir_riaa_towav_20241214210201


図3.表示されたIIRフィルタによるRIAAの特性


これから何曲か処理して、聴き比べをしてみようと思います。はたして聴感的に差があるのか……!?

| |

« RIAAフォノイコライザを FIR デジタルフィルタで実装する ② 係数計算編 | トップページ | RIAAフォノイコライザを FIR デジタルフィルタで実装する ④ 実践&聴き比べ編 »

コメント

コメントを書く



(ウェブ上には掲載しません)


コメントは記事投稿者が公開するまで表示されません。



« RIAAフォノイコライザを FIR デジタルフィルタで実装する ② 係数計算編 | トップページ | RIAAフォノイコライザを FIR デジタルフィルタで実装する ④ 実践&聴き比べ編 »