はみたものの。
オハナシとしては Python の wave モジュール例の殴り書き ～ Python の wave モジュール例の殴り書き(6) からのそのままの流れだけれど、さすがに「wave モジュール」全然関係なくなったので。
ま、やってみたわけさ:
1 import numpy as np
2 import matplotlib.pyplot as plt
3 import matplotlib.ticker as ticker
4 from tiny_wave_wrapper import WaveReader, WaveWriter
5
6
7 _SCALES_L = ["C", "C#", "D", "D#", "E", "F", "F#", "G", "G#", "A", "A#", "B"]
8
9 def _nn2scale(d):
10 return int(d / 12) - 2, _SCALES_L[int(d) % 12]
11
12 def _nn2freq(d):
13 return np.pow(2, (d - 69) / 12.) * 440
14
15 def _freq2nn(f):
16 return 69 + 12 * np.log2(f / 440.)
17
18 def _nnformat(f, pos):
19 if f > 0:
20 nn = _freq2nn(f)
21 oc1, sc1 = _nn2scale(int(nn))
22 oc2, sc2 = _nn2scale(int(np.ceil(nn)))
23 return r"""%.2f
24 [%s(%d) $\sim$ %s(%d)]""" % (
25 nn, sc1, oc1, sc2, oc2)
26 return ""
27
28
29 if __name__ == '__main__':
30 import argparse
31 parser = argparse.ArgumentParser()
32 parser.add_argument("mode", choices=["show", "pdf"])
33 parser.add_argument("-s", "--step", type=int)
34 parser.add_argument("target")
35 args = parser.parse_args()
36
37 with WaveReader(args.target) as fi:
38 nchannels, width, rate, nframes, _, _ = fi.getparams()
39 raw = np.fromstring(fi.readframes(nframes), dtype=np.int16)
40 channels = raw[::2], raw[1::2]
41
42 if args.step:
43 step = args.step
44 else:
45 step = rate // 8
46 fig, ax = plt.subplots()
47 if args.mode == "pdf":
48 fig.set_size_inches(16.53 * 4, 11.69 * 2)
49 for chn in range(len(channels)):
50 ax1 = plt.subplot(2, 1, chn + 1)
51
52 F = np.array([])
53 channel = channels[chn]
54 nframes = len(channel)
55 nframes -= (nframes % step) # drop the fraction frames
56 for i in range(0, nframes, step):
57 f = np.abs(np.fft.fft(channel[i:i + step]))
58 f = f / f.max()
59 F = np.vstack((F, f)) if len(F) else f
60
61 freq = np.fft.fftfreq(F.shape[1], 1./rate)
62 X = np.arange(F.shape[0]) / (rate / step) * (2 / width)
63 Y = freq[:len(freq) // 2]
64 Z = F.T[:len(freq) // 2,:]
65 ax1.contour(X, Y, Z, cmap='jet')
66 ax1.set_ylabel("in Hz")
67 ax1.grid(True)
68
69 ax2 = ax1.twinx()
70 ax2.contour(X, Y, Z, cmap='jet')
71 ax2.yaxis.set_major_formatter(ticker.FuncFormatter(_nnformat))
72 ax2.set_ylabel("in scale (approx)")
73 ax2.grid(True)
74
75 if args.mode == "pdf":
76 plt.savefig(args.target + ".pdf", bbox_inches="tight")
77 else:
78 plt.subplots_adjust(right=0.75)
79 plt.show()
周波数と音階の関係とかは Python の wave モジュール例の殴り書き ～ Python の wave モジュール例の殴り書き(6) でやったことそのまま、なので、コードが長くなっちゃうんでそこいらに元は付けてた docstring は省いちゃった。
割とすっきりしてるんで苦労のあとを感じないかもしれないけれど、実際はチマチマとハマってた。右側のその音階を見せる部分の値が合わなくて、根本を揺るがすこと(件の「殴り書き」からしての間違い)が起こってるのかとビビったけど、「ax1 でも ax2 でも同じ引数で contour を呼び出す」と最初しなかったことが問題だった。オクターブが 5～6 で作った wav のつもりなのになぜか 8～9 あたりを指しちゃってて困った困った。ほかも説明するまでもないような程度のハマり(凡ミスみたいなやつ)はちょいちょい。でこういうの繰り返してるとさ、「核心部(つまり実際に fft してるとこ)が間違えてるんでは?」と不安になって困る。結果最初から最後まで「ここは」間違えてなかったんだけれど。
実際に適用すればこんな具合:
まぁ要するに「あんまし使いやすくない」のよね。「楽譜的」に見えたりせんかなぁと思ったけれど、そうするにはかなりもっと工夫しないといけないだろう。
ただこれ、「理解を深める」のにはとっても役にたった。それに「使いやすくない」と言ったって、Python の wave モジュール例の殴り書き で作ったような「作為的な音声」を「正しく作れてるか」であるとか、Python の wave モジュール例の殴り書き(6) で試みた bandpass フィルタの検証には嬉しいわけよな、こんなでも。
