tb = Take[tbl, {1, 3*sr}];
g = ListPlot[tb, PlotRange -> All]
sr から3秒分のデータを取り出してプロットしてみた。
(サンプリング周波数が sr = 44.1 kHz なので、1 から 3∗sr で、3 秒分のデータということにな
る。Take[] はリストから、指定した範囲のデータを取り出す関数である。)
tb = Take[tbl, {62800+1, 62800+sr}];
g = ListPlot[tb, Joined -> True, PlotRange -> {{1, 1600}, {-0.3, 0.3}}]
音が鳴り始めるのは62800番目辺りからなので、そこから1秒分 (sr×1 s = 44100 個のデータを)取 り出して、1600個分 (1600/44100≒0.036 秒分) プロットしてみた。ここは色々試してみると良い。
ListPlay[tb, SampleRate->sr]
とすると、取り出したデータ tb の音を鳴らすことが出来る。
c = Fourier[tb];
ListPlot[Abs[c], Joined->True, PlotRange->All]
tbの離散Fourier変換cを求め、絶対値をプロットした。これから周波数の分布が読み取れる(はず)。
Abs[]の代わりにRe[],Im[]としてみたり(分かりにくいけれど、Cn =CN−nの関係が見える？)。
(* n1〜n2 の範囲で |c[[n]]| をプロットする。 *)
graph[c_, n1_, n2_] := ListPlot[Abs[c], Joined -> True, PlotRange -> {{n1, n2}, {0, Max[Abs[c]]}}]
graph[c, 1, 1600]
graph[c, 120, 140]
範囲を区切って表示することで、ピークを探してみた(左右対称なので左側だけで探す)。素朴に 目で見て探したが、プログラムを書いて自動化することも難しくないであろう。
ピークは 130 番目である。つまり |C129| が最大ということである。(リスト c の1番目の要素 c[[1]] は C0 であり、リストの要素の番号とFourier 係数のインデックスが1 ずれていることに注 意する。)
これはこのギターの音の基本周波数が 129 Hz (ドの周波数 131 Hz に近い — ぴったりでないの は、録音をした人 (研究室の学生)のチューニングが完璧ではなかったから) であることを意味する。
手近の楽器の音 (ピアノ，リコーダー，…) を自分で録音出来る人は、試してみよう。
逆Fourier変換して元に戻るか？
tb2=InverseFourier[c];
Norm[tb-tb2]
実数データにするため、tb2=Re[InverseFourier[c]]; とする方が良いかもしれない。
余裕があれば、http://nalab.mind.meiji.ac.jp/~mk/fourier/mathematica-memo/node9.html も見てみよう(高い周波数の音をカットする実験)。
4.2 PCM による音のデジタル信号表現
音とは空気中を伝播する縦波である。音があるとき、気圧が音がないときの圧力を基準に時間変化 する。基準圧力からの変位を音圧と呼ぶ。音圧の時間変化を記録することで音を記録(録音)出来る。
PCM (pulse code modulation, パルス符号変調)とは、アナログ信号(連続変数の関数)をデジタ ル信号(離散変数の関数— 数列)で表現するための1つの方法であり、音楽用CD,コンピューター のデジタル・オーディオ,デジタル電話等で標準的な形式となっている。具体的には、次の二つに基 づく。
(a) 一定の時間間隔で信号の値を評価する (サンプリングする、という) (b) 信号の値は有限桁の数で表現する (量子化する、という)
特に値の属する区間を等間隔に区切って、もっとも近い値に丸めることで実現するとき、LPCM (linear PCM) という。
パソコンで音を記録 (録音) する場合は、電気的な信号に変換された音声を数値化して(アナログ 信号からデジタル信号に変換するのでAD 変換(analog-to-digital conversion)という)、数値を記録 することになる。
サンプリング周波数とは、1秒間に何回サンプリングするかを意味している。サンプリング周波 数が高いほど、より高い周波数の音が記録できるようになる。
音楽用CD (1980年にSONYとPhillipsにより規格化された4)では、サンプリング周波数44.1 kHz が採用された。
サンプリング周波数が 44.1kHz とは、1秒間に 44100 回のデータを記録するということになる。
これが採用された理由は主に次の理由による。
(a) 人間が普通聞くことが出来る音の周波数は 20 Hz∼20 kHz と言われている。
(b) サンプリング周波数は、再生したい最も高い音の周波数の2倍以上にする必要がある(これは後で 解説するサンプリング定理(定理5.0.1)を根拠とする、と説明されるが、定理3.2.4(2016/11/16 現在まだ説明していない) で説明することも出来るだろう。)。
つまり、人間が普通に聞くことが出来る音を記録・再生出来るためには、サンプリング周波数は 2×20 kHz = 40 kHz より高くする必要がある。
コンピューターで処理することを考えると、「有限桁の数」は、「2進法の有限桁の数」というこ とになるが、その際の桁数 (ビット数)を量子化ビット数と呼ぶ。
8ビットの場合は 28 = 256段階、16ビットの場合は 216= 65536段階で表現することになる。音 楽用CD では、量子化ビット数として、16ビットが採用された。
サンプリング周波数と量子化ビット数の値が大きいほど、元の信号により忠実なデータが得られ るが、もちろんデータの量はそれだけ増大する。
余談 4.2.1 (CDプレーヤー発売当時のパソコン技術の相場) 音楽用CDでは、ステレオ (2 ch) が 普通なので、1秒間あたり、
44.1 k×16 b×2 = 1411.2 kb = 176.4 kB = 172.266 KB
のデータが流れることになる。1分間では、その60倍の10.0937 MB (ここで 1 MB = 1024 KB と いう意味) が流れることになる。1982年当時に普及していたリムーバブルな外部記憶媒体であるフ ロッピー・ディスクの容量は 1.2 MB 程度だったので、1分の音声信号を記録するのに、10枚近い
4音楽用のCD プレーヤーが初めて販売されたのは1982年である。当時ようやく16ビットCPU を用いたパソコ ンが市販された頃であった。外部記憶装置としては、1.2 MBの容量のフロッピー・ディスクが広く使われていた。
フロッピー・ディスクが必要だったことになる。これでは全然実用的ではない。CD という新しい メディア (74 分程度記録出来るようにするため、740 MB の容量となった) が必要になったのは当 然のことである。
それにしても、2016年現在、フロッピーディスクはもちろん、音楽用 CDさえもあまり目にしな くなってしまった。この種の技術の切り替わりの速さは著しい。
その後普及した MP3では、圧縮によって、データのサイズを 1
10 程度まで小さくすることが出 来るようになった。
4.3 結果の分析
4.3.1 一般論の復習
(これまで独立変数は x と書くことが多かったが) 音声信号を考えるとき、独立変数は時刻なの で、ここではt と書くことにする。
信号の値そのものは x で書いてある、つまり信号を x(t) としている本が多いので、ここでもそ れに従う。
周期 T の周期関数 x: R→C は次のようにFourier 級数展開出来る。
x(t) =
∑∞ n=−∞
cnei2πntT (t∈R), (4.1)
cn= 1 T
∫ T
0
x(t)e−i2πtT dt (n ∈Z).
(4.2)
基音の周波数は、周期の逆数 f = 1
T である。第 n 項の周期は T
|n|, 周波数は |n|f. n0 倍音の周 波数n0f に対応するのは、n =±n0 の項である。
一周期区間[0, T]に N 回測定 (サンプリング)すると、サンプリング周期Ts=T /N,サンプリン グ周波数 fs = NT でサンプリングすることになる。
区間 [0, T]の N 等分点 tj :=jTs でのx の値xj =x(tj)を用いる。このとき離散Fourier係数 Cn は次式で与えられる。
(4.3) Cn = 1
N
N∑−1 j=0
xjω−nj (n ∈Z), ω=e2πi/N.
離散フーリエ係数{Cn}は周期数列であるので、連続するN 項 {Cn}Nn=0−1 だけ求めれば十分である。
離散フーリエ係数 {Cn}Nn=0−1 は
Cn=∑
p≡n
cp を満たす。
{Cn} から{xj} を求めるには、逆離散Fourier 変換すれば良い:
(4.4) xj =
N∑−1 n=0
Cnωjn (j = 0,1, . . . , N −1).
4.3.2 今回の実習では
サンプリング周波数 fs = 44.1 kHz でサンプリングしたデータから、T = 1 s 分の信号 (N =
fsT = 44100個の数値) を取り出して、離散フーリエ変換した。
周期 T = 1 s の周期信号とみなして Fourier 級数展開したことになる。
4.3.3 | C
n| (1 ≤ n ≤ N − 1) は左右対称
u が実数値関数なので、cn =c−n が成り立つ。特に |c−n|=|cn|. 実際 cn = 1
T
∫ T
0
x(t)e−in2πTt dt= 1 T
∫ T
0
x(t)e−in2πTtdt = 1 T
∫ T
0
x(t)e−i(−n)2πTtdt =c−n. (Cf. 実数値関数 f の Fourier 変換fbに対して、fb(ξ) =fb(−ξ) が成り立つ。)
同様に離散フーリエ係数については、
Cn=C−n =CN−n, |Cn|=|C−n|=|CN−n| が成り立つ。
横軸 n (1≤n ≤N),縦軸 |Cn| でプロットすると、左右対称になるのはこのためである。
(細かいことをいうと、実数値データを Fourier[] で離散フーリエ変換するのは、無駄な作業を していることは否めない。それを避けるための手段 (離散Fourier余弦変換など)はあるが、この講 義では大事な原理が分かればそれで良い、という方針なので、そういう話題は省略する。)
4.3.4 第 n 項の周波数は | n | /T
(4.1) の第n項の周期は T
|n|,言い換えると周波数は|n|
T =|n| Hz である。
c1, c−1 は 1 Hz の成分 c2, c−2 は 2 Hz の成分
... ...
|C129|=|CN−129|が最大ということは、周波数が129 Hzの成分が最大ということを意味している。
その次に大きいのは |C258| = |CN−258| であった。258 Hz の成分がその次に大きいことを意味し ている。一番低い 129 Hz の整数倍になっているのは、ギターが1次元の振動現象であることから、
もっともである。
参考 4.3.1 (1次元の弦の振動) 長さ L の弦の(微小な) 振動は、次の波動方程式の初期値境界値問 題をモデルに持つ。
1
c2utt(x, t) = uxx(x, t) (0< x < L, t >0) u(0, t) =u(L, t) = 0 (t >0)
u(x,0) =ϕ(x), ut(x,0) = ψ(x) (0≤x≤L).
u=u(x, t)は、釣り合いの位置xにあった点の時刻tにおける変位(基準の位置からのずれ)を表す。
T を張力、ρ を線密度(単位長さあたりの質量) として、cはc=√
T /ρ で与えられるが、これは 実は弦を伝わる波の速さに等しい。ヴァイオリンやギターなどでは、張力を変えることで音の高さ を調整することが出来る。
この問題の解は次式で与えられる。
(♯) u(x, t) =
∑∞ n=1
sinnπx L
(
ancoscnπt
L +bnsincnπt L
) ,
an= 2 L
∫ L
0
ϕ(x) sinnπx
L dx, bn= 2 cnπ
∫ L
0
ψ(x) sin nπx L dx.
(♯)の第n項の基本周期は 2L
nc,その周波数は nc
2L. これらは(n = 1に対応する)一番低い周波数 c の整数倍である。 2L
余談 4.3.1 (音階) 1オクターブ高い音の周波数は2倍である。西洋の音階では、1オクターブは、
半音12個分に相当する(C, C#, D, D#, E, F, F#, G, G#, A, A#, B)。平均律では各音の周波数が 等比数列になる。よって、半音高いと周波数は 21/12= 1.05946· · · 倍になる。ピアノの鍵盤の中央 付近にあるアーA (ラ) の音は、普通440 Hz に調律されるので、その下のツェーC (ド)の音 (半音 9個分低 い) の周波数は
440
29/12 = 261.6255653· · · Hz.
上の実験のギターの音は、これより1オクターブ低い、261.6
2 = 130.81· · · を目安に調律したので あろう。実際に|C129|が最大になったのは、まあまあのチューニングだったのであろうか？(私はこ のあたり詳しくないので…知っている人は教えて下さい。)
4.3.5 より精密に
(話が細かいので、実習日にはカットする。)
実際には、周波数は自然数に限られるわけではない。その場合は T = 1 s は、その信号の周期に ならない可能性がある。
簡単のため、周波数 f の信号 u(t) =e2πif t を考える。0≤t ≤T で記録して、(周期T の周期関 数としての) Fourier 係数を求めると、
cn= 1 T
∫ T 0
u(t)e−in2πTtdt = 1 T
∫ T 0
e2πi(f−n/T)tdt= 1 T
∫ T 0
eiAntdt= 1 iAnT
(eiAnT −1) . ただし An:= 2π(f−n/T). これから
|cn|= sincAnT 2 .
T = 1 s, f = 130.813 Hz のとき、n= 125,· · · ,135 の範囲で sinc(AnT /2)を調べると、
(おっと、書きかけだ。)
4.4 Mathematica での音の取り扱い
実は Mathematicaの Fourier[ ] がデフォールトで計算するのは Cn′ = 1
√N
N∑−1 j=0
xjωnj
である。Cn を計算させるには、FourierParameters->{-1,-1}というオプションを与えれば良い。
c=Fourier[tb, FourierParameters->{-1,1}];
ところで、信号が実数値であることから、
Cn = 1
√N Cn′
という関係が成り立つ5。ゆえにパワースペクトルについては、|Cn|2 = N1 |Cn′|2 であるから、周波 数を求めたりする場合は(パワースペクトルが大きくなる n はどこか調べる)、デフォールトのまま 使っても良い。
最近のコンピューターを用いて音データを取り扱うための情報については、(かなり雑ではあるが)
「音の取り扱いに関するメモ」6 が参考になると思われる。Import[], ListPlay[], Fourier[] など の使い方を説明してある。
5xjωnj=xjωnj =xjω−nj であることに注意せよ。
6http://nalab.mind.meiji.ac.jp/~mk/labo/text/memo-sound/