ウラギンシジミの日向ぼっこ

きょうは2024年の大晦日です。
しかもきょうは今月2度目の新月で、同月2度目の満月を指す「ブルームーン」ほど知られていませんが、「ブラックムーン」と呼ばれることもあるそうです。次に大晦日がブラックムーンとなるのは19年後の2043年。遠い未来です。

きょうは日中13度ほどまで気温が上がってポカポカ陽気でした。昼過ぎに高円寺天祖神社に年越しの大祓の人形を預けに行って、帰り道で蝶が飛んでいたので見ていたら、陽の当たる路面に止まって日光浴を始めました。近づいてよく見るとウラギンシジミのメスでした(^-^)

ウラギンシジミは、シジミとはいうものの他のシジミチョウとは異なる特徴を持った蝶で、シジミチョウ科ではなくウラギンシジミ科に分類するべきだという議論が、ぼくが知る限り50年くらい前から続いています。ウラギンシジミは、子供の頃、秋が深まって、なったまま熟れた柿の汁を吸いに集ってくるところを捕まえた記憶があるので、晩秋に発生する蝶だと思い込んでいましたが、実際は5月から10月くらいに発生する蝶で、越冬した成虫が春に見られることもあってほぼ一年中生息している蝶のようです。

2024年もよい年でした(^-^)
ことしはSSDAC基板を応用してレコード再生用のRIAAイコライザを実装し、とくに12月になってからRIAAイコライザをIIRとFIRで実装して聴き比べができたのが印象的でした。アナログレコードはまた人気が出てきているようで新譜も出ているようですね。また世界的に見て日本は音楽売上げに占めるCDの割合がとても高く、これは音質の点から喜ばしいことだと思います。

ことしお会いできたみなさまへ。
ことしも大変お世話になりました。また来年もよろしくお願いします。

よいお年をお迎えください。

写真1．きょう2024年大晦日に目撃したウラギンシジミ（メス）



写真2．2019年2月7日に目撃した越冬中のウラギンシジミ

RIAAフォノイコライザを FIR デジタルフィルタで実装する ④ 実践＆聴き比べ編

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

RIAAフォノイコライザを FIR デジタルフィルタで実装する ① 予告編
RIAAフォノイコライザを FIR デジタルフィルタで実装する ② 係数計算編
RIAAフォノイコライザを FIR デジタルフィルタで実装する ③ シミュレーション編

前回までで、実際にダイレクトリッピングしたwavファイルを入力して、FIRおよびIIRフィルタでRIAAイコライジング処理をする具体的な方法を説明しました。

今回は私がもっている何枚かのアナログレコードからパソコンでダイレクトリッピングしたサンプルソースを紹介し、実際に聞き比べた結果をお知らせしたいと思います。

まずはダイレクトリッピングする際に使用した機材です。図1がそのブロック図、図2は使用したフラットアンプの回路です。


図1．ダイレクトリッピングブロック図



図2．フラットアンプ回路（電源は±12V）


実際のリッピング作業では、SoundForgePro11を使ってファイルフォーマットは96kHz24bitのwavで行い、今回の検証用に、同じSound ForgePro11を使って48kHz16bitにダウンサンプリングしました。
ダウンサンプリングした理由は、今回3種類のフィルタタイプによる比較を行うにあたって、どちらかというとロースペックにした方が差が出やすいのではないかと思ったからです。

今回試聴した曲を何曲か紹介します。以下はすべてダイレクトリッピングしてダウンサンプリングした48kHz16bitフォーマットですので、実際にOctaveでRIAAイコライザを通して聴くことができます。

かっこ内はアーティスト名と収録DISKです。

Minor Swing（Stephane Grappelli/Live in Sanfransisco）

Whisper Not（阿川泰子/Journey)

Autumn Leaves（Bill Evans/PORTRAIT IN JAZZ）

夏服のイブ（松田聖子/Seiko Town）

ロックンルージュ（松田聖子/Seiko Town）

Take Five（Dave Brubeck/TAKE FIVE(Single))

Mambo Mongo（Belmonte & His Afro Latin7/OLE)

NIGHT IN TUNISIA（HOT SALSA/Hot Salsa meets Swedish Jazz）


さて、実際に試聴してみた印象ですが、思ったほど違いはありませんでした。というかほぼ同じです。
ただ、ほんのわずかにですが、FIRの直線位相タイプは女性ボーカルのサ行やパーカッションのハイハット、トランペットの音、バイオリンの弦がこすれる音や、音が細く掠れる場面などが少し物足りないというか、リアリティに欠けるように聞こえる印象がしました。
FIRの非直線位相とIIRはほぼ違いがなく、FIRの直線位相に比べるとリアルな音、という印象です。
ピアノの音や、大編成の演奏の音はほとんど差を感じませんでした。

もっともこれはブラインドテストではなく、自分でどのEQでかけているかわかっていますから、先入観でそう感じた可能性は多いにあります。
機会があれば、少し広めの会場でたくさんの人に聴いてもらって比較をしてみたいと思います。

みなさんもぜひ試してみてくださいね(^-^)

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のように表示されます。
実行ボタンを押したらすぐにコマンドウィンドウでエラーが出ていないか確認してください。エラーが出ているといくら待っても何も起こりません(^-^;



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


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



図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に表示された特性を示します。





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


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

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特性を確認し、実際に音楽を再生します。

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＝１１．７２（Ｈｚ）

となります。

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