信号の包絡線を計算する (Common Lisp)

信号の包絡線 (signal envelope) は「信号の頂点にシーツをかぶせたときのような概形」を持つ波形です。

MatlabRで書いた信号包絡線の計算プログラムCommon Lispに移植しました。信号包絡線の計算方法は振幅の絶対値に対して、(1)低域フィルタをかける、(1')Savitzky-Golayフィルターがよい、(2)ピーク検出をしてスプラインでつなぐ、(2')いやいや秋間補間がよい、(3)ヒルベルト変換を使う、(3')結局それに低域フィルタかけるじゃろ、(4)スペクトル包絡には線形予測符号も使われてるよね、などなどいくつも提案できるのですが、ここではヒルベルト変換を用いて解析信号を計算する方法でいきます。解析信号の絶対値を計算すると信号包絡線になります。

まずは解析信号の実装。フーリエ変換をして、負の周波数成分をゼロにして、逆フーリエ変換します。

(defun hilbert (x)
  "Compute analytic signal of real signal vector X using Hilbert transform."
  (let* ((xx (napa-fft:fft (zeropad-pow2 x)))
         (yy (make-array (length xx) :initial-element #C(0 0)))
         (zz nil)
         (out (make-array (length x))))
    (setf (aref yy 0) (aref xx 0))
    (loop for n from 1 below (/ (length xx) 2)
          do (setf (aref yy n) (* 2 (aref xx n))))
    (setf (aref yy (/ (length xx) 2)) (aref xx (/ (length xx) 2)))
    (loop for n from (1+ (/ (length xx) 2)) below (length xx)
          do (setf (aref yy n) 0))
    (setq zz (napa-fft:ifft yy))
    (loop for n from 0 below (length out)
          do (setf (aref out n) (aref zz n)))
    out))

Matlabhilbert()と比較してみましょう。

>> hilbert([1 2 -3 4])
ans =
   1.0000 + 1.0000i   2.0000 + 2.0000i  -3.0000 - 1.0000i   4.0000 - 2.0000i
CL-USER> (hilbert #(1 2 -3 4))
#(#C(1.0d0 1.0d0) #C(2.0d0 2.0d0) #C(-3.0d0 -1.0d0) #C(4.0d0 -2.0d0))

同じ結果になっていますね。ヨシ。

次に信号包絡線を計算します。

(defun envelop (x)
  "Calculate amplitude envelop of a signal vector X."
  (let* ((xh (hilbert x))
         (xe (make-array (length xh))))
    (dotimes (n (length xh))
      (setf (aref xe n) (abs (aref xh n))))
    xe))

こちらもMatlabと比較してみましょう。

>> abs(hilbert([1 2 -3 4]))
ans =
    1.4142    2.8284    3.1623    4.4721
CL-USER> (envelope #(1 2 -3 4))
#(1.4142135623730951d0 2.8284271247461903d0 3.1622776601683795d0
  4.47213595499958d0)

本物の音響信号でやってみましょう。スネアドラムを一発叩いた音の録音があるので、信号包絡線を計算してみます。wavreadは以前作ったプログラムです。(→Common Lisp用のwavread - 丸井綜研

(ql:quickload :clgplot)
(let* ((snd (wavread "snare04-short.wav"))
       (x (first snd))
       (fs (second snd))
       (xe (envelop x))
       (xel (filt (biquad-lpf 100 (sqrt 1/2)) xe))
       (tt (list->vector (mapcar #'(lambda (v) (/ v fs)) (range (length xe))))))
  (clgp:plot x :x-seq tt)
  (clgp:plot xel :x-seq tt))

結果は以下のようになりました。左側が元の信号で右側が包絡線です。横軸は時間(秒)で縦軸は振幅(プラスマイナス1の範囲に正規化)。少しスムーズにするために包絡線にカットオフ周波数100 Hzの低域フィルタをかけてあげています。

Aseprite買った

以前から欲しいと思っていたAsepriteというドット絵作成ツールが半額になっていたので、思わず購入してしまった。半額セールはたぶん2月13日まで。

www.aseprite.org

で、初めて描いてみたのが以下の絵。16×16ピクセルBrüel & Kjær Type 4128-Cさんを想像で描きました。

本物の写真を見てみたら、チョーカーみたいなのつけてるし、目はなかった。ドット絵のほうは緑色のゾンビか、サボテンみたい。

昔ね、素焼きの陶器でできた人形の頭かなんかにカイワレの種が埋まってて水をかけるとニョキニョキ育つというおもちゃがあったのよ。カイワレの芽が髪の毛みたいに見えて、ヘアカットするとカイワレが収穫できるというもの。それを思い出してしまいましたよ。(後日、調べてみたらChia PetとかChia Headというものだったらしい)

www.youtube.com

本日の給油(スーパーカブ110プロ/JA07)

給油日 オドメーター (km) 給油量 (L) 単価 (円/L) 燃費 (km/L) 距離単価 (円/km)
2022-07-27 13405.0 2.21 157.92 57.29 2.76
2022-08-31 13546.7 2.92 156.16 48.53 3.22
2022-11-05 13699.4 3.21 157.01 47.57 3.30
2023-02-06 13848.1 3.26 157.06 45.61 3.44
2023-03-02 13973.9 2.60 156.92 48.38 3.24
2023-04-01 14095.9 2.99 156.86 40.80 3.84
2023-07-13 14218.0 2.73 163.00 44.73 3.64
2023-09-25 14380.5 3.60 176.94 45.14 3.92
2023-11-25 14529.5 2.93 162.12 50.85 3.19
2024-02-07 14644.2 2.80 165.00 40.96 4.03

ストップ・アンド・ゴーばかりの短距離運用だったので、燃費は過去最悪に近いです。寒い時期は燃費が悪くなりがちなのはしょうがないですね。

本日の給油(CRF1100Lアフリカツイン/2BL-SD10)

給油日 オドメーター (km) 給油量 (L) 単価 (円/L) 燃費 (km/L) 距離単価 (円/km)
2022-08-27 3401 16.63 157 16.42 9.56
2022-11-03 3727 17.20 158 18.95 8.34
2023-02-04 3987 16.75 157 15.52 10.12
2023-04-20 4165 10.37 157 17.16 9.15
2023-04-30 4563 18.22 160 21.84 7.32
2023-06-17 4982 16.41 159 25.53 6.23
2023-09-02 5296 19.46 178 16.14 11.03
2023-10-01 5700 18.89 180 21.39 8.42
2024-01-28 6003 18.62 165 16.27 10.14

前回アフリカツインに火を入れたのが11月23日のデイキャンプだったので、じつに2か月以上も放置してしまいました。前に乗っていた某大型バイクだと1か月エンジンをかけないとバッテリーが上がるという謎があったのでヒヤヒヤしましたが、さすがは安定のホンダ車。一発始動です。高価なリチウムイオン電池を積んでいると聞いたのと、交換が難しい箇所にあるとのことなので、工賃もかかるバッテリー交換はしたくない……(パーツカタログで調べたらエリーパワーのHY110https://amzn.to/49r3e6H)だそうです。5万円近くしますね……)

前回の単価は180円でしたが、会員価格で今回は165円。

今日は東京都~埼玉県を走り回り、ようやく6000 km超えました。

キャンプ道具を揃えた2023(更新版)

キャンプ道具を揃えたシリーズの2023年更新版です。これまで購入したものは以下。

そして、2023年に購入したものたち。

amzn.toで始まるリンクはアフィリエイト、それ以外のリンクはアフィではありません。避けたい人は気を付けてください。

寝袋

今年に入るまでは寝袋は家族のものを借りていましたが、キャンプに行くようになると自分用が欲しくなります。モンベルの「ホローバッグ#3(https://amzn.to/43GmlGG)」を購入しました。ファミリーバッグ#3と迷ったのですが、マミー型っぽくも使えるということでホローバッグのほうを選びました。使用可能限界温度が摂氏2度ということで、暑がりの僕なら温度がマイナスに行かない範囲の冬キャンプにも使えそうかもと思いましたが、真冬に使うには寒いと思いますよ……。

もうホローバッグはキッズ用しかカタログに載っていないので、ちょうど終売のタイミングで購入できたようです(アウトレット品はまだ入手可能かもしれません)。

webshop.montbell.jp

バッグ

雨の中のバイク移動に備えて、キャンプ道具を揃えた2022(バッグと小物類)でも書いていた防水バッグを買いました。トロリーシートバッグ DH-746https://amzn.to/3RCGvME)の上に60 Lの防水シートバッグを載せてもアフリカツインなら楽々です。モールシステムなどは何もついていないシンプルなものですが、容量は大きくて良いです。できるだけ荷物を減らしたいのですが、だんだんと増えてしまいますね。

www.daytona.co.jp

このバッグを購入したこともあって、アフリカツインのパニアケースにも食指が動いたのでした。

marui.hatenablog.com

飯盒のための色々

僕はとくにキャンプ地で何か特別な活動をしたい人ではありません。屋外炊爨してぼーっと過ごすということがしたいので、調理器具についつい手が伸びてしまいます。できるだけコンパクトに収納できる高機能なものを探していたのですが、戦闘飯盒2型が良さそうに思えました。すでにフェニックスライズの戦闘鉄板チタン FB-2 Tihttps://amzn.to/3L4u4H1)は持っていましたし、飯盒にあわせた大きさのFUZZチタン Ti-02https://amzn.to/3DdlRyh)もCRUNCH VOXhttps://amzn.to/3Ozf9YR)も揃えていましたので、実のところソラマメ型の飯盒を導入するのはほぼ確定していたのでした。

camphack.nap-camp.com

ソラマメ型の2合炊き飯盒は、ロスコ、エバニュー、ムースルームワークスなど、各社から似たようなものがいくつか販売されています。でも、買うからには自衛隊の仕様書(防衛省仕様書改正票「飯ごう,2形」DSP S 2036C(1))に完全準拠したものがいいです。ロスコの「戦闘飯盒2型」とエバニューの「山岳飯盒弐型」は取っ手の端の形状が仕様と異なりますし、ムースルームワークスの「戦闘飯盒2型」は耳金の塗装がされていません。

いろいろと探していたら、東海理研株式会社の製品情報陸上自衛隊仕様の携帯シャベルと飯ごうの記載があるのを見つけました。一般販売はされていなさそうなのですが、岐阜県関市へふるさと納税すると返礼品として飯盒2型がもらえるという情報を得ました。ふたに「GIFU SEKI」と書かれていてダサい郷土愛を感じるんですが、品物じたいは実際に自衛隊に納入されているものと同等に見えます。すでに関市の返礼品としての扱いは終わっていますが、群馬県太田市が同等品を扱っています。

www.furusato-tax.jp

実際に、この飯盒を入手してからは他の調理なべを使わなくなってしまいました。鍋もフライパンも湯沸かしもこれひとつですべてまかなえます。

アルコールストーブとライター

Amazonで半値以下の販売価格になっていたので、トランギアのアルコールバーナーhttps://amzn.to/3RjjmBm)を衝動買いしてしまいました。真鍮製なので、アルミと接触させないように気を付けないといけません。真鍮は銅と亜鉛の合金ですが、銅とアルミを接触させるとアルミに腐食が起こるらしいです。それは真鍮とアルミでも同じなのだとか。チタン製のアルストを狙っていたのですが、衝動買いは事故のようなものなので、しょうがないですね。

www.iwatani-primus.co.jp

フェニックスライズのCRUNCH VOXに突っ込んで飯盒炊爨に使ったり、ミニ五徳をのせて湯を沸かしたりしています。

火をつける道具はイワタニのアウトドアトーチバーナーとマッチしか持っていませんでした。イワタニのコンパクトジュニアバーナーには着火装置がついているので、火をつける専用品は不要だったのですね。一方、アルコールバーナーには着火装置がついていないので、別途SOTOのスライドガストーチ ST-487https://amzn.to/3PdRflJ)を購入しました。ライター用のガスを注入して使えます。(てっきりCB管から直接ガスを注入できるのかと思っていましたが、差し込み穴の径が異なるようでガスが入りませんでした)

ランタン

小さくて光量があると評判のゴールゼロの価格が落ち着いてきたので導入しました。いざというときに乾電池も使えるレッドレンザー ML4https://amzn.to/43w1K7W)とも迷ったのですが、最終的には充電池の容量が決め手になりました。レッドレンザーの750 mAhに対してゴールゼロは2600 mAhと3倍以上の差があります*1。容量が大きいほど長時間使えるわけなので、ゴールゼロを選びました。一泊二日のキャンプで電池切れになったことはありません。

スマホにもフラッシュ兼用のライトがついていますが、光量の面では専用品のほうが断然良いです。小さくてどこにでも持ち運べるので、キャンプ以外の場面でも重宝しています。

arinomi.co.jp

ペグとハンマーとグローブ

房総と猪苗代でのキャンプで困ったのが、テント付属のアルミペグが地面に刺さらないという問題でした。どちらも友人と一緒に行っていたので鍛造ペグを数本借りて急場をしのぎましたが、ソロのときにペグが刺さらないと危険です。テントが飛ばされないように、最低限の本数でいいので強いペグが欲しくなりました。また、小型のハンマーではうまくいかないことが明らかになったので、ハンマーも専用設計のものを買うことにしました。

チタン好きなので、ペグについてもチタン一択です。ハンマーは鍛造ペグで有名なエリッゼステークのものにしました。バイクで行くのでショート版。ステンレスヘッドのやつが外装に傷があるとかで安くなっていたので、貯まっていたAmazonポイントで購入。

両方そろうと使い勝手の差は明らかでした。ランタンもそうですが、とくにハンマーはなぜ最初からキャンプ専用品にしなかったのか……。

ペグうちなどの作業をするときにはモンベルの軍手を使っていましたが、いつのまにか失くしてしまったので、手袋専門メーカーが出しているイカツイやつにしてみました。熱いものも触れます。最初のうちは革のにおいが気になるので、陰干ししておくといいです。

*1:iPhone SE 3のバッテリー容量は2018 mAhなので、小型スマホよりも容量が大きいことになります。

ピンクノイズを生成するLISP関数

以前PythonやJSFX用のものも作ったのですが、ピンクノイズはいろいろなところで活躍するので、Common Lisp用にも作っておきます。

Voss-McCartney法

DSP generation of Pink (1/f) Noiseには、(1)ホワイトノイズにフィルタをかけてピンクにする方法、(2)Voss-McCartney法、の二つがあると紹介されています。Voss-McCartney法はホワイトノイズを低域に重みをつけて重ね合わせていく方法です。

具体的には、いくつかの乱数生成器を用意して、それらで1, 2, 4, 8, 16, ...サンプルごとに新しい乱数を作らせ、すべてを加算して出力するというものです。たとえば4サンプルごとの生成をする乱数生成器は、4サンプルに一度だけ乱数を更新しますが、それ以外の3サンプルではその前に生成した乱数を返します。より低い頻度(ここでは1/4の標本化周波数)での乱数生成をしていることになりますね。つまり、1, 2, 4, 8, 16, ...サンプルごと生成して加算するということは、オクターブ低いホワイトノイズを次々と作って足し合わせてピンクノイズを作っていくことになります。2nが大きくなればなるほど低い周波数まで伸びるピンクノイズになります。

Common Lispでの実装

とりあえず標本化周波数を決めておきます。

(defparameter +samprate+ 48000)

Voss-McCartney法の本体です。指定した秒数ぶんのピンクノイズを生成します。返り値ですが、多値を使ってピンクノイズ振幅値の列と各サンプルの時刻の列の二つを返しています。どちらも1次元配列(simple-vector)です。

(defun synth-pinknoise (dur &optional (fs +samprate+))
  "Generates pinknoise (wideband noise having -3dB/oct slope) of specified duration (seconds). Uses Voss-McCartney algorithm. http://www.firstpr.com.au/dsp/pink-noise/

Example:
(wavwrite \"pinknoise.wav\" (synth-pinknoise 2.0))"
  (let* ((lowfreq 20)
         (levels (ceiling (log (/ fs lowfreq) 2)))
         (tt (sample-time dur fs))
         (out (make-array (length tt)))
         (x (make-randn-array levels)))
    (dotimes (n (length tt))
      (dotimes (m levels)
        (if (= (1- (expt 2 m)) (mod n (expt 2 (1+ m))))
            (setf (aref x m) (randn)))
        (setf (aref out n) (+ (randn) (vecsum x)))))
    (values (normalize out) tt)))

sample-timeは、指定した秒数の各サンプルの時刻を計算する関数です。横軸に時間・縦軸に振幅をとってグラフを書く時などに使います。標本化周波数が48000 Hzのときに1秒間の時刻列は、#(0 1/48000 2/48000 ... 47999/48000)というようになります。以下のように定義しています。rangePythonの、linspaceMatlabの関数をそれぞれ真似しています。

(defun sample-time (dur &optional (fs +samprate+))
  "Make a vector of time values of each sample in seconds for specified duration."
  (coerce (mapcar #'(lambda (x) (/ x fs))
                  (range (round (* dur fs)))) 'vector))

(defun range (n)
  "Generates a list of integers from 0 to N-1."
  (linspace 0 (1- n) n))

(defun linspace (x1 x2 &optional (n 100))
  "Generates a list of N (100 by default) equally spaced numbers between X1 and X2."
  (let ((x (loop for x from 0 to (1- n)
                 collect (/ x (1- n)))))
    (mapcar #'(lambda (v) (+ v x1))
            (mapcar #'(lambda (v) (* v (- x2 x1)))
                    x))))

randnmake-randn-arrayは標準正規分布に従う乱数を作る関数で、以下のエントリに書いたものです。

marui.hatenablog.com

vecsumはベクトル内の数の総和を計算します。

(defun vecsum (vec)
  (loop as x across vec
        sum x))

ピンクノイズは音波形として使用したいので、出力時には±1の範囲に収めたいです。これをノーマライズといい、以下のnormalizeで実行しています。±1ビチビチに入れてしまうと(0 dBFSにすると)、オーディオ編集ソフトによっては「レベルオーバーしてるよ」という表示になることがあるので-1 dBFSをデフォルトにしていますが、-3 dBFSとかのほうが安全かもしれません。

(defun normalize (x &optional (dB -1.0))
  "Makes a sampled signal normalized to specified level (in dB)."
  (let ((res (make-array (array-dimensions x)
                         :initial-element 0.0))
        (gain (db->amp dB))
        (tmpmax (loop for n from 0 to (1- (length x))
                      maximize (abs (aref x n)))))
    (dotimes (n (length x))
      (setf (aref res n) (* gain (/ (aref x n) tmpmax))))
    res))

(defun db->amp (dB &optional (p0 1.0))
  (* (expt 10 (/ dB 20)) p0))

できあがったピンクノイズをオーディオ編集ソフトに取り込んで30分にまとめたのが以下の動画です。YouTubeにアップロードするときに聴覚符号化で圧縮されてしまうので、本来の音とは少し違いますが……。

youtu.be

同じスペクトルを持つのにクレストファクタが違う波形(Julia Advent Calendar 2023)

Julia Advent Calendar 2023の22日目の記事です。毎年Juliaで音に関係したことを書いていて、これまでにこのブログで書いたアドカレ記事は以下の通りです。

今年は、以下の論文を参考にして「同じスペクトルを持つのにクレストファクタが違う波形」を作ってみます。

  • Schroeder, M. R.. “Synthesis of Low-Peak-Factor Signals and Binary Sequences with Low Autocorrelation.“ IEEE Transactions on Information Theory, pp.85–89. January 1970. DOI: 10.1109/TIT.1970.1054411.

クレストファクタ

クレストファクタ(波高率)は「振幅と実効値の比」で計算されます。振幅の幅はゼロから波形ピーク値で測られます(文献によっては負側のピークから正側のピークまでとしているものもあります)。一方で実効値は波形の時間平均に近いもので、以下の式 (root mean square) で計算できます。音の波形のような交流波において振幅値の算術平均を計算すると正負の値が打ち消しあってしまうので、それを防ぐために二乗してから平均をとり平方根でもとの単位に戻すという流れで計算しています。

 \displaystyle
x_\mathrm{rms} = \sqrt{\frac{1}{N}\sum_{i=1}^{N}\left(x_i\right)^2}

もし波形が±1の範囲に正規化された正弦波(サイン波)だったとすると、振幅値のピークは1、実効値は 1/√2 になります。クレストファクタは振幅と実効値の比なので 1 ÷ 1/√2 = √2 です。正規化波形に対しての実効値は必ず1以下になるので、クレストファクタは1以上の値になります。最もクレストファクタが小さくなる(1になる)のは矩形波のように-1と+1の信号しか含まない信号です。

一般に、定常的な持続音(オルガンやクラリネットなど)のほうがクレストファクタは小さくなりがちで、逆に、大きな振幅値が時々しか現れない音(太鼓のリズムとか)はクレストファクタが大きくなりがちです。正規化波形で考えるので実際に音を聞いたときに感じる音の大きさとは異なりますが、クレストファクタは音ファイルから音量感を計算する指標の一つとして使うこともできます(たとえば以下の論文参照)。

パルス列の生成

例として余弦波の重ね合わせでパルス列を作っていくことを考えます。

 \displaystyle
r(t) = \sum_{k=1}^{N}\left(\frac{p_{k}}{2}\right)^{1/2} \cos\left(\frac{2\pi k t}{T} + \theta_k\right)

このとき、各余弦波の位相をあわせて(つまり 𝜃ₖ をすべて同じ値にして)重ね合わせると、高いクレストファクタ(論文中ではpeak-factorと書いてある)になります。ここで、𝑝ₖ は 𝑘 番目の倍音のパワーで、∑ 𝑝ₖ = 1を満たす必要があります。以下のコードでは論文中の数式に従って、第16倍音まで計算しています。

# number of harmonics
N = 16

# relative power of each harmonics
p = zeros(N)
for n in 1:N
    p[n] = 1/8 * (sin(π * (2n - 1) / 32))^2
end
p ./= sum(p)

# phase angles
ϕ = zeros(N)
for n in 1:N
    ϕ[n] = -π/2
end

# signal generation (Eq.1)
fs = 1000
T = 1
t = range(-π/2, stop=π/2, step=1/fs)
r = zeros(size(t))
for k in 1:N
    r .+= (p[k] / 2)^(1/2) * cos.(2π * k * t / T .+ ϕ[k])
end

plot(t, r, xlabel="Time", ylabel="Amplitude", legend=false, ylim=(-2.5, 2.5))

位相を調整する

測定信号や通信用の信号を生成したいときには、一般にクレストファクタは低ければ低いほどよいです。許された振幅を十分に使えるということは、信号・雑音比(S/N)を低くすることができるし、たとえばFM変調で伝送するときを考えると周波数帯域幅をまんべんなく使えるということですからね。

できるだけ全体にわたってピーク値が均等になるように、各倍音の位相を倍音の強さに応じて調整します。

 \displaystyle
\phi_n = \phi_1 - 2\pi\sum_{l=1}^{n-1}\left(n-l\right)p_l

それをプログラムにすると以下のようになります。

# phase angles
ϕ = zeros(N)
ϕ[1] = -π/2
for n in 2:N
    tmp = 0
    for l in 1:(n-1)
        tmp += (n - l) * p[l]
    end
    ϕ[n] = ϕ[1] - 2π * tmp
end

# signal generation (Eq.1)
fs = 1000
T = 1
t = range(-π/2, stop=π/2, step=1/fs)
r = zeros(size(t))
for k in 1:N
    r .+= (p[k] / 2)^(1/2) * cos.(2π * k * t / T .+ ϕ[k])
end

plot(t, r, xlabel="Time", ylabel="Amplitude", legend=false, ylim=(-2.5, 2.5))

パワースペクトルを調べると、各倍音の強さは両方とも同じになっています。位相がずれているので異なる信号なのですが、定常音についてヒトは位相に鈍感なので、実際にスピーカから再生すると同じ音に聞こえます。波形はまったく違うのに、同じ音に聞こえるって不思議ですね。