Common LispでFFTライブラリを使ってみる

Common Lispを使っていろいろと車輪の再実装を続けています。今回はフーリエ変換ライブラリを試してみます。Common Lispにはフーリエ変換のライブラリは入っていないので、Quicklispでfftをキーワードに探していきます。すると以下の4つがヒットしました。

CL-USER> (ql:system-apropos "fft")
#<SYSTEM bordeaux-fft / bordeaux-fft-20150608-http / quicklisp 2019-12-27>
#<SYSTEM fft / fft-20180711-git / quicklisp 2019-12-27>
#<SYSTEM napa-fft3 / napa-fft3-20151218-git / quicklisp 2019-12-27>
#<SYSTEM pfft / fft-20180711-git / quicklisp 2019-12-27>

これらに関係しそうなのは以下のページのようです。fftとpfftは同じライブラリに含まれるもののようです。他にもf2clというプロジェクトにもFFTは含まれているようでした。

• https://github.com/ahefner/bordeaux-fft
• http://nklein.com/software/cl-fft/
• https://github.com/pkhuong/Napa-FFT3

どれが良いのかわからないので、なんとなく最初に出てきたBordeaux-FFTを使ってみました。（Napa-FFT3も高機能そうで良いかもです）

インストールと読み込みはQuicklispを使えば簡単でした。（Windows 10とSBCLを使っています）

(ql:quickload 'bordeaux-fft)

長さ16の正弦波（周波数は3 Hzとします）を作って、それにFFTをかけました。sfftというのはsimple fftの略あるいはstupid fftの略らしく、なにも考えずに配列のFFTを計算してくれるものです。fftは配列要素の型を指定したりする必要があるのですが、今回はシンプルなほうを使います。

;; 長さを16とする（2のべき乗にする）
(defparameter *N* 16)

;; 正弦波を入れる配列を作る
(setq sinusoid (make-array *N*))

;; 3 Hzの正弦波を要素ごとに計算
(loop for ii from 0 to (1- *N*)
      do (setf (aref sinusoid ii)
               (sin (* 3 2 pi (/ ii *N*)))))

;; FFTをかけてfftresにバインド
(setq fftres (bordeaux-fft:sfft sinusoid))

この結果、fftresには以下のような16要素の複素数配列が入りました。

#(#C(-2.4322862379963016d-16 0.0d0)
  #C(6.107235716236265d-16 1.4695673325073857d-15)
  #C(2.4963710534292284d-16 4.373545779112952d-16)
  #C(-4.5862415333198d-15 -8.0d0)
  #C(1.588639366831878d-15 -7.216449660063518d-16)
  #C(7.669517010585668d-17 0.0d0)
  #C(4.851509741454889d-16 -6.728684467138615d-16)
  #C(5.6559388869686625d-16 1.8962923073275474d-15)
  #C(-1.4644739508873025d-15 0.0d0)
  #C(-1.3244021561850481d-17 -2.357745752207511d-15)
  #C(4.851509741454891d-16 6.728684467138614d-16)
  #C(7.428289848809507d-16 -1.3322676295501878d-15)
  #C(1.588639366831878d-15 7.216449660063518d-16)
  #C(-1.255572459444331d-15 8.0d0)
  #C(2.496371053429228d-16 -4.373545779112951d-16)
  #C(9.20064081065035d-16 -1.1993546792729692d-16))

最初の#C(-2.43d-16 0.0d0)はDCオフセットの値で、2番目の要素以降が各周波数ビンの値です。4番目の#C(-4.59d-15 -8.0d0)と14番目の#C(-1.26d-15 8.0d0)が3 Hzの振幅値に対応しています。絶対値が8になっているのは正の周波数と負の周波数に分割されているからで、合計すると16になります。これを信号の長さで正規化すると（各要素を信号長16で除算すると）1.0になります。

CL-USER> (abs #C(-4.5862415333198d-15 -8.0d0))
8.0d0
CL-USER> (abs #C(-1.255572459444331d-15 8.0d0))
8.0d0

配列の最初から、DC成分、正の周波数成分、ナイキスト成分、負の周波数成分という順に値が入っているので、正負の周波数成分をまとめます。このあたりは以前Julia版を書いたときに説明したので、そちらを参照してください。

marui.hatenablog.com

ここでは、まずampspecという振幅スペクトルを入れる配列を作ります。この配列の長さはN/2+1となります（今回は16/2+1で長さ9です）。次に、FFT結果（fftres）の配列の0番目（DC成分）とN/2-1番目（ナイキスト成分）をampspecにコピーします。そして、それらの間の成分について、負の周波数成分の複素共役を正の周波数成分に加算します。

(setq ampspec (make-array (1+ (/ (array-dimension fftres 0) 2))))

(setf (aref ampspec 0) (aref fftres 0))

(setf (aref ampspec (1- (array-dimension ampspec 0)))
      (aref fftres (1- (array-dimension ampspec 0))))

(loop for ii from 1 to (- (array-dimension ampspec 0) 2)
      do (setf (aref ampspec ii)
               (+ (aref fftres ii)
                  (conjugate (aref fftres (- *N* ii))))))

結果、ampspecに得られるのが複素数の周波数スペクトルです。

#(#C(-2.4322862379963016d-16 0.0d0)
  #C(1.5307876526886614d-15 1.5895028004346827d-15)
  #C(4.992742106858457d-16 8.747091558225903d-16)
  #C(-5.841813992764131d-15 -16.0d0)
  #C(3.177278733663756d-15 -1.4432899320127035d-15)
  #C(8.195241549868074d-16 1.3322676295501878d-15)
  #C(9.70301948290978d-16 -1.3457368934277227d-15)
  #C(5.523498671350158d-16 4.254038059535058d-15)
  #C(-1.4644739508873025d-15 0.0d0))

振幅スペクトルにするには、配列の各要素にabsを適用して計算するのですが、配列にはmapcarが使えないのでリストに変換してからmapcarをかけてあげます。（配列にmapをかけるarray-mapという関数がStack Overflowにあったので、それを使ってもいいかもです）

CL-USER> (mapcar #'abs (coerce ampspec 'list))
(2.4322862379963016d-16 2.2067691293413003d-15 1.0071697199260125d-15 16.0d0
 3.489725774217968d-15 1.5641473323680597d-15 1.6590640907420562d-15
 4.289747100668858d-15 1.4644739508873025d-15)

さらに信号長で正規化すれば、所望の振幅スペクトルが得られます。ちゃんと3 Hzに対応した場所に振幅1.0があり、それ以外の要素はゼロ近傍になっていますね。

CL-USER> (mapcar #'(lambda (x) (/ x *N*))
                 (mapcar #'abs (coerce ampspec 'list)))
(1.5201788987476885d-17 1.3792307058383127d-16 6.294810749537578d-17 1.0d0
 2.18107860888623d-16 9.775920827300373d-17 1.0369150567137852d-16
 2.6810919379180364d-16 9.15296219304564d-17)

次は、これを関数化するのと、WAVファイルから読み込んだ音に適用するのをやっていきたいです。

本日の給油（CRF1100Lアフリカツイン／2BL-SD10）

youtu.be

給油日	オドメーター (km)	給油量 (L)	単価 (円/L)	燃費 (km/L)	距離単価 (円/km)
2021-10-24	10	9.69	155	-	-
2021-10-24	48	7.62	155	-	-
2021-10-30	267	13.46	162	16.27	9.96
2021-11-06	507	15.29	157	15.70	10.00
2021-11-14	737	13.01	157	17.68	8.88
2021-11-27	924	13.07	157	14.31	10.97
2021-12-25	1188	15.90	156	16.60	9.39
2022-01-10	1398	13.71	153	15.32	9.99
2022-02-05	1740	18.22	158	18.77	8.42
2022-03-07	2007	16.46	162	16.22	9.99
2022-04-10	2283	14.23	163	19.40	8.40
2022-06-26	2628	18.97	166	18.19	9.13
2022-07-18	2860	12.74	160	18.21	8.78
2022-07-31	3128	13.63	157	19.66	7.99
2022-08-27	3401	16.63	157	16.42	9.56

首都圏のなかで少し遠出をしたりはしたのですが、渋滞がちな道を通ることが多かったので、燃費は悪くなってしまいました。それでも24リットルのタンクなので、無給油でゆうに400 kmは走れます。遠くに行くときに給油計画をきっちりしなくていいのが良いですね。

1/3オクターブ周波数のリストを得る

目次

標準数とは

1/3オクターヴ幅での周波数リストが欲しい時があります。20, 31.5, 40, 50, 63, 80, 100, 125, 160, 200, 250, 315 Hz... というやつです。実はこのリスト、つまり1/3オクターヴごとの周波数の値の決め方についてはJIS Z8601「標準数」やISO 3 "Preferred Numbers"に規定があります。

ja.wikipedia.org

その前に、オクターヴと周波数の関係についておさらいしましょう。物理的な1オクターヴは周波数2倍と同等ですが、1/3オクターヴの周波数が欲しいとき、たとえば基準周波数を1000 Hzとしたら、それに2^1/3、2^2/3、2^3/3を掛け算すると、それぞれ1259.9、1587.4、2000 Hzとなります。1000 Hzの1/3オクターヴ上が1259.9 Hz、さらにその1/3オクターヴ上が1587.4 Hzという具合です。さて、オクターヴと周波数は指数の関係にありますので、1000、1259.9、1587.4、2000という値は対数軸上で等間隔になっています。

オクターヴに似たものにディケードがあります。1ディケードは周波数10倍で、上記と同様にして10^1/10、10^2/10、10^3/10を掛け算すると、それぞれ1258.9、1584.9、1995.3 Hzになります。

オクターヴとディケードのどちらを採用するかは、分野によるのではないかと思いますが、ディケードを採用すると前述のように1000 Hzを基準に周波数がぴったり2倍の2000 Hzは出てきませんし、オクターヴを採用すると1000 × 2^10/3 = 10079.4と、周波数10倍が10000ちょうどになりません。

計算上でぴったりとした周波数が必要であれば上記のように計算すればよいのですが、見た目がきれいな値でも良い場合があります。そのときに使えるのが標準数です。扱いやすい数列を用意しておいて、その2倍や10倍などを計算することで最初に挙げたようなキリの良い（良さそうに見える）数列を得るものです。JIS規格で「扱いやすい数列」にはR5, R10, R20, R40, R80などが用意されています。Rのあとの数は1ディケード（10倍）を対数軸上で何分割したかをあらわしていて、たとえばR5なら1.00, 1.60, 2.50, 3.00, 6.30の5つの数が使われます。3つ隣の値がオクターヴ関係になるという性質がありがたいからか、聴覚の臨界帯域と1/3オクターヴの相性がよいからか、音の世界ではR10（1.00, 1.25, 1.60, 2.00, 2.50, 3.15, 4.00, 5.00, 6.30, 8.00）が使われることが多いです。

Common Lispでの実装

この標準数を任意の周波数範囲で欲しいときに使えるプログラムをCommon Lispで作ってみました。let*のなかで計算に必要な数たちを計算しています。R10は前述の通りで、N1とN2は欲しい数がおよそ10の何乗のオーダーなのかを計算します。たとえば20～20000 Hzの間の数が欲しい時には、その範囲は10¹～10⁵に含まれますので、N1は1、N2は5にします。ffは何をしているかというと、R10の内容を（N1とN2をもとに計算した）10¹, 10², 10³, 10⁴, 10⁵それぞれと掛け合わせ、それをリストにしています。全ての組み合わせを掛け算するために2重ループを使い、1つのリストにするためにreduce #'appendを使っています。最後のループでは、指定した範囲の値だけをcollectで集めています。

;;; preferred frequencies
;;;
;;; (pref-freqs)
;;; => (20.0 25.0 31.5 40.0 50.0 63.0 80.0 100.0 125.0 160.0 200.0
;;;     250.0 315.0 400.0 500.0 630.0 800.0 1000.0 1250.0 1600.0
;;;     2000.0 2500.0 3150.0 4000.0 5000.0 6300.0 8000.0 10000.0
;;;     12500.0 16000.0 20000.0)
;;;
;;; (pref-freqs :min 10 :max 100)
;;; => (10.0 12.5 16.0 20.0 25.0 31.5 40.0 50.0 63.0 80.0 100.0)
;;; 
(defun pref-freqs (&key (min 20) (max 20000))
  "Calculate R10 series of preferred frequencies (JIS Z8601)."
  (let* ((R10 '(1.00 1.25 1.60 2.00 2.50 3.15 4.00 5.00 6.30 8.00))
         (N1 (floor   (log min 10)))
         (N2 (ceiling (log max 10)))
         (ff
           (reduce #'append
                   (loop for sc from N1 to N2
                         collect
                         (loop for fr in R10
                               collect (* (expt 10 sc) fr))))))
    (loop for f in ff
          if (and (<= min f) (<= f max))
            collect f)))

実行例

コメントにも使い方が書いてありますが、SBCLでの実行結果は以下のようになりました。

CL-USER>  (pref-freqs)   ; デフォルトでは20～20000 Hzの範囲
(20.0 25.0 31.5 40.0 50.0 63.0 80.0 100.0 125.0 160.0 200.0 250.0 315.0 400.0
 500.0 630.0 800.0 1000.0 1250.0 1600.0 2000.0 2500.0 3150.0 4000.0 5000.0
 6300.0 8000.0 10000.0 12500.0 16000.0 20000.0)

CL-USER> (pref-freqs :min 10 :max 100)   ; minとmax片方だけの指定も可能です
(10.0 12.5 16.0 20.0 25.0 31.5 40.0 50.0 63.0 80.0 100.0)

似たようなことをJuliaでやってみると……

だいぶ前に似たようなことをJuliaでも作っていました。上記の2重ループを使った方法とはアプローチが違うので紹介します。（Julia対Lispという比較をしたいわけではなく、あくまでも「アプローチの違いを味わいたい」という目的です。）

最低限の標準数の計算をするだけのプログラムは次のようなものになります。1行目で10のべき乗の列ベクトルと、R10の行ベクトルの積をとって、すべての組み合わせの積が入った行列を作ります。それを（値の並べ替えのために）転置してから[:]で1列にするというものです。使い勝手を良くするための周波数範囲指定などは入っていませんが、言語仕様に行列計算が組み込まれているのでここまで簡単にできます。おそらくCommon Lispでも行列計算ライブラリを使えば同じことができるでしょう。

julia> pref_freqs = 10.0 .^ range(1, 4) * [1.00 1.25 1.60 2.00 2.50 3.15 4.00 5.00 6.30 8.00]
4×10 Matrix{Float64}:
    10.0     12.5     16.0     20.0     25.0     31.5     40.0     50.0     63.0     80.0
   100.0    125.0    160.0    200.0    250.0    315.0    400.0    500.0    630.0    800.0
  1000.0   1250.0   1600.0   2000.0   2500.0   3150.0   4000.0   5000.0   6300.0   8000.0
 10000.0  12500.0  16000.0  20000.0  25000.0  31500.0  40000.0  50000.0  63000.0  80000.0

julia> pref_freqs = pref_freqs'[:]
40-element Vector{Float64}:
    10.0
    12.5
    16.0
    20.0
    25.0
    31.5
    40.0
    50.0
    63.0
     ⋮
 16000.0
 20000.0
 25000.0
 31500.0
 40000.0
 50000.0
 63000.0
 80000.0

周波数から音名を得る

Common Lispの勉強をしていて、簡単なプログラムを再実装しています。数年前にJuliaで書いた「周波数を与えると音名・セント値が返ってくる関数」をLispに移植してみました。Julia版は以下の記事に書いてあります。

marui.hatenablog.com

やっていることはJulia版と同じですし、あまり難しいことはしていないのでコードだけ張り付けておきます。（アルゴリズムの説明はJulia版の記事を読んでください）

(defun freq-to-midicent (frequency &optional (concert-pitch 440))
  "Convert frequency (Hz) to MIDI cent value (MIDI note number times 100
   for representing partial pitch). `(floor (/ (freq-to-midicent 262) 100))`
   can be used to get MIDI note number and cent value separately."
  (+ 6900 (* 1200 (log (/ frequency concert-pitch) 2))))

(defun freq-to-notename (freq &optional (concert-pitch 440) (flat nil))
  "Convert frequency (Hz) to note name with standard octave number (center C = 4).
   The boundary between octaves lies between the notes B and C."
  (let* ((notenames-flat '("C" "Db" "D" "Eb" "E" "F" "Gb" "G" "Ab" "A" "Bb" "B"))
         (notenames-sharp '("C" "C#" "D" "D#" "E" "F" "F#" "G" "G#" "A" "A#" "B"))
         (mc (freq-to-midicent freq concert-pitch))
         (notenum (round (/ mc 100)))
         (notecent (round (* (- (/ mc 100) notenum) 100)))
         (noteoct (floor (1- (/ notenum 12)))))
    (format nil "~A~A~@D"
            (nth (mod notenum 12) (if flat notenames-flat notenames-sharp))
            noteoct
            notecent)))

下のように使います。ピアノ鍵盤の中央のC音のオクターブ番号を4としています（アメリカ音響学会などのピッチ表記に合わせています）。また、デフォルトのコンサートピッチは440 Hzで、シャープ記号を使った音名表示にしています。

CL-USER> (freq-to-notename 1000)
"B5+21"
;; 1000 Hzの音名はB、オクターブ番号は5、ちょうどピッタリから21セント上

CL-USER> (freq-to-notename 270)
"C#4-45"
;; 270 Hzの音名はC#、オクターブ番号は4、ちょうどピッタリから45セント下

CL-USER> (freq-to-notename 270 442)
"C4+47"
;; コンサートピッチを442 Hzにしたとき

CL-USER> (freq-to-notename 270 440 t)
"Db4-45"
;; コンサートピッチを440 Hzにして、シャープではなくフラット表示にしたとき

本日の給油（CRF1100Lアフリカツイン／2BL-SD10）

給油日	オドメーター (km)	給油量 (L)	単価 (円/L)	燃費 (km/L)	距離単価 (円/km)
2021-10-24	10	9.69	155	-	-
2021-10-24	48	7.62	155	-	-
2021-10-30	267	13.46	162	16.27	9.96
2021-11-06	507	15.29	157	15.70	10.00
2021-11-14	737	13.01	157	17.68	8.88
2021-11-27	924	13.07	157	14.31	10.97
2021-12-25	1188	15.90	156	16.60	9.39
2022-01-10	1398	13.71	153	15.32	9.99
2022-02-05	1740	18.22	158	18.77	8.42
2022-03-07	2007	16.46	162	16.22	9.99
2022-04-10	2283	14.23	163	19.40	8.40
2022-06-26	2628	18.97	166	18.19	9.13
2022-07-18	2860	12.74	160	18.21	8.78
2022-07-31	3128	13.63	157	19.66	7.99

渡良瀬遊水地は群馬・栃木・埼玉・茨城の4県が出会っている場所にある大きな遊水池で、100年以上前に起きた公害の影響を低減させるために作られたものです（足尾鉱毒事件とか田中正造とか、小学校の社会の時間に習った記憶があります）。今は公園として整備され、魚釣りやデイキャンプなど自然の中の遊びができるようになっています。暑い中でしたが、今日は渡良瀬遊水地まで行って、ラーメン作って食べました。（火が使えるのは子供広場ゾーン内の一部で、焚き火台、BBQグリル、カセットコンロなどが使えます。直火は禁止。）

前回・今回と長距離のんびり走行が多かったので、これまでの中で最高の燃費になりました。信号待ちの時間が少ない道を見つけたからかもしれません。とはいえ、くねくねした道をSモードで走行したりもしてたので、抑えて走れば燃費はもう少し上げられそう。「目指せ、20 km/L！」です。

本日の給油（スーパーカブ110プロ／JA07）

給油日	オドメーター (km)	給油量 (L)	単価 (円/L)	燃費 (km/L)	距離単価 (円/km)
2021-07-31	12417.3	2.62	148.09	48.40	3.06
2021-09-25	12563.7	2.81	148.04	52.10	2.84
2021-11-04	12711.0	2.44	156.97	60.37	2.60
2022-01-03	12814.3	2.18	153.21	47.39	3.23
2022-02-23	12910.1	2.55	158.04	37.57	4.21
2022-03-06	13058.0	2.83	160.07	52.26	3.06
2022-05-26	13182.1	2.46	158.10	50.45	3.09
2022-06-22	13278.4	1.99	165.83	48.39	3.43
2022-07-27	13405.0	2.21	157.92	57.29	2.76

少しガソリン単価が下がったのかな。それでも距離単価2円を下回るまではだいぶありそう。ガソリン単価120円/Lを下回り、燃費が70 km/Lくらいになると距離単価が2円/kmを下回ってくる感じがします。

この8月で10年目を迎えるので、例年の定期点検に加えて前後輪ともにタイヤとチューブを新品に交換。本当はBridgestone TW8のようなブロックパターンのものを履かせたかったけど、メーカー在庫なしとのことで、IRCの純正品に。110プロはホイールサイズがやや特殊なのでしょうがないですね。他にも後輪ブレーキシューやらスパークプラグやらを交換してもらい、これからまた10年間たのしく乗れるようにしてもらいました。

sfplay~のディスクバッファサイズ設定

Max 8の[mc.sfplay~]を使って、29チャンネルのWAVファイルを3つ同時再生するというパッチを作っていたのですが、普通に組んだら、どうしても最初の1秒くらいしか再生されません。詳しい人に質問したら、8チャンネルを超えるファイルを再生するときには[mc.sfplay~]にバッファサイズの指定をしないといけないとのことでした。

バッファの大きさは様々なバイト長に対応できるように、「8の階乗÷2（つまり20160）」の整数倍にするのが基本みたいです。今回は96000 Hz / 24 bit / 29 chのファイルを再生したいので、[mc.sfplay~ 29 584640]と指定ところ問題なくファイルを読み込んで再生することができました（20160×29=584640です）。