p = 1, 2, cos 2n + p)πj = cos 2nπj 2n + p)πj, sin = sin 2nπj 7.1) f j = a ) 0 + a p + a n+p cos 2nπj p=1 p=0 1 + ) b n+p p=0 sin 2nπj 1 2 a 0 +

(1)

7 離散フーリエ変換

音源データに離散フーリエ変換を施し，エネルギースペクトルの分布を調べます．エネルギースペクトルの分布によって楽器の音色が特徴付けられるのではないか，と思っています．

7.1 音源データ

音楽や音声などのアナログ波形は，一定時間ごとに標本を採ってデジタル化されます．１秒間に採る標本の数をサンプリング周波数といいます．ダウンロードしたフォルダ sound_wav_files の中に含まれる flu00.wav などの音源データ（.wav ファイル）のサンプリング周波数は 44.1 kHz です．つまり，これらの音源データはアナログ波形から 1/44100 秒ごとに標本を採取してデジタル化されたものです．音源データは spwave を使えばテキストファイルに変換できます．メニューバーの＜ファイル＞から＜名前を付けて保存＞を選びます．≪ファイルの種類 ≫を Text with Time にして適切な名称で保存すれば，第１列に時刻，第２列に音の強さが記録されたテキストファイルが作成されます．テキストファイルに記録された音は spwave を利用して聴くことができます． T > 0 を周期，∆t をサンプリング周波数の逆数とします．サンプリング周 波数が 44.1 kHz のときは ∆t = 1 44100 です．時刻 j∆t での音の強さを fj とすれば，音源データは {f0, f1, f2,· · · , fN−1} で表現されます（N = T ∆t）．

7.2 離散データの三角関数による表現

音源データのような周期 T の離散データ {f0, f1, f2,· · · , fN−1} を周期 T の連 続関数 f (t) に拡張しましょう．つまり， fj = f (j∆t) (j = 0, 1, 2,· · · , N − 1; ∆t = T N) となるように周期関数 f (t) を選びましょう．f (t) は (6.1) のようにフーリエ級 数展開 f (t) = 1 2a0+ ∞ ∑ n=1 ( ancos 2nπt T + bnsin 2nπt T ) されますが，格子点 j∆t（j = 0, 1, 2,· · · , N − 1）に限定すれば fj = f (j∆t) = 1 2a0+ ∞ ∑ n=1 ( ancos 2nπj∆t T + bnsin 2nπj∆t T ) = 1 2a0+ ∞ ∑ n=1 ( ancos 2nπj N + bnsin 2nπj N ) (7.1)

(2)

となります．ここで，p = 1, 2,· · · について cos2(n + N p)πj N = cos 2nπj N , sin 2(n + N p)πj N = sin 2nπj N であることに注目すれば，(7.1) は fj = 1 2a0+ ∞ ∑ p=1 aN p+ N∑−1 n=1 ( _∞ ∑ p=0 an+N p ) cos2nπj N + N∑−1 n=1 ( _∞ ∑ p=0 bn+N p ) sin2nπj N と書き直せます．1 2a0 + ∞ ∑ p=1 aN p を改めて a0 とおきましょう． ∞ ∑ p=0 a1+N p，· · · ， ∞ ∑ p=0 aN−1+Np を，それぞれ，改めて a1，· · · ，aN−1 とおき，さらに， ∞ ∑ p=0 b1+N p， · · · ， ∞ ∑ p=0 bN−1+Np も，それぞれ，改めて b1，· · · ，bN−1 とおけば，上式は fj = a0+ N∑−1 n=1 ( ancos 2nπj N + bnsin 2nπj N ) (j = 0, 1, 2,· · · , N − 1) (7.2) となります．

7.3 離散フーリエ変換

離散データ{f0, f1, f2,· · · , fN−1} に対して (7.2) を満たす a0，a1，a2，· · · ，aN−1， b1，b2，· · · ，bN−1 を対応させる変換を離散フーリエ変換と呼びます．演習７−１【離散フーリエ変換の基礎】 N を正の整数とし，θ = 2π/N と おく．また，i を虚数単位（i =√−1），n を整数，初項（第１項）が 1，公比 r = einθ の等比数列の初項から第 N 項までの和を S とする． (1) n = 0 のとき S を求めなさい． (2) n6= 0 のとき S を求めなさい．

(3) j を整数とするとき，cos njθ− cos (N − n)jθ と sin njθ + sin(N − n)jθ を 計算しなさい． (4) 複素数の列 f0，f1，f2，· · · ，fN−1 と整数 n に対して an = 1 N N∑−1 k=0 fkcos nkθ, bn= 1 N N∑−1 k=0 fksin nkθ とする．このとき an− aN−n と bn+ bN−n を計算しなさい．

(3)

{f0, f1, f2,· · · , fN−1} の離散フーリエ変換が                an = 1 N N−1 ∑ k=0 fkcos 2πnk N (n = 0, 1, 2,· · · , N − 1) bn = 1 N N−1 ∑ k=0 fksin 2πnk N (n = 1, 2,· · · , N − 1) (7.3) で与えられることを証明しましょう．【証明】θ = 2π/N とおきます．i を虚数単位（i =√−1）とし，cn= an− ibn

とおきましょう．cos nkθ− i sin nkθ = e−inkθ を利用すれば，(7.3) は

cn = 1 N N∑−1 k=0 fke−inkθ (n = 0, 1, 2,· · · , N − 1) (7.4) にまとめられます（b0 = 0）．また，このような複素表現を採用すれば，(7.2) の右辺は a0+ N∑−1 n=1 (ancos jnθ + bnsin jnθ) = N∑−1 n=0 (ancos jnθ + bnsin jnθ) = Re N∑−1 n=0 cneijnθ (j = 0, 1, 2,· · · , N − 1) となります．ここに，Re z は複素数 z の実部を表します． N∑−1 n=0 cneijnθ に (7.4) を代入すれば N∑−1 n=0 cneijnθ = N∑−1 n=0 ( 1 N N∑−1 k=0 fke−inkθ ) eijnθ = 1 N N∑−1 k=0 fk (_N₋₁ ∑ n=0 ei(j−k)nθ ) となります．ここで， N∑−1 n=0 ei(j−k)nθ は初項 1，公比 r = ei(j−k)θ の等比数列の初 項から第 N 項まで和ですから，演習７−１ (1)，(2) より N∑−1 n=0 cneijnθ = 1 N N∑−1 k=0 fk (_N₋₁ ∑ n=0 ei(j−k)nθ ) = fj が得られます．これは (7.3) が {f0, f1, f2,· · · , fN−1} の離散フーリエ変換であることを示しています．【証明終了】 以下では N は偶数であるとしましょう，N = 2M ．演習７−１ (3)，(4) よ り，n = 1, 2,· · · , M − 1 に対して an = aN−n, ancos 2nπj N = aN−ncos 2(N − n)πj N

(4)

bn=−bN−n, bnsin 2nπj N = bN−nsin 2(N− n)πj N となります．また，sin2M πj N = sin πj = 0 ですから，(7.2) は fj = (2a0) 2 + M∑−1 n=1 ( 2ancos 2nπj N + 2bnsin 2nπj N ) + aMcos 2M πj N (j = 0, 1, 2,· · · , N − 1) (7.5) とも表現されます．

7.4 エネルギースペクトル

三角関数の合成          2ancos 2nπj N + 2bnsin 2nπj N = 2sncos ( 2nπj N − ψn ) sn = √ an2+ bn2, cos ψn = an sn , sin ψn= bn sn を適用すれば，(7.5) は fj = (2a0) 2 + M∑−1 n=1 2sncos ( 2nπj N − ψn ) + aMcos 2M πj N (j = 0, 1, 2,· · · , N − 1) (7.6) と書き換えられます．2nπj N = 2π n T(j∆t)ですから，2snは波{f0, f1, f2,· · · , fN−1} における周波数 n T の強さを表します．定数分の強さ|2a0|，周波数 1 T， 2 T，· · · ， M − 1 T の強さ，および，周波数 M T の強さ|aM| を並べた {|2a0|, 2s1, 2s2, · · · , 2sM−1, |aM|} (7.7) をエネルギースペクトルと呼びます．

Figure 2.1 はこのようにして得られた ﬂute の octave 3 の Do のエネルギースペクトルを表します．同様に，Figures 2.2，2.3 は，それぞれ，harp と piano の octave 3 の Do のエネルギースペクトルを表します．エネルギースペクトルに各楽器の特徴が表れています．

7.5 高速フーリエ変換

離散フーリエ変換のプログラムは (7.3) をプログラム言語で表現すれば得られ ます．しかし，(7.3) そのままでは N が大きいときに計算時間が長くなってし まいます．例えば，サンプリング周波数 44.1 kHz，記録時間 0.5 秒のときは

(5)

N = 22050 になり，私のノートパソコンでは約 1 分 30 秒でようやく結果が得られます． 1965 年頃に N を 2 のべき乗（N = 2k，k は正の整数）に選び，三角関数 の性質を十分に利用して計算の順序を工夫する高速フーリエ変換（Fast Fourier Transform，略称：FFT）が提案されました．FFT の応用範囲は大変広いので多数の研究者が様々な FFT パッケージを作ってきました．みなさんには，大浦拓哉さん（京都大学数理解析研究所）が作成した FFT パッケージを利用したプログラムを提供します．ちなみに，NASA（アメリカ航空宇宙局）も大浦さんの FFT を使っているそうです．また，(7.3) そのままでは約 1 分 30 秒を要した計算も大浦さんの FFT では 2 ∼ 3 秒で終わってし まいました．大浦さんの FFT パッケージは http://www.kurims.kyoto-u.ac.jp/~ooura/fft-j.html からダウンロードできます．みなさんに提供するプログラム fft_h_order.c では簡易版である fft4g_h.c を利用しています． ≪ ﬀt h order.c の解説≫ このプログラムは，テキストファイルから離散データを読み込み，fft4g_h.c を利用して離散フーリエ変換を行います．さらに， エネルギースペクトルを計算し，n = 1, 2,· · · , M − 1 をエネルギースペクトル の降順にソーティングし，結果をファイルに書き込みます． (1) ファイル名の指定このプログラムは，拡張子が .txt であるテキストファイルに離散データが記録されていることを仮定しています．たとえば，ファイル名が flu04.txt の場合は flu04 と入力しなさい． (2) 離散データの読み込み データ数 N ，周期 T の値も判断します．周期性を 確定するために f0 と fN の平均を新たな f0 ≡ fN にします． (3) 離散データ数のチェック FFT を適用するためには，データ数 N は 2 の べき乗である必要があります．そうでないときは，N < 2k _{となる最小の整数 k} を求め，(4) の段階で 2k を新たなデータ数とします． (4) 離散データ数を 2 のべき乗にする N < 2k _の場合は fj = f0 （j = N + 1, N + 2,· · · , 2k）によって離散データを拡張します．また，周期とデータ数を新しいものに置き換えます．さらに，拡張された離散データをファイルに記録します．ファイル名は，たとえば (1) で flu04 と入力した場合 flu04_ext.txt になります． (5) FFT による離散フーリエ変換 rdft(N, 1, f) が FFT を行う関数です． rdft(N, 1, f)は 1367 行ありますが，大浦さんが作成した FFT パッケージの中で最も簡易なものです．rdft(N, 1, f) を呼ぶとき，配列 f[0]，f[1]，f[2]， · · · ，f[N-1] には離散データ f0，f1，f2，· · · ，fN−1を入れます．rdft(N, 1, f) から戻ってきたとき，配列 f[0]，f[1]，f[2]，· · · ，f[N-1] は (7.5) の a0，a1， a2，· · · ，aM，b1，b2，· · · ，bN−1 につぎのように対応します： f[2n] = N an（n = 0, 1, 2,· · · , M − 1），f[1] = NaM f[2n+1] = N bn（n = 1, 2,· · · , M − 1）．

(6)

/* fft_h_order.c */ #include <math.h> #include <stdio.h> #include <string.h> #define NMAX (32768) #define NMAXH (16387)

void rdft(int, int, double *);

//2^12 = 4096, 2^13 = 8192, 2^14 = 16384, 2^15 = 32768, 2^16 = 65536 int main()

{

int N, Ntmp, M, i, j, k, n, flag, order[NMAXH+1]; double T, f[NMAX+1], es[NMAXH+1];

FILE *fp;

char file_name[100], file_in[100], file_out[100]; double zero = 0.0, tmp;

// (1) ファイル名の指定

printf("filename = "); scanf("%s", file_name);

strcpy(file_in, file_name); strcat(file_in, ".txt"); fp = fopen(file_in, "r");

// (2) 離散データの読み込み j = 0; flag = 1;

while (1) {

flag = fscanf(fp, "%lf %lf", &T, &f[j]); if (flag == EOF) break;

j = j+1; } N = j - 1; fclose(fp); printf("ORIGINAL: N = %d T = %f\n", N, T); f[0] = (f[0] + f[N])/2.0; f[N] = f[0]; // (3) 離散データ数のチェック k = log(N)/log(2.0); if (pow(2, k) < N) k = k + 1; Ntmp = pow(2, k); // (4) 離散データ数を 2 のべき乗にする if (Ntmp > N) { for(j = N + 1; j <= Ntmp; j++) f[j] = f[0]; T = (T*Ntmp)/N; N = Ntmp; printf("ADDED: N = %d T = %f\n", N, T);

strcpy(file_out, file_name); strcat(file_out, "_ext.txt"); fp = fopen(file_out, "w"); for(j = 0; j <= N; j++) fprintf(fp, "%f %f\n", j*T/N, f[j]); fclose(fp); } // (5) FFT による離散フーリエ変換 rdft(N, 1, f); // a[n] = f[2*n]*2.0/N (n = 0, 1, 2, ..., M-1) // a[M] = f[1]/N // b[n] = f[2*n+1]*2.0/N (n = 1, 2, 3, ..., M-1) // b[0] = b[M] = 0

(7)

// (6) エネルギースペクトルの計算 M = N/2; es[0] = fabs(f[0])*2.0/N; es[M] = fabs(f[1])/N; for (n = 1; n < M; n++) { es[n] = sqrt(f[2*n]*f[2*n]+f[2*n+1]*f[2*n+1])*2.0/N; } // (7) エネルギースペクトルの大きさによる並べ替え for (n = 0; n <= M; n++) order[n] = n; for (n = 2; n < M; n++) {

for (k = 1; k < n; k++) if (es[k] < es[n]) break; tmp = es[n]; for (i = n; i > k; i--) { es[i] = es[i-1]; order[i] = order[i-1]; } es[k] = tmp; order[k] = n; } // (8) 計算結果のファイルへの記録

strcpy(file_out, file_name); strcat(file_out, "_FFT.txt"); fp = fopen(file_out, "w");

n = 0;

fprintf(fp, "%d %d %f %f %f %f %f %f\n", N, n, n/T, f[n]*2.0/N, zero, es[n], es[n]/es[1], T);

fprintf(fp, "%d %d %f %f %f %f %f\n", M, M, M/T, f[1]/N, zero, es[M], es[M]/es[1]);

for (k = 1; k < M; k++) { n = order[k]; fprintf(fp, "%d %d %f %f %f %f %f\n", k, n, n/T, f[2*n]*2.0/N, f[2*n+1]*2.0/N, es[k], es[k]/es[1]); } return 0; } (6) エネルギースペクトルの計算配列 es[0]，es[1]，es[2]，· · · ，es[M] と (7.7) の a0，s1，· · · ，sM−1，aM の対応はつぎの通りです．

es[0] =|2a0|，es[n] = 2sn（n = 1, 2,· · · , M − 1），es[M] = |aM|

(7)エネルギースペクトルの大きさによる並べ替え モード n = 1, 2,· · · , N − 1 をエネルギースペクトルが降順になるようにソーティングします． (8) 計算結果のファイルへの記録計算結果をファイルに記録します．ファイル名は，たとえば (1) で flu04 と入力した場合 flu04_FFT.txt になります． 1 行目にはデータ数 N と周期 T に加えて n = 0 に関する情報が記録され ます．左から順に N 0 0 T 2a0 0 |2a0| |2a0| 2smax T が記録されます．ここに，smax は s1，s2，· · · ，sN−1 の最大値です． 2 行目には n = M に関する情報が記録されます．左から順に

(8)

M M M T aM 0 |aM| |aM| 2smax が記録されます． 3 行目からはモード n = 1, 2,· · · , N − 1 に関する情報が記録されます．エ ネルギースペクトルの降順にモード n が並べられます．2s1，2s2，· · · ，2sM−1 の中で大きい方から k 番目のものが 2sn であれば，k+2 行目には，左から順に k n n T 2an 2bn 2sn 2sn 2smax が記録されます．3 列目の n T は周波数を表すことに注意しましょう． 0 0.2 0.4 0.6 0.8 1 0 2000 4000 6000 8000 10000 12000 14000 0 0.2 0.4 0.6 0.8 1 0 2000 4000 6000 8000 10000 12000 14000 0 0.2 0.4 0.6 0.8 1 0 2000 4000 6000 8000 10000 12000 14000 0 0.2 0.4 0.6 0.8 1 0 2000 4000 6000 8000 10000 12000 14000 0 0.2 0.4 0.6 0.8 1 0 2000 4000 6000 8000 10000 12000 14000

Figure 7.1: octave 4の Do のエネルギースペクトルの分布（左上：ﬂute，右上： harp，左中：piano，右中：tuba，下：violin）

以下の課題を行うためには，更新された csc10.tar をダウンロード・展開する必要があります（Linux）．新しいファイル flu04.txt，har04.txt，pia04.txt，

(9)

tub04.txt，vio04.txt は，それぞれ，flute，harp，piano，tuba，violin の octave 4 の Do の音源データ（音一つ分）をテキストファイルの形で記録したものです． fft_h_order.c の実行形式ファイル fft_h_order は ./run.sh とすることによって得られます． ® © ª 課題７−１（Linux）【ﬂute の octave 4 の Do エネルギースペクトル】 fft_h_orderを実行し，ファイル名として flu04 を入力すれば，flu04.txt に記録された音源データの高速フーリエ変換が行われ，フーリエ係数とエネルギースペクトルが計算されます．さらに，これらの計算結果は，エネルギースペクトルの降順に，flu04_FFT.txt に記録されます．

横軸を周波数，縦軸をエネルギースペクトルを 2smax で割ったものとする

Figure 7.1 左上のような ﬂute の octave 4 の Do エネルギースペクトル分布図を gnuplot を利用して描き，印刷しなさい．ただし，描画範囲は横軸：−300 ∼ 15000，縦軸：−0.1 ∼ 1.1 としなさい． ® © ª

課題７−２（Linux）（ＴＡチェック） harp，piano，tuba および violin の octave 4の Do エネルギースペクトル分布図を gnuplot を利用して描き，印刷しなさい（Figure 7.1 右上，左中，右中，下参照）．ただし，描画範囲は横軸： −300 ∼ 15000，縦軸：−0.1 ∼ 1.1 としなさい． 離散データ {f0, f1, f2,· · · , fN−1} に対して (7.2) を満たす a0，a1，a2，· · · ， aN−1，b1，b2，· · · ，bN−1 を対応させる変換 (7.3) を離散フーリエ変換と呼びま した．逆に，a0，a1，a2，· · · ，aM，b1，b2，· · · ，bM−1 に対して (7.2) や (7.5) や (7.6) で定まる{f0, f1, f2,· · · , fN−1} を対応させる変換を逆離散フーリエ変換 と呼びます（N = 2M ，an=aN−n，bn=−bN−n に注意）． (7.6) の右辺ではすべての n = 1, 2,· · · , M − 1 について和をとっています fj = (2a0) 2 + aM cos 2M πj N + M∑−1 n=1 2sncos ( 2nπj N − ψn ) (j = 0, 1, 2,· · · , N − 1) が，smax を s1，s2，· · · ，sN−1 の最大値とし，sn > smax/2 を満たす n につい てのみ和をとったものを fj[1] （j = 0, 1, 2,· · · , N − 1）としましょう： f_j[1] = (2a0) 2 + aMcos 2M πj N + ∑ sn>smax/2 2sncos ( 2nπj N − ψn ) . 同様に m = 2, 3,· · · として，sn > smax/2m を満たす n についてのみ和をとっ たものを fj[m] （j = 0, 1, 2,· · · , N − 1）としましょう： f_j[m] = (2a0) 2 + aM cos 2M πj N + ∑ sn>smax/2m 2sncos ( 2nπj N − ψn ) .

(10)

≪ idft.c の解説≫ サンプルプログラム idft.c は { f₀[m], f₁[m],· · · , f_N[m]₋₁ } が m が大きくなるにつれて (7.2) や (7.5) や (7.6) の{f0, f1,· · · , fN−1} に収束することを確認するために， { f₀[m], f₁[m],· · · , f_N[m]₋₁ } （m = 1, 2,· · · , 9）を計算し， ファイルに記録するものです． (1) ファイル名の指定離散フーリエ変換の結果が fft_h_order.c の (8) 計算結果のファイルへの記録で述べたものと同じ方法で記録されていることが仮定されています．たとえば，元々の離散データが記録されたファイル名が flu04.txt の場合は flu04_FFT と入力しなさい． (2)フーリエ係数とエネルギースペクトルの読み込み配列 a[0]，a[1]，a[2]，

· · · ，a[M]，b[1]，b[2]，· · · ，b[M-1] と (7.5) の a0，a1，a2，· · · ，aM，b1，b2，

· · · ，bM−1 の対応はつぎの通りです． a[n] = 2an（n = 0, 1, 2,· · · , M − 1），a[M] = aM b[n] = 2bn（n = 1, 2,· · · , M − 1） /* idfc.c */ #include <stdio.h> #include <math.h> #include <string.h> #define NMAX (32768) #define NMAXH (16387) int main () {

int N, M, i, j, k, n, flag, order[NMAXH+1];

double T, f[NMAX+1], a[NMAXH+1], b[NMAXH+1], es[NMAXH+1]; FILE *fp;

char file_name[100], file_in[100], file_out[100], file_add[100]; int ktmp1, ktmp2;

double dt, dtmp1, dtmp2, dtmp3, dtmp4, dtmp5; // (1) ファイル名の指定

printf("filename = "); scanf("%s", file_name);

strcpy(file_in, file_name); strcat(file_in, ".txt"); fp = fopen(file_in, "r");

// (2) フーリエ係数とエネルギースペクトルの読み込み fscanf(fp, "%d %d %lf %lf %lf %lf %lf %lf",

&N, &ktmp2, &dtmp1, &a[0], &dtmp2, &dtmp3, &dtmp4, &T); fscanf(fp, "%d %d %lf %lf %lf %lf %lf",

&ktmp1, &ktmp2, &dtmp1, &dtmp2, &dtmp3, &dtmp4, &dtmp5); M = ktmp1;

a[M] = dtmp2; k = 1; flag = 1; while (1) {

flag = fscanf(fp, "%d %d %lf %lf %lf %lf %lf",

&ktmp1, &n, &dtmp1, &dtmp2, &dtmp3, &dtmp4, &es[k]); if (flag == EOF) break;

order[k] = n; a[n] = dtmp2; b[n] = dtmp3; k = k+1; } fclose(fp);

(11)

// (3) フーリエ係数による関数の展開と結果のファイルへの記録 dt = T/N;

for (j = 0; j <= N; j++) f[j] = a[0]/2.0 + a[M]*cos(2.0*M*M_PI*j/N); k = 1; dtmp1 = 1.0; for (i = 0; i <= 8; i++) { dtmp1 = dtmp1/2.0; printf("i = %d: ", i, k); while (es[k]/es[1] >= dtmp1) { n = order[k]; for (j = 0; j <= N; j++) f[j] = f[j] + a[n]*cos(2.0*n*M_PI*j/N) + b[n]*sin(2.0*n*M_PI*j/N); k = k + 1; if (k > M) break; } strcpy(file_out, file_name);

sprintf(file_add, "_IDFT%1d%05d.txt", i+1, k-1); strcat(file_out, file_add);

fp = fopen(file_out, "w");

for (j = 0; j <= N; j++) fprintf(fp, "%f %f\n", j*dt, f[j]); fclose(fp);

printf("%s includes until %d\n", file_out, k-1); } return 0; } (3) フーリエ係数による関数の展開と結果のファイルへの記録たとえば，(1) でファイル名を flu04_FFT と指定した場合には， sn> smax/2である n がちょうど 9 個， sn> smax/22 である n がちょうど 21 個， sn> smax/23 である n がちょうど 40 個， · · · · sn> smax/29 である n がちょうど 827 個 あることが計算されます．このことに対応して { f₀[1], f₁[1],· · · , f_N[1] } がファイル flu04_FFT_IDFT100009.txt に， { f₀[2], f₁[2],· · · , f_N[2] } がファイル flu04_FFT_IDFT200021.txt に， { f₀[3], f₁[3],· · · , f_N[3] } がファイル flu04_FFT_IDFT300040.txt に， · · · · { f₀[9], f₁[9],· · · , f_N[9] } がファイル flu04_FFT_IDFT900827.txt に 記録されます．ただし，f_N[m] = f₀[m] としています（m = 1, 2,· · · , 9）． idft.c の実行形式ファイル idft は cc -O2 idft.c -lm -o idft

(12)

とすることによって得られます．

®

© ª

課題７−３（Linux & Windows）【ﬂute の octave 4 の Do の近似】 Linux で idft を実行し，ファイル名として flu04_FFT を入力すれば，ﬂute の octave 4 の Do がフーリエ級数展開によって近似されます．近似の結果は 9 つのファイル flu04_FFT_IDFT100009.txt，flu04_FFT_IDFT200021.txt，・・・， flu04_FFT_IDFT900827.txtに記録されます．この 9 つのファイルと元々の音源データ flu04.txt をフラッシュメモリーを介して Windows に移動し，フーリエ級数展開の収束の様子を spwave で音を聴くことによって確認しなさい． spwave で flu04_FFT_IDFT100009.txt などのテキストファイルを読む込むには，＜ファイル＞から＜開く＞を選び，ファイルを指定して≪開く≫ボタンを押した後，＜ファイルフォーマット＞を Text with Time に＜サンプリング周波数＞を44100 にして≪ＯＫ≫ボタンを押します．再生する前に，＜サウンド＞から＜ループ再生を有効に＞を選んでおくと便利でしょう． ® © ª

課題７−４（Linux & Windows）（ＴＡチェック） Linux で idft を実行し，ファイル名として har04_FFT を入力すれば，harp の octave 4 の Do がフーリエ級数展開によって近似されます．入力ファイルを pia04_FFT，tub04_FFT， vio04_FFTにすれば，それぞれ，piano，tuba，violin の octave 4 の Do がフーリエ級数展開によって近似されます．このようにして得られたすべてのファイルと元々の音源データをフラッシュメモリーを介して Windows に移動し，フーリエ級数展開の収束の様子を spwave で音を聴くことによって確認しなさい．以下のの課題に取り組む前に，sound wav ﬁles cats dogs をダウンロードしなさい（Windows）．その方法については第２章の≪音楽データ（楽器）のダウンロード≫を参照しなさい．

sound_wav_files や sound_wav_files_cats_dogs （Windows）に格納されている音源データ（.wavファイル）を解析しましょう．以下では音源データ flu02.wav （ﬂute の octave 3 の Mi が 9 音）を例として手順を説明します．注意１：(1)∼(3) および (8) は Windows で行う 注意２：(5)∼(6) は Linux で行う 注意３：(4)，(7) は Windows と Linux の間のデータの移動 (1) spwaveを起動し，flu02.wav を読み込み，再生する． (2) 音一つ分，あるいは，鳴き声一回分を切り出す．切り出す部分が 0.743 秒を越えないよう注意する．

(3)切り出した部分を別の名前で，Text with Time という種類のファイルとして保存する．拡張子は.txt とする．

(13)

(4) 作成したテキストファイルをフラッシュメモリーを介して Linux へ移動する． (5) fft_h_orderによって高速フーリエ変換（FFT）を行う．さらに，FFT の結果を利用してエネルギースペクトル分布図を描き，印刷する． (6) idft を利用して，音源データをフーリエ級数展開で近似する． (7) (6) で得られたファイルをフラッシュメモリーを介して Windows へ移動する． (8) フーリエ級数展開の収束の様子を，spwave で音を聴くことによって確認する． ® © ª

課題７−５（Windows，Linux & Windows）（ＴＡチェック） 指定された音源データに対して，上記の (1)∼ (8) を行いなさい．ただし，エネルギース ペクトル分布図を描く際の描画範囲は横軸：−300 ∼ 15000，縦軸：−0.1 ∼ 1.1 としなさい．扱うべき音源データを受講学生ごとにつぎのように指定します：学籍番号の下 1 桁が 0： flu01.wav，flu02.wav，flu03.wav の 3 個学籍番号の下 1 桁が 1： flu05.wav，flu06.wav，flu07.wav の 3 個学籍番号の下 1 桁が 2： har01.wav，har02.wav，har03.wav の 3 個学籍番号の下 1 桁が 3： har05.wav，har06.wav，har07.wav の 3 個学籍番号の下 1 桁が 4： pia01.wav，pia02.wav，pia03.wav の 3 個学籍番号の下 1 桁が 5： pia05.wav，pia06.wav，pia07.wav の 3 個学籍番号の下 1 桁が 6： tub01.wav，tub02.wav，tub03.wav の 3 個学籍番号の下 1 桁が 7： tub05.wav，tub06.wav，tub07.wav の 3 個学籍番号の下 1 桁が 8： vio01.wav，vio02.wav，vio03.wav の 3 個学籍番号の下 1 桁が 9： vio05.wav，vio06.wav，vio07.wav の 3 個 ® © ª

課題７−６（Windows，Linux & Windows）（ＴＡチェック）

sound_wav_files_cats_dogsに格納されている猫の鳴き声 cat1.wav，cat2.wav， cat3.wav，cat4.wav （cat1.wav∼ cat4.wav の内の 3 つは猫の鳴き声を録音したもの，1 つは人間が猫の鳴き声をまねたもの）に対して，上記の (1) ∼ (8) を行いなさい．ただし，エネルギースペクトル分布図を描く際の描画範囲は横

軸：−200 ∼ 10000，縦軸：−0.1 ∼ 1.1 としなさい．

なお，体験７−１で求められるレポートにエネルギースペクトル分布図を取り込みたい場合は，画像を保存する際に

set terminal postscript set output "filename.ps" replot

だけではなく，

(14)

課題７−７（Windows，Linux & Windows）（ＴＡチェック）

sound_wav_files_cats_dogs に格納されている犬の鳴き声のうち，指定されたものに対して，上記の (1) ∼ (8) を行いなさい．ただし，エネルギースペク トル分布図を描く際の描画範囲は横軸：−200 ∼ 10000，縦軸：−0.1 ∼ 1.1 とし なさい．扱うべき犬の鳴き声を受講学生ごとにつぎのように指定します：学籍番号の下 1 桁が 0： dog1.wav，dog2.wav，dog3.wav の 3 個学籍番号の下 1 桁が 1： dog4.wav，dog b1.wav，dog b2.wav の 3 個学籍番号の下 1 桁が 2： dog b3.wav，dog b4.wav，dog b5.wav の 3 個学籍番号の下 1 桁が 3： dog b6.wav，dog b7.wav，dog b8.wav の 3 個学籍番号の下 1 桁が 4： dog b9.wav，dog bss1.wav，dog bss2.wav の 3 個学籍番号の下 1 桁が 5： dog bss3.wav，dog bss4.wav，dog1.wav の 3 個学籍番号の下 1 桁が 6： dog2.wav，dog3.wav，dog4.wav の 3 個

学籍番号の下 1 桁が 7： dog b1.wav，dog b2.wav，dog b3.wav の 3 個学籍番号の下 1 桁が 8： dog b4.wav，dog b5.wav，dog b6.wav の 3 個学籍番号の下 1 桁が 9： dog b7.wav，dog b8.wav，dog b9.wav の 3 個なお，dog3.wav は成犬の鳴き声，dog4.wav は幼犬の鳴き声です．また dog bss1.wav

∼ dog bss4.wav は springer spaniel という犬種の鳴き声です．

®

© ª

体験７−１（Windows，Linux & Windows）【レポート提出】 主としてエネルギースペクトルの分布という視点から，猫の鳴き声と犬の鳴き声について 1000 字以上 2000 字以内で論じなさい．つぎのようなことがらに対する観察，解析と見解を述べることが期待されています．・猫の鳴き声に見られる共通点・犬の鳴き声に見られる共通点・猫と犬の鳴き声の共通点と相違点・本物の猫と人間が真似をした鳴き声の判別可能性レポートに図を入れることは大変良いことですが，字数としてはカウントされません．提出日時は１１月２２日１１時とします．