読者です 読者をやめる 読者になる 読者になる

遺伝的アルゴリズム的な何かでフィルターを設計してみる(でもあまりうまくいかない)

アストロノーカで遊んで興味を持ったのが遺伝的アルゴリズム。その考え方を使ってフィルタ設計でもできないかと思ったので、とりあえずやってみました。中に出てくる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

*1:あてはまりの良さとしては、とりあえず周波数帯域ごとのデシベル値の差分を取って二乗和を計算しています