月別アーカイブ: 2020年12月

マイクロプロセッサで立ち上がり・減衰時間を減衰率に変換する方法

前回の記事で減衰時間から減衰率を求める計算式を考えましたが、K = e-0.01/T とかなり変な式、グラフを引いてみると、こんな風に極端な形で、この減衰率をマイクロプロセッサで求めるのはこれまた単純にはゆかなさそうです。

上記の式は非線形の関数なので、プロセッサで都度計算するのは時間がかかり、処理が追い付かなくなりそうです。近似値でよいので素早く計算したいです。手っ取り早いのはテーブル参照と補間・補外を使う方法、それなら実用的な速さで計算できそうです。

入力は、1024 ステップの整数値でこれを欲しい減衰時間に線形に対応させます。出力はその減衰時間のための減衰率です。テーブルには主要の入力値にたいする減衰率を覚えさせます。

ただここで一つ問題が。参照テーブルの値の精度をどれぐらい取る必要があるか試行錯誤しましたが、今の僕の技術ではどうしても 32bit 必要です。これを 1024 点全部覚えさせるには、1024 x 4 = 4KB の SRAM が必要です。今プロセッサに AVR の ATmega168 を使っていますが、内部 SRAM が 1KB しかありません。思い切りあふれています。テーブルの点数を 256 点にして補間しましたが変化が急激すぎて時間分解能が足りません。それでもまだメモリがあふれています。テーブルの作りに工夫が必要です。

その前に、まず時間変化の範囲をどれぐらいに取る必要があるか決めました。これは音を出して決めるのが一番、パラメータを変え、音を聞きを繰り返して、多分時定数を 1ms から 5s の範囲で動かせれば十分なようです。いずれもっと長い立ち上がりと減衰が欲しくなるかもしれませんが、まあ将来の改良点ということにしてほっときます。

で、もう一度上のグラフを見ると、カーブが急激に動くのはせいぜい 1 秒ぐらいまでで、あとはほぼ横ばいです。なので、時間の短いところと長いところでテーブルの時間分解能を変えるのはどうでしょうか?目標は 128点。これで 512B のサイズで、メモリには余裕ができます。

これまた試行錯誤して、最初の 64 点 (0.3秒ぐらいまで) は入力値とテーブルを一対一で対応して、残りの 1024 – 64 = 960 点は時間分解能がそれほど必要ないのでテーブルの後半のもう 64 点でカバーしてみると、かなり良い精度で近似できました。後半はテーブル一点で入力値の15点をまかなう計算になります。それだと値がカクカクするから線形に補完します。

テーブルの値の計算をするプログラムはこちら。ちょっと雑なプログラムですけどまあさっと動かしてしまいたかったので。テーブルの値は、Q0.32 の固定小数点数です。減衰率なので値は絶対に 1 未満なのでこれで大丈夫です。

#include <math.h>
#include <stdint.h>
#include <stdio.h>

uint32_t table[128];

int main(int argc, char *argv[]) {
  for (int i = 0; i < 1024; ++i) {
    if (i >= 64 && (i - 64) % 15 != 0) {
      continue;
    }
    double x = -0.0001 / (0.001 + 5.0 * i / 1024.0);
    uint32_t value = ((uint32_t) (exp(x) * (float) ((uint64_t) 1 << 32)));
    int index = (i < 64) ? i : ((i - 64) / 15 + 64);
    table[index] = value;
    printf("%12u,", table[index]);
    if (index % 4 == 3) {
      printf("\n");
    }
  }
  return 0;
}

このプログラムで作ったテーブルをマイクロプロセッサに埋め込んで減衰率を求める関数を書きました。今のところ良好です。

/*
 * Exponential decay/grow ratio lookup table.
 *
 * The table stores data points to give decay/grow ratio for lookup keys ranging
 * from 0 to 1023 (inclusive). The key 0 receives the ratio for time constant
 * 0.001s (1ms) approximately. The key 1023 receives the ratio for time constant
 * about 5s. The values in th table are Q0.32 unsigned fixed point numbers.
 *
 * In order to save data memory space, the values have different key resolutions
 * depending on the position in the table.
 * Table values with indexes from 0 to 63 has 1-1 key mapping to index,
 * i.e.
 *   ratio = kTransientRatios[key]
 * But for the indexes from 64 to 127, the table keeps only values for every 15 keys,
 * i.e.
 *   ratio = kTransientRatios[(key - 64) / 15 + 64]
 * where (key - 64) % 15 == 0
 *
 * Ratios for keys that are not on these indexes need to be calculated
 * by interpolation or extrapolation. See the method GetRatio().
 */
uint32_t kTransientRatios[128] = {
  3886247118,  4222575580,  4255256816,  4267608186,
  4274098987,  4278100538,  4280814394,  4282775974,
  4284259997,  4285421933,  4286356375,  4287124175,
  4287766262,  4288311174,  4288779419,  4289186113,
  4289542645,  4289857756,  4290138268,  4290389583,
  4290616034,  4290821137,  4291007774,  4291178333,
  4291334804,  4291478865,  4291611934,  4291735225,
  4291849776,  4291956486,  4292056132,  4292149393,
  4292236865,  4292319070,  4292396469,  4292469473,
  4292538445,  4292603710,  4292665560,  4292724255,
  4292780031,  4292833101,  4292883656,  4292931872,
  4292977906,  4293021905,  4293063999,  4293104310,
  4293142949,  4293180018,  4293215611,  4293249813,
  4293282706,  4293314362,  4293344850,  4293374234,
  4293402573,  4293429921,  4293456330,  4293481846,
  4293506515,  4293530379,  4293553475,  4293575840,
  4293597508,  4293856889,  4294033677,  4294161903,
  4294259161,  4294335461,  4294396917,  4294447478,
  4294489805,  4294525758,  4294556676,  4294583547,
  4294607117,  4294627960,  4294646522,  4294663159,
  4294678155,  4294691742,  4294704109,  4294715414,
  4294725787,  4294735340,  4294744166,  4294752345,
  4294759946,  4294767027,  4294773641,  4294779831,
  4294785639,  4294791097,  4294796237,  4294801085,
  4294805667,  4294810002,  4294814111,  4294818011,
  4294821717,  4294825243,  4294828603,  4294831807,
  4294834867,  4294837792,  4294840590,  4294843270,
  4294845839,  4294848303,  4294850670,  4294852944,
  4294855131,  4294857236,  4294859264,  4294861218,
  4294863103,  4294864922,  4294866678,  4294868376,
  4294870017,  4294871604,  4294873141,  4294874628,
  4294876070,  4294877467,  4294878823,  4294880138,
};

uint32_t GetRatio(uint16_t key) {
  const int kTableBoundary = 64;
  if (key < kTableBoundary) {
    return kTransientRatios[key];
  }
  const int kResolution = 15;
  int index = (key - kTableBoundary) / kResolution + kTableBoundary;
  int ramainder = (key - kTableBoundary) % kResolution;
  register uint32_t temp = kTransientRatios[index];
  if (ramainder > 0) {
    register uint32_t temp2;
    if (index == 127) {
      temp2 = temp - kTransientRatios[index - 1];
    } else {
      temp2 = kTransientRatios[index + 1] - temp;
    }
    temp2 *= ramainder;
    temp2 /= kResolution;
    temp += temp2;
  }
  return temp;
}

減衰係数はどうやって計算すれば良いか?

ちょっと間違っているかもしれませんが後で見直すためのメモでもあります。

まず放電とかの減衰曲線の時定数の定義が何だったか忘れました。調べたら行き着くのは当然 Wikipedia、もとの値の約 37% ということで、なんでそんな半端な数なのかも思い出せませんが、とにかく、37% ね。ちなみに時定数を英語で何というかも知りませんでした。普通に time constant なのですね。あっ英語ページには 1 / e ≒ 36.8% て書いてある。と思ってよく見たら日本語ページにも書いてあるじゃないか。

と、与太話はこれぐらいにして、今設計している EG モジュールの PWM は、約10ミリ秒ごとに値を更新するので、T 秒で 1 / e に落とすなら、T / 0.01 回の更新で値が 1 / e にまで減衰すれば良いのだから、

K = (1/e)0.01/T = e-0.01/T

プロットしてみました

なんというか、不安を誘う曲線ですね。精度とか大丈夫だろうか?横軸を対数に取ってみたらどうでしょうか?

うん、これならまだ何とかなるかもしれません。手持ちのボリュームに B カーブのものがなくてDカーブのものを使っているのがとっても気になりますが、もうハードウェアは作ってしまったので仕方がない。でもいずれにせよこの曲線をリアルタイムで計算するのはちょっと厳しそうです。実装にまた工夫がいります。先は長い、かも

マイクロプロセッサでの減衰曲線の実装

前回、参照テーブルを使わずに減衰曲線を生成する方法を考察しましたが、その方法をマイクロプロセッサで試してみました。やはり一筋縄ではいかないようです。

実際のプログラムは記事に載せるには長すぎますが、関係のある部分を抜粋して、こんな感じのプログラムを組みました。一見動きそうだし、これを普通に PC で実行すれば動作するはずですが、AVR ではうまく動作しません。

volatile uint32_t g_value;  // value holder, Q10.22 fixed point number
volatile uint32_t g_decay_factor;  // amount of decay,  Q0.32 fixed point
volatile uint8_t g_update_ready;

ISR(Time1_OVF_vect) {
  OCR1A = g_value >> 22;  // reflect the value to PWM、Q10.22 to integer
  g_update_ready = 1;
}

void UpdateValue() {
  if (!update_ready) {
    return;
  }
  // reduce the value by multiplying the decay factor
  uint64_t temp = g_value;
  temp *= g_decay_factor;
  g_value = temp >> 32;
  g_update_ready = false;
}

int main() {
  g_value = 1023 << 22; // integer to Q10.22
  g_decay_factor = 0.99309249543703590153321021688807 * 0x100000000;
  g_update_ready = 0;

  while (true) {
    UpdateValue();
  }
}

どうも AVR 用のコンパイラ avr-gcc は、C 言語標準に厳密には従っていないようです。どこかにドキュメントがあるかもしれませんが見つけていません。レジスタを細かく制御したいための仕様なのかもしれません。レジスタを意識して、データ領域からレジスタにどんな風にデータをロードしてレジスタ上でどんな演算をするのか意識して書くと動くようです。上記のプログラムは以下のように書き換えれば動くはずです。

volatile uint32_t g_value;  // value holder, Q10.22 fixed point number
volatile uint32_t g_decay_factor;  // amount of decay,  Q0.32 fixed point
volatile uint8_t g_update_ready;

ISR(Time1_OVF_vect) {
  register uint32_t temp = g_value;  // copy the Q10.22 value to a variable
  temp >>= 22;  // Q10.22 to integer
  OCR1A = temp;  // reflect the value to PWM
  g_update_ready = 1;
}

void UpdateValue() {
  if (!update_ready) {
    return;
  }
  register uint64_t temp = g_value;
  // copy the Q10.22 value to a 64-bit int
  temp *= g_decay_factor;
  // multiply by the Q0.32 decay factor
  temp >>= 32;  // shift 32 bit to adjust scale, now it's Q10.22 again
  g_value = temp;
  // copy back to the value holder
  g_update_ready = false;
}

int main() {
  register uint32_t temp = 1023 << 22;   // this is calculated by the compiler
  g_value = temp;
  g_decay_factor = 0.99309249543703590153321021688807 * 0x100000000;
  g_update_ready = 0;

  while (true) {
    UpdateValue();
  }
}

temp は register 変数として宣言しなくても動くようです。コンパイラが良きに計らってくれるようです。

要点は、ビットを伸ばしたり縮めたりする演算をするときには、以下の点を気を付けると動くようです

  • データをどこにロードさせるかはっきりさせる (一時変数を宣言する)
  • 演算がどこで行われるかはっきりさせる
    • 二項演算子 +, -, *, >>, << のかわりに +=, -=, *=, >>=, <<= を使う。
    • 一行一演算、一行に二つの演算を入れない

「ようです」が多いのが恐縮です。このことが書かれているドキュメントを読んだわけではないので、コンパイラの挙動から推測しています。

指数関数曲線を固定小数点演算で求める

エンベロープジェネレータのプログラムですが、現状は指数関数テーブルを参照して出発点の値に時間ごとに減衰する指数関数を掛け合わせて減衰曲線を作っています。

https://github.com/naokiiwakami/vceg

これは精度をあまり気にせずにざっと計算できるので手っ取り早いのですが、プログラムはどうしても煩雑になります。また、値の更新周期は約 10ms なのですが、減衰が遅いと 10ms の間には値が動かず、時間軸を補間しないといけません。これはまだ実装していませんがプログラムはもっと煩雑になります。

ちょっと嫌だなと考え始めたのですが、そもそも指数関数の減衰は

V(i+1) = V(i) * k

where 0 <= k < 1

のような式で簡単に作れそうなものなのですが、実はそう甘くはありません。例えば、1秒で半分に減衰する曲線を、10ms の更新周期で計算しようとすると、k の値は

k = 0.5 ^ (1 / 100) = 0.99309249543703590153321021688807

限りなく 1 に近いけれども 1 にしてしまったら減衰しません。変数として float を使えばいいじゃない?とも思いますが、例えば AVR で float を static 領域に持つとなんだかわからないけれどもデータがばんばん壊れてしまいます(バグかもしれませんが)。それにそもそもふつうのマイクロプロセッサでの浮動小数点演算は遅いです。ということで、固定小数点演算で曲線が求まらないか考えてみました。

ブランクがあまりにも長いのでやり方を忘れてしまいました。ので調べました。以下の解説がとてもわかりやすかったです

https://www.allaboutcircuits.com/technical-articles/multiplication-examples-using-the-fixed-point-representation/

いつか上記リンク先がなくなってしまうかもしれないので解説を書いとくのも良いのですが、何分もう夜分遅いので、書いたプログラムだけだっと貼っちゃいます。

#include <stdio.h>
#include <stdint.h>

int main(int argc, char *argv[]) {
  // fixed point integer 10bit.6bit
  uint16_t value = (1023) << 6;
  // fixed point integer .16bit
  uint16_t decay_factor = 0.99309249543703590153321021688807 * 65536;
  for (uint16_t i = 0; i < 1000; ++i) {
    // multiplication
    uint32_t temp = value * decay_factor;
    value = temp >> 16;
    // retrieve the integer part
    uint16_t output = value >> 6;
    printf("%d %d\n", i, output, 10);
  }
  return 0;
}

結果もさくっと貼っちゃいます

これはいけそうな気がします。

電圧制御充放電回路

電圧制御のエンベロープジェネレータを開発中ですが、今のところすべてマイクロプロセッサで計算して値を出力する予定。現実的ですがもうちょっとアナログでやれないのかとも考えました。以前から何度か挑戦しているんですが。アナログ式は難しいです。

前回検討した時には、簡易的な指数関数回路を使いましたが、あまりうまくゆきませんでした。確か、エンベロープの遷移の最短と最長の差があまり大きくならないのが最大の問題でした。どういうことなのか。シミュレータでちょっと調べてみました。いろいろとパラメータを調整してみましたが直線性と変化のレンジのバランスはこの辺が精いっぱいです。最小と最大で数百倍というところでしょうか。これでは全然足りず、充放電回路もシミュレーションしてみましたがやはり実用にならなさそうです。

ということで、もっときちっとした指数関数回路を使ってみましょう。

今度は直線性がずっと良いし制御電流のダイナミックレンジも1000倍を超えています。もっと大きなダイナミックレンジも取れましたが大きすぎるのもダメなのでの1000倍です。この回路ならいけそうな感じがします。ではこれを充放電回路と組み合わせてみましょう。制御するならマイクロプロセッサでやると思うので、制御電圧も充放電入力も 0 – 5V で設計しています。以下が制御電圧を最大 (5V) にした時の遷移

こちらは、制御電圧を中天にした場合 (2.5V)

こちらは最小にした場合です。

これなら使えそうな気がします。でもこれをマイクロプロセッサで制御するとなると結構面倒な回路になるしプログラミングはそれほど複雑にはならなさそうですが、PWM 出力を 2ch 占有するのはどうにもよろしくないような気もします。あとアタックからディケイへのフェーズ切り替えに比較機を外付けしないといけないし(プロセッサにも内蔵していますが AD 変換器と共有なのでタイミング制御にはちょっと使いたくない)、エンベロープをもっと複雑にしたかったらその比較機にも制御信号をぶち込まないといけないわけで、どうにも中途半端な感じがします。むしろアナログシーケンサみたいなものとつなげたほうがおもしろくならないかな?

ということで、回路を組むところまではやらないと思いますが検討にかなりの時間がかかったので記録に残しておきます。何か別のモジュールを思いついた時のために。

ウインドコントローラの進捗

遅々としていますが進んではいます。

大まかにファームウェアを書いてみたところで、どう作ろうか大体の方針が決まってきたので回路を決めて基板を作る準備をしています。

出力は二系統にして、VCF と VCA 別々に入れるようにしようと考えています。なので、パラメータ数はボリュームの数4 よりも多いので、操作性を工夫しないといけません。スイッチは一個、表示用の LED は二個。電源は Analog2.0 の接続方法をとっています。ユーロラックに入れることも考えているのですが Analog2.0 のライフラインバスをそのまま使うか悩ましいところです。今はとりあえず保留。製作時間が非常に限られているので手を止めないのがとっても重要です。

細かい動きはすべてファームウェアで決めることになります。そちらに時間がかかるのでハードウェアはとっとと作ってしまいたいです。