第4回 Tuning #4

( 変更履歴:
v0.3 ['16/10/11] セント値のグラフからインハーモニシティ値を加え無くしました。
v0.2 ['11/10/30] c2fの計算式を変更しました。)

Javaの run部分 tuningの計算を取り出してみました。
mwdファイルとうなり数を

>> tuning('design.mwd', 1.0)

とする事で tuning後のセント値のグラフが表示されます。

function毎に別ファイルにする事で計算速度は早くなります。

827.98 -> 39.638[sec]

function tuning(fname, beat)

global SM KM NG WM AA DC W head is_fortepiano;
SM = 7.85;          # 芯線密度 g/cm^3
KM = 6.9906;        # 銅と空気の混合密度 g/cm^3
NG = 9.80665204821; # 重力加速度 kgw
WM = 9.0;           # 巻線密度 g/cm^3
DC = 0.1;           # Δcent値

if (nargin < 1 || nargin > 2)
  usage("tuning('mwd_file_name' [, beats (1.0)])");
elseif (nargin < 2)
  beat = 1.0;
end

[head, wire] = ReadWire(fname);
W = Wire(wire);

if (length(head) > 3)    # フォルテピアノ
  is_fortepiano = true;
  start_key = head(4)+1;
  AA = 415;
else                     # ピアノ
  is_fortepiano = false;
  start_key = 1;
  AA = 440;
end

lenw = length(wire);
bee = zeros(1, lenw);

SetWire(start_key, AA);

k2k = 24;             # 2octave tune
inte = Interval(k2k);

wi36 = k2w(36);
wi37 = k2w(37);
wi61 = k2w(61);
wi62 = k2w(62);
up = sa = sb = ds = 0;

do
  ds = up/(wi61-wi37);
  set2octave(up, wi37, wi61);

  tune(true, wi62, inte, beat);
  sa = W(wi62).cnt-W(wi61).cnt;

  tune(false, wi36, inte, beat);
  sb = W(wi37).cnt-W(wi36).cnt;
  if (sa-ds > 0.3)
    d = DC*10;
  else
    d = DC;
  end
  up = up+d;

until (sa <= ds && sb >= ds)
up = up-d;
set2octave(up, wi37, wi61);

for i=wi36:-1:start_key
  tune(false, i, inte, beat);
end
for i=wi62:lenw
  tune(true, i, inte, beat);
end

cent = zeros(1, lenw);
for i=start_key:lenw
  cent(i) = W(i).cnt;
end

plot(cent, '@;Cent;');
end
#-----

function set2octave(up, i_start, i_end)
  global W;

  ds = up/(i_end-i_start);
  dd = up/2;
  n = 0;
  for i=i_start:i_end
    W(i).cnt = ds*(n++)-dd;
    ResetWire(i);
  end
end  # 2octave分のセント値を加えます
#-----

function beat = tune(sw, wi, inte, max)
  global W DC;

  if (sw == true) # 上行
    dc = DC;
    wh = wi;
    wl = itvl2wire(wi, -inte(1));
  else            # 下行
    dc = -DC;
    wl = wi;
    wh = itvl2wire(wi, inte(1));
  end

  sa = 0.2; % 差がsa以内でdcを1/10に変更
  beat = GetWbeats(wh, wl, inte);
  while (beat < max)
    if (max-beat > sa)
      d = dc*10;
    else
      d = dc;
    end;

    W(wi).cnt = W(wi).cnt+d;
    ResetWire(wi);
    beat = GetWbeats(wh, wl, inte);
  end

  while (beat > max)
    if (max-beat > sa)
      d = dc*10;
    else
      d = dc;
    end
    W(wi).cnt = W(wi).cnt-d;
    ResetWire(wi);
    beat = GetWbeats(wh, wl, inte);
  end
end  # 一定のうなりになるまでセント値を増減します
#-----

function beats = GetWbeats(wh, wl, inte)
  global W;

  beats = c2f(W(wh).frq*inte(3), W(wh).inh*inte(3)*inte(3)+W(wh).cnt)-\
      c2f(W(wl).frq*inte(2), W(wl).inh*inte(2)*inte(2)+W(wl).cnt);
end  # 音程間のうなりを計算します
#-----

function [header, wire] = ReadWire(fname)
  wire = dlmread(fname);

  len = length(wire(1, :));
  if (len <= 3)
    header = [wire(1,1) wire(1,2) wire(1,3)];
  elseif (len > 3)
    header = [wire(1,1) wire(1,2) wire(1,3) wire(1,4) wire(1,5)];
  else
    error('Not Wire data.'); 
  end
  wire(1,:) = [];
end  # mwdファイルを読み込みます
#-----

function w = Wire(x)
  [row, col] = size(x);

  for i=1:row
    if (col > 3)
      w(i) = struct('length', x(i,1), 'bante', x(i,2), 'copper', x(i,3),\
                    'density', x(i,4), 'pitch', x(i,5),\
                    'key', 0, 'dia', 0, 'ald', 0, 'den', 0,\
                    'ten', 0,'inh', 0, 'cnt', 0, 'frq', 0);
    else
      w(i) = struct('length', x(i,1), 'bante', x(i,2), 'copper', x(i,3),\
                    'density', 0, 'pitch', 0,\
                    'key', 0, 'dia', 0, 'ald', 0, 'den', 0,\
                    'ten', 0, 'inh', 0, 'cnt', 0, 'frq', 0);
    end
  end
end  # 読み込んだmwdからwireデータを作成します
#-----

function SetWire(start_key, aa)
  global SM W is_fortepiano;

  for i=start_key:length(W)
    W(i).key = w2k(i);               # key.no
    W(i).dia = getDia(W(i).bante);   # dia
    W(i).ald = GetAlld(W(i).bante, W(i).copper); # all-d
    W(i).den = GetDensity(W(i).dia, W(i).copper, W(i).ald,\
                          W(i).density, W(i).pitch); # den
    W(i).frq = c2f(GetFreq(W(i).key, aa), W(i).cnt); # freq.
    [gten, W(i).ten] = GetTension(W(i).frq, W(i).length, W(i).bante,\
                                  W(i).copper, W(i).ald, W(i).den); # ten
    if (is_fortepiano)
      sm = W(i).density;
    else
      sm = SM;
    end
    W(i).inh = GetInha(W(i).frq, W(i).length, W(i).bante, W(i).copper,\
                       sm, W(i).ten, W(i).ald); # inha
  end
end  # wireにtuningに必要な値いを設定します
#-----

function ResetWire(wi)
  global SM W is_fortepiano;

  freq = c2f(W(wi).frq, W(wi).inh+W(wi).cnt); #
  [gten, W(wi).ten] = GetTension(freq, W(wi).length, W(wi).bante, W(wi).copper,\
                                 W(wi).ald, W(wi).den); # ten
  if (is_fortepiano)
    sm = W(wi).density;
  else
    sm = SM;
  end
  W(wi).inh = GetInha(freq, W(wi).length, W(wi).bante, W(wi).copper,\
                      sm, W(wi).ten, W(wi).ald); # inha
end  # セント値の変更後に wireの値を再計算します
#-----

function ad = GetAlld(bante, copper)
  global is_fortepiano;

  if (is_fortepiano)
    ad = bante;
  else
    ad = copper*0.94*2+GetDia(bante);
  end
end  # 巻線を加えた時の弦径を計算します
#-----

function den = GetDensity(dia, copper, alld, density, wpitch)
  global WM KM SM is_fortepiano;

  mm = 20;

  if (copper > 0)
    if (is_fortepiano)
      den = (density*(dia/mm)^2+WM*(copper/mm)^2*\
            ((dia+copper)/mm)*pi*10/wpitch)/(dia/mm)^2;
    else
      j = (alld/2)^2;
      g = (dia/2)^2;
      den = (KM*(j-g)+SM*g)/j;
    end
  else
    if (is_fortepiano)
      den = density;
    else
      den = SM;
    end
  end
end  # 弦密度を計算します 
#-----

function [dia, aa] = GetDia(bante)
  ba = 12:0.5:25;
  dm = 1.25:0.05:1.500;
  da = 0.725:0.025:1.225;
  di = [da dm];

  dia = -1;
  aa = 440;

  if (bante < 12 && bante > 0.1)
    dia = bante;
    aa = 415;
  else
    x = find(ba == bante);
    if (x > 0)
      dia = di(x);
    end
  end
end  # 番手から弦径を得ます 
#-----

function inh = GetInha(freq, leng, ban, cop, sm, ten, alld)
  me = 1000;
  cm = 10;
  [di, aa] = GetDia(ban);

  if (cop > 0)
  inh = 3.3*cm^13*(di/cm)^2/(leng/cm)^4/\
        ((ten/pi/(sm*me))^(1/2)/(leng/me)/(di/me))^2*\
        (1+alld/di/cm);
  else
    inh = 3.3*cm^13*(di/cm)^2/(leng/cm)^4/freq^2;
  end
end  # インハーモニシティ値を計算します
#-----

function [gten, ten] = GetTension(freq, leng, ban, cop, alld, den)
  global SM NG is_fortepiano;

  mm = 1000;
  [dia, pitch] = GetDia(ban);

  if (is_fortepiano)
    di = dia;
  else
    if (cop > 0)
      di = alld;
    else
      di = dia;
    end
  end

  if (cop > 0)
    de = den;
  else
    if (is_fortepiano)
      de = den;
    else
      de = SM;
    end
  end

  ten = pi*de*mm*freq^2*(leng/mm)^2*(di/mm)^2;
  gten = ten/NG;
end  # 弦張力を計算します
#-----

function freq = GetFreq(key, pitch)
  if (nargin < 1)
    usage("GetFreq([key | key_list], [pitch (440)])");
  elseif (nargin < 2)
    pitch = 440;
  end

  freq = pitch.*2.^((key-49)./12);
end  # キーから周波数を計算します
#-----

function f = c2f(freq, cent)
  f = freq.*2.^(cent/1200);
end  # セント値から周波数を計算します
#-----

function wi = k2w(k)
  global head;
  ob = head(1);
  db = head(2);
  sb = head(3);

  if (k < ob)
    wi = k-1;
  elseif (k >= sb)
    wi = k-sb+db;
  else
    wi = (k-ob)*2+ob-1;
  end
  wi = wi+1;
end  # キーからwire番号を得ます
#-----

function key = w2k(w)
  global head;
  ob = head(1);
  db = head(2);
  sb = head(3);

  if (w < db)
    if (w >= ob)
      key = round(ob+1+(w-ob)/2);
    else
      key = w+1;
    end
  else
    key = w-db+sb;
  end

  key = key-1;
end  # wire番号からキーを得ます
#-----

function x = itvl2wire(wi, inte)
  nk = w2k(wi);
  x = wi+inte;
  xk = w2k(x);

  while(nk+inte > xk)
    xk = w2k(++x);
  end
  while(nk+inte < xk)
    xk = w2k(--x);
  end
end  # 音程からwire番号を得ます
#-----

function [inter, intst] = Interval(inte)
  intn = [1 16 15; 2 9 8; 3 6 5; 4 5 4; 5 4 3; 6 7 5;\
          7 3 2; 8 8 5; 9 5 3; 10 7 4; 11 15 8; 12 2 1;\
          16 5 2; 19 3 1; 24 4 1; 28 5 1; 31 6 1; 36 8 1];

  ratio = {'mi.2nd' 'Mj.2nd' 'mi.3rd' 'Mj.3rd' 'P.4th' 'Dim.5th'\
           'P.5th' 'mi.6th' 'Mj.6th' 'mi.7th' 'Mj.7th' 'Octave'\
           '10th' '12th' '2octave' 'Oct.10th' 'Oct.12th' '3octave'};

  inter = -1;
  intst = '';

  x = find(intn(:, 1) == inte);
  if (x > 0)
    inter = intn(x, 1:3);
    intst = ratio{x};
  end
end  # 音程からキー間隔｜倍音数を得ます