第10回 離散フーリエ変換とスペクトル分解(6月16日)


*第二学生実習室にて行います。


フーリエ変換とは


フーリエ変換

逆フーリエ変換(フーリエ積分)

この2つの式によって、時間領域 と周波数領域 を自由に行き来することが可能になる。時間領域では、波の波形や周期などを知ることができ、周波数領域ではどのような周波数がどれだけ含まれているかを知ることができる。



標本化(サンプリング)と標本化(サンプリング)定理


 アナログ信号をデジタル信号に変換する仕組みをサンプリングという。音楽を記録する最初の機械は1877年のエジソンによる蓄音機とレコードである。レコードでは、連続的な波ををそのまま溝に刻み込んでいたので、蓄音機で再生するたびにレコードの溝は少しずつ削られて音が劣化してしまう弱点があった。

 

音の劣化問題を解決したのはCD (Compact Disk)であった。CDでは、一定の時間間隔で波形を読み取り、離散的なデータとして音を記録している。CDが劣化に強いのは、アナログ信号が実数値をとるのに対し、デジタル信号は整数値しかとらないためである。(図で実数値を取っているところは、デジタル信号では、四捨五入され整数値として処理される)

 

サンプリングは、ある時間間隔(サンプリング周期、標本化周期)でアナログ信号を読み取る標本化と、標本化したアナログ信号を数として記録する量子化の2つの手順でアナログ信号をデジタル信号へ変換する処理である。

 前回学んだように、サンプリング周波数(1秒間に何回データを取るか)は、とても重要である。では、一体どれくらいの数を取ればいいのか。その指標を表すものが「サンプリング(標本化)定理」である。

左図に示すように、正弦波をデジタル信号として適切に表現しようとするとき、1周期あたり少なくとも2個のデータが必要である。1個のデータでは、山と谷が表現できない。つまり、サンプリング周期 と波の周期 には、

 

という関係が成り立っていなくてはならない。ここで、周期ではなく周波数で考えると、周波数は周期の逆数なので、サンプリング周波数 と波の周波数 では、



という関係式となる。この式はつまり、周波数 の正弦波をサンプリングするには、サンプリング周波数 を の2倍以上に設定しなければならない、という意味である。

 さて、実際にはサンプリング周波数 は、波に含まれる最大の周波数 の2倍以上に設定する必要があるが、最大の周波数 は未知である場合が多い。そこで、音の解析では、人間の聴覚の最大周波数である20 kHzの2倍以上に設定する。例えば音楽CDではサンプリング周波数は44.1 kHzに設定されている。

 アナログ信号をデジタル信号に変換すると、サンプリング周波数の1/2を中心として本来の周波数成分とその鏡像に当たる周波数成分が対称に出現する。この鏡像をエイリアス(alias)と呼ぶ。




離散フーリエ変換(Discrete Fourier Transform)


 離散フーリエ変換は、フーリエ変換の時間を離散化したものである。  秒観察された連続的な波 に対して、  秒の一定間隔でサンプリングされたN 個の数値列(離散的な波)を

 

とする。この離散的な波をフーリエ変換はできないので、区分求積法を使って近似されたフーリエ変換 を求めることができる。ここで は周波数分解能(基本周波数)である。




Mathematicaで離散フーリエ変換をしてみよう


解析の手順

❶ 音源を取り込む

・Workフォルダにある、音源サンプルをホームディレクトリに置く。(デスクトップではないので注意)

・念のため、定義すべてをクリア

Clear["Global`*”]

・作業ディレクトリを設定する

SetDirectory[$UserDocumentsDirectory]

・「ディレクトリにあるファイル」のボタンを押して、音源ファイルがあることを確認する。

・音源を取り込んで、名前をつける。ここではguitar5.wavを取り入れ、sound1とした。

sound1 = Import["guitar5.wav”]


❷ 音源データの中にはいくつかデータがあるので、音の離散データ(数列)とサンプリング周波数を抜き出し、それぞれ名前をつける。(ステレオデータなら左右の音のデータがある)

・sound1の1番(左)の音のデータを取り出しsound1listとした。

sound1list = Part[sound1, 1, 1, 1]

(2番(右)データの場合は sound1listr = Part[sound1, 1, 1, 2] とする。 )

・サンプリング周波数をデータから抜き出し、samrateと名前をつけた。

samrate = Part[sound1, 1, 2]

・データ数の確認

Length[sound1list]


❸ データ数を減らす。サンプリング周波数と同じデータ数にすると正規化された値で出てくるので、そのデータ数だけ抜き出す。

・soundlist1の中から、64001番目から64000+samrate番目のデータを抜き出し、data1とした。

data1 = Take[sound1list, {64000 + 1, 64000 + samrate}]

・data1をプロットして確認(プロット範囲を変えてみてみる)

ListLinePlot[data1, PlotRange -> {{1, samrate}, {-0.3, 0.3}}]

・data1を鳴らして確認(SampleRate をsamrateより小さくするとどうなるか試してみる)

ListPlay[data1, SampleRate -> samrate]


❹ 離散フーリエ変換する。

・data1を離散フーリエ変換しcを名前をつけた。

c = Fourier[data1]

・振幅の絶対値(パワースペクトル)を見たいので、cの絶対値Absをとってグラフで確認

ListPlot[Abs[c], PlotRange -> {{0, 1600}, {0, Max[Abs[c]]}}, Joined -> True]

・横軸は(ほぼ)周波数なので、ピークになっている数値を見ると、含まれている周波数成分がわかる。実際には、周波数成分は0からで、データは1からになるので、1つずれる。

例:130, 259, 389, 777, 647にピークがあった場合、含まれている周波数成分は

129 Hz, 258 Hz, 388 Hz, 656 Hz, 776 Hz

である。


課題1 その他の音源も取り込んで周波数成分を解析してみよう。



サンプリング周波数>音のデータの場合


課題2 recorder.wavを取り込んで周波数分析してみよう。(途中までは上と同じ。)

56001個目のデータからsamrate個取り出してみると、3つ音の波形が入ってしまっている。

s1 = ListLinePlot[partdata2list,  PlotRange -> {{5000, 5000 + 3 44100/440}, {-1, 1}}]

s2= ListLinePlot[partdata2list, PlotRange -> {{13000, 13000 + 3 44100/440}, {-1, 1}}]

s3 = ListLinePlot[partdata2list, PlotRange -> {{25000, 25000 + 3 44100/440}, {-1, 1}}]

❶ このまま離散フーリエ変換すると、どの音の周波数成分なのかわからない。そこで、解析したい音の部分だけ取り出す。

・samrateのデータ数だけ取り出したpartdata2listの中から、1501番目から10320番目のデータを取り出し、data21と名前をつけた。

data21 = Take[partdata2list, {1501, 10320}]

・データ数確認

Length[data21]

❷ データ数がsamrateより少ないと、出てきた横軸が周波数にならない(少し計算が必要)。そこで、samrateと同じデータ数になるように調整する。調整する方法はたくさんあるが、いびつなつなぎ目を作らずに余分な周波数成分を入れないことが重要である。

・最も簡単な方法(残りのデータを0にする)

・足りない分だけ0のデータを用意する

point0 = Table[0, {4*8820}];

・data21のデータに0のデータをつなげる

data22 = Join[data21, point0]

・データ数確認

Length[data22]

・それぞれ音で確認

ListPlay[data21, SampleRate -> 44100]

ListPlay[data22, SampleRate -> 44100]

・data22を離散フーリエ変換(Chopは0に近い値を0にするコマンド)

data22f = Chop[Fourier[data22]]

・パワースペクトルの確認(3音入っていたものと比較してみよう)

ListPlot[Abs[data22f],  PlotRange -> {{0, 3600}, {0, Max[Abs[data22f]]}}, Joined -> True]