アストロノーカで遊んで興味を持ったのが遺伝的アルゴリズム。その考え方を使ってフィルタ設計でもできないかと思ったので、とりあえずやってみました。中に出てくるbiquad_*は、Robert Bristow-Johnsonによるフィルター設計レシピの実装で、中心周波数やゲインを入力するだけでデジタルフィルター係数が出てくるというもの。tfcとtmxはそれぞれターゲットとしている周波数と、そこでの振幅値です。
まずフィルタをランダムに作ります。ローシェルフ×1、ハイシェルフ×1、ピーキングEQ×N個の、合計N+2個のフィルタが1セットです。それをMセット作ります。ランダムに作ったので、最初のうちはターゲットとするフィルタとは似ても似つかないものです。
ターゲットフィルタとMセットのフィルタを比較して、あてはまりの良いもの*1だけを選抜します。それらのフィルタどうしを交配させて次世代のMセットをつくり、ターゲットと比較します。再度あてはまりの良いものだけを残して交配させて・・・ということを50世代ほど繰り返すと、あら不思議、なんとなくターゲットに似た形状のフィルタができあがります。
運が良いと当てはまりの良いフィルタができますが、運が悪いとそこそこのものしかできません。交配の仕方によってもだいぶ結果に違いが出てきます。探せば、同様の方法でフィルタの最適化をしている例はたくさん出てくると思いますし、僕みたいな思いつきでやってみたものよりも良い結果になっているでしょう。アストロノーカ面白かった→遺伝的プログラミングの解説本を読んでみた→プログラム書いてみる、という道筋はあまりにも乱暴でしたねー。失礼しました。
% This program attempts to calculate the best fitting filter from a % combination of a low-shelf filter, a high-shelf filter, and multiple peaking % filters. All basic filters are biquad filters described by Robert % Bristow-Johnson. %[b,a] = biquad_peak(fc, gain, q, fs) %[b,a] = biquad_lowshelf(fc, gain, s, fs) %[b,a] = biquad_highshelf(fc, gain, s, fs) N = 6; % use up to N biquad peaking filters M = 16; % number of filters to compete against eath other fs = 44100; % example target filter tfc = [ 0.0010 0.0012 0.0016 0.0020 0.0025 0.0031 0.0039 0.0050 0.0063 0.0079 0.0099 0.0125 0.0157 0.0198 0.0250 0.0315 0.0397 0.0500 0.0630 0.0794 0.1000 0.1260 0.1587 0.2000 0.2520 0.3175 0.4000 0.5040 0.6350 0.8000 1.0079 1.2699 1.6000 2.0159 ] * 10000; tmx = db2mag([ -26.6999 -21.4950 -17.1127 -14.2987 -11.8276 -8.7456 -6.2109 -4.3881 -4.0955 -3.6969 -4.1158 -4.8928 -5.7010 -6.3959 -7.1127 -7.3416 -7.4848 -7.1744 -5.9094 -3.8359 -2.1256 -0.0946 3.1514 7.2099 9.8199 7.7844 4.1823 4.9488 2.0716 2.1853 3.5547 -6.4740 -15.6247 -18.1134 ]); %% initial randomization of filters filt(M) = struct('fc', [], 'gain', [], 'q', []); fcmax = 16000; fcmin = 1/fcmax*1000; gainmax = +20; gainmin = -gainmax; qmax = 100; qmin = 1/qmax; for n=1:length(filt) % create N peaks, 1 low-shelf, and 1 high-shelf biquad filter parameters filt(n).fc = sort(rand(N+2,1) * (fcmax-fcmin) + fcmin); filt(n).gain = rand(N+2,1) * (gainmax-gainmin) + gainmin; filt(n).q = rand(N+2,1) * (qmax-qmin) + qmin; filt(n).q(1) = 1; filt(n).q(N+2) = 1; end %% test filters for kkk=1:50 score = zeros(length(filt),1); for n=1:length(filt) f = filt(n); A = zeros(N+2, 3); B = zeros(N+2, 3); [B(1,:),A(1,:)] = biquad_lowshelf(f.fc(1), f.gain(1), f.q(1), fs); [B(N+2,:),A(N+2,:)] = biquad_highshelf(f.fc(N+2), f.gain(N+2), f.q(N+2), fs); for k=2:N+1 [B(k,:),A(k,:)] = biquad_peak(f.fc(k), f.gain(k), f.q(k), fs); end x = [1 ; zeros(fs-1,1)]; for k=1:N+2 x = filter(B(k,:), A(k,:), x); end [f3,mx3] = magnitude(x,fs,3); score(n) = sqrt(sum((mag2db(tmx)-mx3).^2)); end [tmp, ind] = sort(score); fprintf(1, '%3d %f\n', kkk, min(score)); %% improve filters ff = filt(ind(1)); for n=2:length(filt) filt(n).fc = (rand(N+2,1)*0.6+0.7) .* ff.fc; filt(n).gain = (rand(N+2,1)*0.6+0.7) .* ff.gain; filt(n).q = (rand(N+2,1)*0.6+0.7) .* ff.q; end filt(ind(1)) = ff; filt(ind(M)).fc = rand(N+2,1) * (fcmax-fcmin) + fcmin; filt(ind(M)).gain = rand(N+2,1) * (gainmax-gainmin) + gainmin; filt(ind(M)).q = rand(N+2,1) * (qmax-qmin) + qmin; end %% display results for n=1:length(filt) f = filt(n); A = zeros(N+2, 3); B = zeros(N+2, 3); [B(1,:),A(1,:)] = biquad_lowshelf(f.fc(1), f.gain(1), f.q(1), fs); [B(N+2,:),A(N+2,:)] = biquad_highshelf(f.fc(N+2), f.gain(N+2), f.q(N+2), fs); for k=2:N+1 [B(k,:),A(k,:)] = biquad_peak(f.fc(k), f.gain(k), f.q(k), fs); end x = [1 ; zeros(fs-1,1)]; for k=1:N+2 x = filter(B(k,:), A(k,:), x); end subplot(4,4,n); magnitude(x,fs,3); hold on; semilogx(tfc, mag2db(tmx), 'k--'); hold off; [f3,mx3] = magnitude(x,fs,3); score(n) = sqrt(sum((mag2db(tmx)-mx3).^2)); end [tmp,ind] = sort(score); disp([kkk min(score)]); for n=1:length(filt) subplot(4,4,ind(n)); title(num2str(n)); set(gca, 'ylim', [-12 12]); end
