🐥

『Pythonではじめる音のプログラミング』をRustではじめてみたときのメモ

2023/01/05に公開

この本を買った。ちょうどいい機会なので、Rustの練習がてらやってみて、半分くらい読んで力尽きた。ここでは気付いたことをメモっておく。

https://www.ohmsha.co.jp/book/9784274228995/

コードはGitHubに置いたのでもし興味あれば。汚いですが...

https://github.com/yutannihilation/sound-programming-practice

Rustでやる場合に使うcrate

信号処理

本では、NumPyを使って信号の演算をしていく。Rustの場合、信号処理はずばりdaspというcrateがある。もっと低レベルから自分でやらないと勉強にならない?、という気もしつつ、Rustでやりたいと思った理由の半分はdaspに慣れることだったので、とりあえず使ってみる。

https://github.com/RustAudio/dasp

I/O

本では、wavファイル書き出しにSciPy(SciPyにそんな機能があるの知らなかった...)を使っている。Rustだと、wavファイルの読み書きはhoundというcrateが定番っぽい。

https://github.com/ruuda/hound

ただし、直接音を出すならcpalを使えばいいらしい。本がwav書き出しにした理由はおそらく、オーディオの入出力はややこしくてトラブルが多そう、というだけのような気がするので、直接音を出せるならそれでよさそう。とりあえずexampleをcargo runしてみたらトラブルなく動いたのでこっちでいく。

https://github.com/RustAudio/cpal

MIDI

MIDI関連はどれが決定版なのかよくわからなかった。daspとcpalを擁するRustAudioには、wmidirimdというcrateがあったりするけど、スターの数的にはこれがいちばんメジャー?

https://github.com/Boddlnagg/midir

ちなみに本では、自力でMIDIをパースしてるので、nomとかで自分で実装しないと勉強にならないのかもしれない。

第2章

サイン波を鳴らす

ちょうどいい例がdaspのexamplesにあるので、これを少し改変してサイン波を鳴らす。コードはここに置いた。

https://github.com/yutannihilation/sound-programming-practice/blob/master/examples/ch2-sine-wave.rs

いきなり長いけど順番に見ていく。

main()

let host = cpal::default_host();
let device = host.default_output_device().unwrap();
let config = device.default_output_config()?;

ここは、どういう処理なのかよくわからないけど、なんかいい感じにデフォルトのデバイスと設定の情報を取得してくれるらしい。便利。

match config.sample_format() {
    cpal::SampleFormat::F32 => run::<f32>(&device, &config.into())?,
    cpal::SampleFormat::I16 => run::<i16>(&device, &config.into())?,
    cpal::SampleFormat::U16 => run::<u16>(&device, &config.into())?,
}

ここもよく理解してないけど、デバイスによってサンプルのフォーマットが浮動小数点だったり整数だったりするらしいので、それに合わせた型を選ぶ必要があるらしい。

run()

run()は、以下のシグネチャになっている。

fn run<T>(device: &cpal::Device, config: &cpal::StreamConfig) -> Result<(), anyhow::Error>
where
    T: cpal::Sample,

cpal::Samplef32i16u16に対して定義されていて、要は from() ってやるとそのフォーマットに変換してくれる、フォーマットの違いを吸収するための trait らしい。

let sine = signal::rate(config.sample_rate.0 as f64)
    .const_hz(440.0)
    .sine();

let mut frames = sine.take(config.sample_rate.0 as usize);

sineSineという型のオブジェクトで、SineSignal traitを実装しているので信号を取り出せる。Signalには色々メソッドが生えていて、ディレイをかけたり周波数を変えたり増幅・減衰したりできるっぽい。sine自身はイテレーターではないが、next()で次のフレームに進んだり、take()で指定した数だけフレームを取ってきたりできる。ここでは、サンプルレートの数だけサンプルを取ってきている、つまり、サンプルレートとは1秒間のサンプルの数なので、1秒間分の信号が入ることになる。

そんなSineオブジェクトをつくるには、まずサンプルレートとベースの周波数の設定が必要になる。signal::rate().const_hz()がそれ。const_がないhz()だとFMがかけられるっぽい。

let (complete_tx, complete_rx) = mpsc::sync_channel::<()>(1);

let channels = config.channels as usize;
let stream = device.build_output_stream(
    config,
    move |data: &mut [T], _: &cpal::OutputCallbackInfo| {
        write_data(data, channels, &complete_tx, &mut frames);
    },
    |err| eprintln!("{err}"),
)?;

stream.play()?;

complete_rx.recv().unwrap();
stream.pause()?;

mpscが登場していきなり難しくなる...。ここは要は、cpal::Streamをつくって、さっきつくった信号を流し込んでいる。cpal::Streamは、自動では実行されない(こともある?)ので、明示的にplay()する必要がある。また、遅延とか出力の区切りを考えると?どこで終わるかわからないので、チャンネル越しに終了が知らされるのを待ってからpause()しないといけない、みたいな感じなんだと思う。あんまり理解できてないけど、とりあえずこの本の内容とは関係が薄そうなので、こういうものだと思って黙って写経しておく。

fn write_data<T>(
    output: &mut [T],
    channels: usize,
    complete_rx: &mpsc::SyncSender<()>,
    frames: &mut dyn Iterator<Item = f64>,
) where
    T: cpal::Sample,
{
    for frame in output.chunks_mut(channels) {
        let sample = match frames.next() {
            Some(sample) => sample.to_sample::<f32>(),
            None => {
                complete_rx.try_send(()).ok();
                0.0
            }
        };
        let value: T = cpal::Sample::from::<f32>(&sample);
        for sample in frame.iter_mut() {
            *sample = value;
        }
    }
}

実際にデータを書き込むところ。進むのは1フレームずつ進むが、出力のチャネルの数(たぶん多くの場合は左右の2つ?)だけ値をコピーする必要がある点に注意。
もしframeが終わりに達した(next()Noneを返した)ら、チャネル越しに終了を知らせて残りの出力は 0.0 で埋める。

クリックノイズ対策

上ので音は鳴るけど、サイン波の鳴り始めと鳴り終わりにクリックノイズが入ってしまう。本に書かれているように、音の前後にフェードイン、フェードアウトをつければクリックノイズが低減される。

これには、mul_amp()が使える。mul_amp()は、1つのSignalに別のSignalを掛け合わせる処理で、要はVCAみたいなものらしい。具体的には、sineからframeを取り出す前に.mul_amp(env)を付け加えることになる。

    let mut frames = sine.mul_amp(env).take(config.sample_rate.0 as usize);

で、ここのenvは、envelope generatorは用意されていないので自分で作るしかない。(注:後述のように、Iterator traitではなくSignal traitを直接実装した方がいいとあとで気付いた。ここは記録用に残しておく)しばらく悪戦苦闘してたけど、Iterator traitを実装しているstructがあれば、from_iter()Signalに変換できるのでそれほど難しくなかった。具体的には、以下のattack → sustain → releaseのenvelope generatorをつくった。アタックもリリースも直線というこのざっくりさ。

struct Env {
    cur_frame: usize,
    total_frames: usize,
    attack_frames: usize,
    release_frames: usize,
}

impl Env {
    fn new(total_frames: usize, attack_frames: usize, release_frames: usize) -> Self {
        Self {
            cur_frame: 0,
            total_frames,
            attack_frames,
            release_frames,
        }
    }
}

impl Iterator for Env {
    type Item = f64;

    fn next(&mut self) -> Option<Self::Item> {
        self.cur_frame += 1;

        // already ended
        if self.cur_frame > self.total_frames {
            return None;
        }

        // release phase
        if self.cur_frame > self.total_frames - self.release_frames {
            return Some((self.total_frames - self.cur_frame) as f64 / self.release_frames as f64);
        }

        // attack phase
        if self.cur_frame <= self.attack_frames {
            return Some(self.cur_frame as f64 / self.attack_frames as f64);
        }

        // sustain phase
        Some(1.0)
    }
}
const ATTACK: usize = 1000;
const RELEASE: usize = 1000;

// ...snip...

let env = signal::from_iter(Env::new(total_frames, ATTACK, RELEASE));

ただ、これだとまだ終わりのクリックノイズが取れない。本に書かれているもう一つの対策、無音のフレームを1000フレームほど足したら大丈夫だった。無音のフレームはsignal::equilibrium()で作れる。もしかして、鳴り終わりのタイミング間違えてるせい...?という気もするけどとりあえず気にしないことにする。

let mut frames = sine
    .mul_amp(env)
    .take(total_frames)
    // To prevent click noise at the end, fill some silence
    .chain(signal::equilibrium().take(1000));

第3章

色々試行錯誤したけど、どうやらSignal traitを実装するのがいいらしい。Signalは、Iteratorとほぼ同じだけど、Optionではなく必ず値が返る点が異なる。上に書いたようにIteratorfrom_iter()に渡してもSignalは作れるけど、ゼロから作るものなら初めからSignal traitを実装した方がいい。

8音分のエンベロープを生成する

最初は、envを渡しておいて、env.trigger()みたいなのであとからトリガーすればいい、と思ってたけど、cpalで出力するときに別スレッドになるのでちょっとややこしくなりそうで諦めた。とりあえず、boolの配列を受け取って、trueなら1拍分音を鳴らす、みたいなエンベロープを実装する。合ってるのかわからないけどこんな感じ?

struct Env {
    seq: Vec<bool>,
    cur_frame: usize,
    note_on: bool,
    step_length: usize,
    attack_frames: usize,
    release_frames: usize,
}

impl Env {
    fn new(
        seq: Vec<bool>,
        step_length: usize,
        attack_frames: usize,
        release_frames: usize,
    ) -> Self {
        Self {
            seq,
            cur_frame: 0,
            note_on: false,
            step_length,
            attack_frames,
            release_frames,
        }
    }
}

impl Signal for Env {
    type Frame = f64;

    fn next(&mut self) -> Self::Frame {
        self.cur_frame += 1;

        // proceed to the next step
        if self.cur_frame > self.step_length {
            self.cur_frame -= self.step_length;
            self.note_on = self.seq.pop().unwrap_or(false);
        }

        if !self.note_on {
            return 0.0;
        }

        // release phase
        if self.cur_frame > self.step_length - self.release_frames {
            return (self.step_length - self.cur_frame) as f64 / self.release_frames as f64;
        }

        // attack phase
        if self.cur_frame <= self.attack_frames {
            return self.cur_frame as f64 / self.attack_frames as f64;
        }

        // sustain phase
        1.0
    }
}
const SEQ: [bool; 8] = [true; 8];

// ...snip...

let step_length = config.sample_rate.0 as usize;
let env = Env::new(SEQ.to_vec(), step_length, ATTACK, RELEASE);

8音を順番に鳴らす

2章の例では一定の音だったのでconst_hz()を使った。ここでは、周波数を途中で変えられるhz()を使う。といっても、hz()は、あとで外から周波数を変えることはできない。hz()は引数にSignalを渡す必要があり、そのSignalの値に従って周波数が設定される。なので、まずは1拍ごとに指定した周波数の値を返すSignalをつくる必要がある。

さっきとだいたい同じでこんな感じ。

struct Track {
    seq: Vec<f64>,
    step_length: usize,
    cur_frame: usize,
    note: f64,
}

impl Track {
    fn new(seq: Vec<f64>, step_length: usize) -> Self {
        Self {
            seq,
            step_length,
            cur_frame: 0,
            note: 0.0,
        }
    }
}

impl Signal for Track {
    type Frame = f64;

    fn next(&mut self) -> Self::Frame {
        self.cur_frame += 1;

        // proceed to the next step
        if self.cur_frame > self.step_length {
            self.cur_frame -= self.step_length;
            self.note = self.seq.pop().unwrap_or(0.0);
            println!("note: {}", self.note);
        }

        self.note
    }
}
const TRACK1: [f64; 8] = [659.26, 587.33, 523.25, 493.88, 440.00, 392.00, 440.00, 493.88];

let track1 = signal::rate(config.sample_rate.0 as f64)
    .hz(Track::new(TRACK1.to_vec(), step_length))
    .sine();

2音を同時に鳴らす

2つのSignalを加算するにはadd_amp()を使うといいらしい。エンベロープを掛けるのとあわせて、こんな感じ。

let mut frames = track1
    .add_amp(track2)
    .mul_amp(env)
    .take(step_length * SEQ.len())
    .chain(signal::equilibrium().take(1000));

第5章

Biquad filter

ここは、コードがやってることの説明がなくてちょっと不親切な気がした。個人的には、サンプリング周波数で中心周波数を割るあたりとかがさっぱりわからず、ググって出てきたこちらのサイトの方を参考に実装した。

https://webaudio.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html

方針としては、Signalを受け取ってSignalを返すstructを実装する。ring bufferが用意されているので無駄にそれも使ってみる。

struct Lpf<S: Signal<Frame = f64>> {
    signal: S,
    fs: f64, // sampling rate
    fc: f64,
    q: f64,
    before: dasp::ring_buffer::Fixed<[f64; 2]>,
    after: dasp::ring_buffer::Fixed<[f64; 2]>,
}

impl<S: Signal<Frame = f64>> Lpf<S> {
    fn new(signal: S, fs: f64, fc: f64, q: f64) -> Self {
        println!("central frequency: {fc}");
        println!("Q: {q}");

        Self {
            signal,
            fs,
            fc,
            q,
            before: dasp::ring_buffer::Fixed::from([0.0; 2]),
            after: dasp::ring_buffer::Fixed::from([0.0; 2]),
        }
    }
}

impl<S: Signal<Frame = f64>> Signal for Lpf<S> {
    type Frame = f64;

    // c.f. https://webaudio.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html
    fn next(&mut self) -> Self::Frame {
        let orig = self.signal.next();

        let pi = std::f64::consts::PI as Self::Frame;
        let omega0 = 2.0 * pi * self.fc / self.fs;
        let alpha = omega0.sin() / 2.0 / self.q;

        // Since `push()` pushes on to the back of the queue,
        //
        //   - before[1]: input of 1-step before
        //   - before[0]: input of 2-step before
        //   - after[1]:  output of 1-step before
        //   - after[0]:  output of 2-step before
        //
        let mut out = (1.0 - omega0.cos()) / 2.0 * orig
            + (1.0 - omega0.cos()) * self.before[1]
            + (1.0 - omega0.cos()) / 2.0 * self.before[0]
            - (-2.0 * omega0.cos()) * self.after[1]
            - (1.0 - alpha) * self.after[0];
        out /= 1.0 + alpha;

        self.before.push(orig);
        self.after.push(out);

        out
    }
}

ちなみに、Biquad filterは、daspにも実装の提案がされているが、どれくらいの抽象度のものを入れるかで議論があり(あるいは単にメンテナの時間がないので?)マージされていないらしい。まあ、実装してみて思ったけど、こういうエフェクト部分ってどういう抽象度でどこに置くかあんまり自明じゃないですよね。

https://github.com/RustAudio/dasp/pull/148

第6章

Poly BLEP

まあ、この辺りはきっとdasp側で実装されてるしやることないでしょ、と思ってスキップしようとしたけど、コードを覗いてみるとそういう考慮はなにもなかった...。

https://github.com/RustAudio/dasp/blob/6b15274b471835e586089e54228e54601f92d391/dasp_signal/src/lib.rs#L1756-L1767

daspの実装は、実際に波形を生成するSawSquare[0,1] の位相だけを知ればいいようになっている。この位相はPhaseから渡される。この実装の仕方だと、Poly BLEPにはあまり都合がよくない。Poly BLEPではステップ間の時間を知る必要があるが、ステップの長さや周波数を直接知ることはできない。Phaseじゃなくて周波数もパラメーターとして取るやつをつくるか迷ったけど、前回の位相を持つようにして毎ステップ直前からの位相差を計算することにした。効率は悪そう。

とりあえず、のこぎり波だけやってみる。

pub struct PolyBlepSaw<S> {
    phase: Phase<S>,
    prev_phase: f64,
}

impl<S: Step> PolyBlepSaw<S> {
    fn new(phase: Phase<S>) -> Self {
        Self {
            phase,
            // TODO: The initial phase is not always 0.0?
            prev_phase: 0.0,
        }
    }
}

impl<S: Step> Signal for PolyBlepSaw<S> {
    type Frame = f64;

    fn next(&mut self) -> Self::Frame {
        let phase = self.phase.next_phase();
        let mut out = phase * -2.0 + 1.0;

        let delta = if phase > self.prev_phase {
            phase - self.prev_phase
        } else {
            // if the phase decreased, it should be because the phase got wrapped at 1.0.
            1.0 + phase - self.prev_phase
        };

        if phase < delta {
            let t = phase / delta;
            out += -t * t + 2.0 * t - 1.0;
        } else if phase > 1.0 - delta {
            let t = (phase - 1.0) / delta;
            out += t * t + 2.0 * t + 1.0;
        }

        self.prev_phase = phase;

        out
    }
}

この実装はDaisySPの実装を参考にした。ちなみに、DaisySPだと最後-1をかけてるけど、これはなぜなのかわからなかった。。

https://github.com/electro-smith/DaisySP/blob/b9f9048e064b80ddca8d6e81dc94d035830d190b/Source/Synthesis/oscillator.cpp#L34-L39
https://github.com/electro-smith/DaisySP/blob/b9f9048e064b80ddca8d6e81dc94d035830d190b/Source/Synthesis/oscillator.cpp#L69-L86

FM

FMは、const_hz()の代わりにhz()を使えばいい。そしてhz()の引数にモジュレーターのSignalを渡せば、その値に従ってキャリアの周波数が変化する。問題は、通常のSignal[-1,1] の範囲になっていてFMに使うのに適切なレンジになっていないこと。これは、スケールするにはscale_amp()、オフセットをかけるにはoffset_amp()で合わせればいい。

ちなみに、ここでさっきのPoly BLEPのありがたみが分かった。素のsawだと音が汚すぎる...

let base_hz = 440.0 * 8.0;
let ratio = 3.5;
let depth = 400.0;

let modulator = PolyBlepSaw::new(
    signal::rate(config.sample_rate.0 as f64)
        .const_hz(base_hz * ratio)
        .phase(),
    )
    .scale_amp(base_hz)
    .offset_amp(depth);

let carrier = PolyBlepSaw::new(
    signal::rate(config.sample_rate.0 as f64)
        .hz(modulator)
        .phase(),
);

Karplus-Strong

ここは、本文の説明は悪くない気がするけど、平均を取って差し引く、みたいな本文にない処理がコードにあって混乱した。どうやらこれ?と思ったけど、別にlow-pass filterでもない気がするし、理解できなかった。DCキャンセル? なくても音は鳴ったのでとりあえず無視する。

This is a very energetic excitation, and usually in practice the white noise is lowpass filtered;
(https://www.dsprelated.com/freebooks/pasp/Karplus_Strong_Algorithm.html)

あと、daspで実装するうえでよくわからなかったのが、delay用のバッファの作り方。本文で論じられているようにdelay timeがサンプル長の整数倍になるとは限らないので、隣接するサンプル2つを補間して値を返す必要がある。DaisySPもエルミート補間して返す実装になっていた。

https://github.com/electro-smith/DaisySP/blob/6ad2e5928b1b0af924ae7e6cd23d08967acb8864/Source/Utility/delayline.h#L87-L104

daspのring bufferはただのFIFOで値を1つしか返してくれないので、こういうのやろうと思ったら自分で実装するしかない...? よくわからなかったので、今回はとりあえずdelay timeは固定なので、前ステップのサンプルも別で保存してそれを使った。定数の計算部分は省略して、こんな感じになった。あってるかあまり自信ない...

impl Signal for KarplusStrong {
    type Frame = f64;

    // c.f. https://webaudio.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html
    fn next(&mut self) -> Self::Frame {
        self.cur_frame += 1;

        let cur_delayed_sample = self.delay_line.pop().unwrap_or(0.0);

        let all_passed_feedback = -self.g * self.last_all_passed_feedback
            + self.g * cur_delayed_sample
            + self.last_delayed_sample;

        // trigger once per second with the same lenght as the delay line
        let orig_noise = if self.cur_frame % (self.fs as usize) < self.delay_line_length {
            self.noise_source.next_sample()
        } else {
            0.0
        };

        let out = orig_noise
            + self.c
                * ((1.0 - self.d) * all_passed_feedback + self.d * self.last_all_passed_feedback);

        self.last_all_passed_feedback = all_passed_feedback;
        self.last_delayed_sample = cur_delayed_sample;
        self.delay_line.push(out);

        out
    }
}

とりあえず音は鳴った。

https://twitter.com/yutannihilation/status/1612270915974463488

感想

この本は、図がふんだんにあって悪い本ではないと思うけど、とにかくコードにコメントがないのがつらかった。最初に一気に説明して、最後にコードをドカッと載せる、というスタイルだけど、どの説明がどのコードに対応してるのかが分かりづらい。まあでも、苦しい分勉強になったといえば勉強になった。

あと、Rustでやってみた感想としては、daspなんかしっくりこない...。APIはきれいにまとまってると思うんだけど、隠蔽され過ぎてるというか、例えば、マウスでクリックしてエンベロープをトリガーするみたいなリアルタイムの操作どうやるの?、みたいな疑問が無限に湧いてくる。まあ、ここは自分のRust力が低すぎるだけ、という可能性もある。

Discussion