Waveファイルからウダー譜を作成するPythonスクリプトを公開します。
ウダーとは？
宇田道信さんが開発した螺旋上を指で押さえて演奏する新しい楽器です。1周で1オクターブで、半音より小さい間隔の音も出すことが出来ます。
ウダー譜とは？
宇田さんが考案した、新しい音楽の表記方法です。ウダー同様に1周で1オクターブになっています。軽く探してみましたがまだネット上にスクリーンショットなどはないようです。
ウダー譜作成ソフトUdaFFTとは？
宇田さん本人によるMIDIをウダー譜で表示するプログラムを拝見した筆者が勝手に作ったWAVをウダー譜で表示するプログラムです。
コンセプト
従来の五線譜は直感的ではありません。同じ位置にある音符が、行頭で宣言されているシャープやフラットといった記号の違いで違う音になります。また、同じ和音であっても根音の違いによって五線譜の上での表記は異なったものになります。ウダー譜では同じ和音は同じ形で表されます。また、まだ普及版ウダーが販売されていないので実践は出来ませんが「ウダー譜に表示されているとおりに押さえれば出したい音が出る」という点が五線譜にはないメリットになるでしょう。既存の音楽ファイルをウダー譜に変換するプログラムがあれば、見たとおりに弾くだけで聞いたとおりの音が出るはずです。
まぁ、理由付けはさておき、筆者はウダー譜に美を感じており、ウダー譜でカノンを表示すると美しいのではないか、ウダーでカノンを弾くと楽しいのではないか、と考えています。
UdaFFT Ver0.1
現時点では残念ながら簡単に使っていただけるような形になっていません。まずSoundEngine Freeなどを使って8bitモノラルのWAVファイルを用意します。それをUdaFFT Ver0.1にかけると0.1秒あたり1枚のBMPファイルが出力されます。それをAVIMakerを使って無圧縮AVI形式の動画にし、携帯動画変換君を使ってH.264形式の動画に変換しています。最新のQuickTimeで見ることが出来るはずです。実行にはPython Japan User's Groupと、numarray、Python Imaging Library (PIL)のインストールが必要です。
原理
まずWavファイルの波形のデータから周波数ごとの強さのデータを作るためにフーリエ解析を用います。フーリエ解析は自分で実装してもいいですが、すでにいろいろな言語で高速フーリエ解析(FFT)の実装があるのでそれを参考にするのもいいでしょう。Mathtools.net : Java/FFT。UdaFFT Ver0.1ではnumarrayというライブラリのFFTモジュールを利用しています。38-39行目で、音の波形データから0.1秒分を切り出し、それをフーリエ解析した結果の絶対値を計算して周波数ごとの音の強さを求めています。「0.1秒」という長さは適当に決めました。
13行目:
from numarray import FFT
38-39行目:
f = data[start : start + framePerStep]
s = [abs(v) for v in FFT.fft(f)]
元のWavファイルのサンプリングレートが44KHzであれば0.1秒分で約4400個のサンプルがFFTに渡されることになります。FFTの結果sは、s[1]が「0.1秒で1周するような周期の強さ」、s[2]が「0.1秒で2周するような周期の強さ」になっています。つまり周波数に直すとそれぞれ1秒で10周するので10Hz、1秒で20周するので20Hz…となります。UdaFFT Ver0.1では41行目で周波数を求めた後、42行目でその周波数が200未満の場合と3200以上の場合を思い切って捨てています。目的の音以外のノイズを表示しないように、という意図です。この下限上限は適当に決めました。
41-42行目:
freq = frameRate * j / framePerStep
if freq < 200 or freq > 3200: continue
さて、周波数ごとの音の強さが求まったところで、次はこれをオクターブごとにまとめましょう。1オクターブ上がると周波数は2倍になります。2オクターブなら4倍です。周波数の比から、それが何オクターブ分に当たるかを計算するには対数を使います。UdaFFT Ver0.1では43行目で対数を使って周波数freqが真ん中のラから見て何オクターブ離れたところにあるかを計算しています。式の最後の「+5」は特に重要な意味はなく、負の数を切り捨てた場合の挙動を調べるのが面倒だったので全部正の数になるように足したものです。
43行目:
l = log(freq / FREQ_CENTER_A) / log(2) + 5
UdaFFT Ver0.1はこの後、この値を1オクターブ64個の音階に丸めていますが、これは改善の余地があると思います。
最後に音階から画面上の位置への変換について説明します。先ほど対数で求めた値lの整数部分はオクターブの情報なので、「l - int(l)」とやって小数部分を取り出します。この値が0～1で0度～360度に対応するので、360を掛ける…と言いたいところですが、度ではなくてラジアンを使う言語の方が多いのでここは2π(6.283)を掛けます。この「ラジアン表記された角度」から画面上の位置x, yを求めるには、三角関数を使います。今回は0ラジアンのラが時計で9時の位置に来て時計回りに値が増えるようにします。そうするとラどドの間が3半音なのでドが12時の位置に来ます。0ラジアンのラが9時の位置に来るということは、左上隅が原点のスクリーン座標系では円の中心から見てX軸方向にマイナスでY軸方向に0ということです。sin(0)は0でcos(0)は1なので、X座標の値は「中心のX座標 - 半径 × cos(ラジアン)」となります。次に、0.1ラジアンの場合を考えると、時計回りなのでXが少し増えてYが少し減ります。sin(x)は0のそばではxが増えるとsin(x)も増えます。つまりY座標の値もマイナスがついて「中心のY座標 - 半径 × sin(ラジアン)」となるわけです。
サンプル動画
ドレミ(219KB)、
肉声(1.17MB)
iTunesか最新のQuickTimeをインストールしている必要があるかと思います。
既知の問題点と今後の方針
- 人間の声は想像以上に倍音成分が多いようです。根音から見て7時の位置に二つめの音符が出てしまいます。
- 和音と倍音の区別は難しいように思います。和音の表示をあきらめれば「一番強い音のところにだけ音符を表示する」という解決策も考えられます。
- どのくらいの音の大きさを白(255)にするかを手動で設定しなければならない。
- 「直前のステップで一番強かった音の強さ × 0.8」を上限にすれば柔軟に対応できるかも？
- オクターブの情報を捨てている。
- 宇田さんのソフトでは丸の大きさでオクターブの違いを表現している。
- ウダー譜を得るまでの操作がかなり面倒。
- Javaで実装して宇田さんのソフトに組み込んでもらう？
- WindowsMediaPlayerの視覚エフェクトとして実装するのが楽かも知れない。MSDNによれば、フーリエ解析を視覚エフェクトがやらなくても周波数ごとの強さのデータを取得できるようです。
ソースコード
ウダー及びウダー譜の普及に少しでも役立てるように、以下のソースコードに関しては商用・非商用を問わず複製・改変・頒布を認めます。
# -*- coding: cp932 -*-
#
# UDA-FFT
#
from time import clock
startTime = clock()
import os
import wave
import Image
import ImageDraw
from math import sqrt, pi, log, sin, cos
from numarray import FFT
PROJECT_DIR = r"c:\Home\Projects\UdaFFT"
INPUT_WAVE = "hana-8bitM.wav" # 8-bit モノラル
OUTPUT_DIR = "bmps"
wav = wave.open(os.path.join(PROJECT_DIR, INPUT_WAVE))
(nChannels, sampWidth, frameRate, nFrames, compType, compName) = wav.getparams()
RADIUS = 256
SIZE = RADIUS * 2
print frameRate
framePerStep = frameRate / 10 # 0.1 sec
numStep = nFrames / framePerStep
NOTE_PER_OCTAVE = 64
RENDER_LENGTH = numStep # 通常numStep、デバッグ時は減らして軽くする
FREQ_CENTER_A = 440.0 # 真ん中のラ
data = wav.readframes(RENDER_LENGTH * framePerStep)
print clock() - startTime
print "RENDER_LENGTH:", RENDER_LENGTH
powerMap = {}
for i in range(RENDER_LENGTH):
start = framePerStep * i
f = data[start : start + framePerStep]
s = [abs(v) for v in FFT.fft(f)]
for j in range(1, len(s)):
freq = frameRate * j / framePerStep
if freq < 200 or freq > 3200: continue
l = log(freq / FREQ_CENTER_A) / log(2) + 5
angle = int(NOTE_PER_OCTAVE * (l - int(l))) / float(NOTE_PER_OCTAVE)
th = 2 * pi * angle
v = s[j]
if powerMap.has_key(i) and powerMap[i].has_key(th):
powerMap[i][th].append(v)
else:
if powerMap.has_key(i):
powerMap[i][th] = [v]
else:
powerMap[i] = {th: [v]}
if i % 100 == 0: print i, "analysed"
wav.close()
print clock() - startTime
r = 10
SCALE = 750
LBOUND = 100
INNER_RADIUS = 40
PIXEL_PER_STEP = 15
for start in range(RENDER_LENGTH - 10):
image = Image.new("L", (SIZE, SIZE))
d = ImageDraw.Draw(image)
d.ellipse([(RADIUS - INNER_RADIUS, RADIUS - INNER_RADIUS),
(RADIUS + INNER_RADIUS, RADIUS + INNER_RADIUS)],
LBOUND)
for i in range(10):
rad = INNER_RADIUS + r + i * PIXEL_PER_STEP
for th in powerMap.get(start + i, []):
x = int(RADIUS - rad * cos(th))
y = int(RADIUS - rad * sin(th))
v = powerMap[start + i][th]
light = max(v) / SCALE
if light > LBOUND:
d.ellipse([(x - r, y - r),
(x + r, y + r)],
light)
image.save(os.path.join(PROJECT_DIR, OUTPUT_DIR, "%04d.bmp" % start), "BMP")
if start % 100 == 0: print start, "outputed"
print "OK."
print clock() - startTime
