6
離散フーリエ変換
音源データに離散フーリエ変換を施し，エネルギースペクトルの分布を調べま す．エネルギースペクトルの分布によって楽器の音色が特徴付けられるのでは ないか，と思っています．6.1
音源データ
音楽や音声などのアナログ波形は，一定時間ごとに標本を採ってデジタル化さ れます．１秒間に採る標本の数をサンプリング周波数といいます．ダウンロー ドしたフォルダ sound_wav_files の中に含まれる flu00.wav などの音源デー タ（.wav ファイル）のサンプリング周波数は 44.1 kHz です．つまり，これら の音源データはアナログ波形から 1/44100 秒ごとに標本を採取してデジタル化 されたものです． 音源データは spwave を使えばテキストファイルに変換できます．メニュー バーの＜ファイル＞から＜名前を付けて保存＞を選びます．≪ファイルの種類 ≫を Text with Time にして適切な名称で保存すれば，第１列に時刻，第２列 に音の強さが記録されたテキストファイルが作成されます．テキストファイル に記録された音は spwave を利用して聴くことができます． ∆t をサンプリング周波数の逆数，T = N · ∆t を周期とします（N は正の 整数）．サンプリング周波数が 44.1 kHz のときは ∆t = 1 44100 sec です．時刻 j∆t での音の強さを fj とすれば，音源データは {f0, f1, f2,· · · , fN−1} で表現 されます．6.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) は (5.1) のようにフーリエ級 数展開 f (t) = 1 2a0+ ∞ ∑ n=1 ( ancos 2nπt T + bnsin 2nπt T ) されますが，格子点 j∆t に限定すれば次のようになります： 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 ) (6.1)6.3
離散フーリエ変換
p = 1, 2,· · · について，2(n + N p)πj N = 2nπj N + 2pπj より， cos2(n + N p)πj N = cos 2nπj N , sin 2(n + N p)πj N = sin 2nπj N であることに注目しましょう．(6.1) の右辺に現れる n を N で割ったときの余 りが等しいものにグループ分けすれば，(6.1) は fj = 1 2a0+ ∞ ∑ p=1 aN p + ( ∞ ∑ p=1 aN p+1 ) cos2πj N + ( ∞ ∑ p=1 bN p+1 ) sin2πj N + ( ∞ ∑ p=1 aN p+2 ) cos2· 2πj N + ( ∞ ∑ p=1 bN p+2 ) sin2· 2πj N + · · · · + ( ∞ ∑ p=1 aN p+N−1 ) cos2(N− 1)πj N + ( ∞ ∑ p=1 bN p+N−1 ) sin2(N− 1)πj N と書き直せます．1 2a0+ ∞ ∑ p=1 aN pを改めて a0とおきましょう． ∞ ∑ p=0 aN p+1， ∞ ∑ p=0 aN p+2， · · · ， ∞ ∑ p=0 aN p+N−1 を，それぞれ，改めて a1，a2，· · · ，aN−1 とおき，さらに， ∞ ∑ p=0 bN p+1， ∞ ∑ p=0 bN p+2，· · · ， ∞ ∑ p=0 bN p+N−1も，それぞれ，改めて b1，b2，· · · ，bN−1 とおけば，上式は fj = a0+ N∑−1 n=1 ( ancos 2nπj N + bnsin 2nπj N ) (j = 0, 1, 2,· · · , N − 1) (6.2) となります．離散データ{f0, f1, f2,· · · , fN−1} に対して (6.2) を満たす a0，a1， a2，· · · ，aN−1，b1，b2，· · · ，bN−1 を対応させる変換を離散フーリエ変換と呼 びます． 演習６−１【離散フーリエ変換の基礎】 N を正の整数として θ = 2π/N とおく． (1) i を虚数単位（i = √−1），p を −N + 1 以上，N − 1 以下の整数とする．初項（第１項）が 1，公比 r = eipθ の等比数列の初項から第 N 項までの和を S
とする．p = 0 のとき S を求めなさい．また，p̸= 0 のとき S を求めなさい．
(2) 整数 n，j に対して
cos njθ− cos (N − n)jθ, sin njθ + sin(N − n)jθ をそれぞれ計算しなさい． (3) 複素数の列 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 を計算しなさい． 離散フーリエ変換の計算N を正の整数とします．{f0, f1, f2,· · · , fN−1} の離散 フーリエ変換が an = 1 N N∑−1 k=0 fkcos 2nπk N = 1 N N∑−1 k=0 fkcos nkθ (n = 0, 1, 2,· · · , N − 1) bn = 1 N N∑−1 k=0 fksin 2nπk N = 1 N N∑−1 k=0 fksin nkθ (n = 1, 2,· · · , N − 1) (6.3) で与えられることを証明しましょう．ここに θ = 2π N です． （証明）(6.3) を (6.2) の右辺に代入すると fj になることを示しましょう．(6.3) に は現れませんが，b0 = 0 とします．i を虚数単位（i = √ −1）とし，cn= an−ibn とおきましょう．すると， cn = 1 N N∑−1 k=0 fk(cos nkθ− i sin nkθ) = 1 N N∑−1 k=0 fk(cos(−nkθ) + i sin(−nkθ)) = 1 N N∑−1 k=0 fke−inkθ (6.4) です．(6.2) の右辺は a0+ N∑−1 n=1 (ancos njθ + bnsin njθ) = N∑−1 n=0 (ancos njθ + bnsin njθ) であり，
cneinjθ = (an− ibn)(cos njθ + i sin njθ)
= (ancos njθ + bnsin njθ) + i(ansin njθ− bncos njθ)
ですから (6.2) の右辺 = N−1∑ n=0 Re(cneinjθ ) = Re (N−1 ∑ n=0 cneinjθ )
が成り立ちます．ここに，Re (z) は複素数 z の実部を表します． N∑−1 n=0 cneinjθ は (6.4) を利用して， N∑−1 n=0 cneinjθ = N∑−1 n=0 ( 1 N N∑−1 k=0 fke−inkθ ) einjθ = 1 N N∑−1 k=0 fk (N−1 ∑ n=0 ei(j−k)nθ ) と計算されます．ここで， N∑−1 n=0 ei(j−k)nθ は初項 1，公比 r = ei(j−k)θ の等比数列 の初項から第 N 項まで和です．また，−N + 1 ≤ j − k ≤ N − 1 ですから，演 習６−１ (1) より N∑−1 n=0 cneijnθ = 1 N N∑−1 k=0 fk (N−1 ∑ n=0 ei(j−k)nθ ) = fj です．これは，(6.2) の右辺が fjであること，すなわち，(6.3) が{f0, f1, f2,· · · , fN−1} の離散フーリエ変換であることを示しています．（証明終了）
6.4
エネルギースペクトル
以下では N は偶数である（N = 2M ）として，(6.2) に関する解説を続けましょ う．演習６−１ (2)，(3) より，n = 1, 2,· · · , M − 1 に対して an = aN−n, ancos 2nπj N = aN−ncos 2(N − n)πj N bn=−bN−n, bnsin 2nπj N = bN−nsin 2(N− n)πj N となります．また，sin2M πj N = sin πj = 0 ですから，(6.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) (6.5) とも表現されます．さらに，三角関数の合成 2ancos 2nπj N + 2bnsin 2nπj N = 2sncos ( 2nπj N − ψn ) sn = √ an2+ bn2, cos ψn = an sn , sin ψn= bn sn を適用すれば，(6.5) は fj = (2a0) 2 + M∑−1 n=1 2sncos ( 2nπj N − ψn ) + aMcos 2M πj N (j = 0, 1, 2,· · · , N − 1) (6.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|} (6.7) をエネルギースペクトルと呼びます． 0 1000 2000 3000 4000 5000 0 1000 2000 3000 4000 5000
Figure 6.1: flute の octave 4 の Sol のエネルギースペクトル
0 1000 2000 3000 4000 5000 0 1000 2000 3000 4000 5000
Figure 6.2: harp の octave 4 の Sol のエネルギースペクトル
スペクトルを表します．同様に，Figures 6.2，6.3 は，それぞれ，harp と piano の octave 4 の Sol のエネルギースペクトルを表します． 0 1000 2000 3000 4000 5000 0 1000 2000 3000 4000 5000
Figure 6.3: piano の octave 4 の Sol のエネルギースペクトル
6.5
高速フーリエ変換
離散フーリエ変換のプログラムは (6.3) をプログラム言語で表現すれば得られ ます．しかし，(6.3) そのままでは N が大きいときに計算時間が長くなってし まいます．例えば，サンプリング周波数 44.1 kHz，記録時間 0.5 秒のときは N = 22050 になり，私のノートパソコンでは約 1 分 30 秒でようやく結果が得 られます． 1965 年頃に N を 2 のべき乗（N = 2k，k は正の整数）に選び，三角関数 の性質を十分に利用して計算の順序を工夫する高速フーリエ変換（Fast Fourier Transform，略称：FFT）が提案されました．FFT の応用範囲は大変広いので 多数の研究者が様々な FFT パッケージを作ってきました． みなさんには，大浦拓哉さん（京都大学数理解析研究所）が作成した FFT パッケージを利用したプログラムを提供します．ちなみに，NASA（アメリカ 航空宇宙局）も大浦さんの FFT を使っているそうです．また，(6.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 を利用しています． ≪ fft h order.c の解説≫ このプログラムは，テキストファイルから離散デー タを読み込み，fft4g_h.c を利用して離散フーリエ変換を行います．さらに， エネルギースペクトルを計算し，エネルギースペクトルの降順にソーティングし，結果をファイルに書き込みます． (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] は (6.5) の a0，a1， a2，· · · ，aM，b1，b2，· · · ，bN−1 につぎのように対応します： f[2n] = Nan（n = 0, 1, 2,· · · , M − 1），f[1] = NaM f[2n+1] = Nbn（n = 1, 2,· · · , M − 1）． (6) エネルギースペクトルの計算 配列 es[0]，es[1]，es[2]，· · · ，es[M] と (6.7) の a0，s1，· · · ，sM−1，aM の対応はつぎの通りです．
es[0] =|2a0|，es[n] = 2sn（n = 1, 2,· · · , M − 1），es[M] = |aM|
(7) エネルギースペクトルの大きさによる並べ替え エネルギースペクトルが 降順になるようにソーティングします． /* 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 // (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; } (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 に関する情報が記録されます．左から順に 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 は周波数を表す ことに注意しましょう． ≪実行形式ファイルの作成と実行≫ fft_h_order.c の実行形式ファイル fft_h_order は ./run.sh とすることによって得られます．fft_h_order が必要とする離散データを記録し たファイルとして，flu04.txt，har04.txt，pia04.txt，tub04.txt，vio04.txt が ~/csc に格納されています．これらは，それぞれ，flute，harp，piano，tuba， violin の octave 4 の Do の音源データ（音一つ分）をテキストファイルの形 で記録したものです．
課題６−１（Linux）【flute の octave 4 の Do エネルギースペクトル】 fft_h_order を実行し，ファイル名として flu04 を入力すれば，flu04.txt に 記録された音源データの高速フーリエ変換が行われ，フーリエ係数とエネルギー スペクトルが計算されます．さらに，これらの計算結果は，エネルギースペク トルの降順に，flu04_FFT.txt に記録されます． 0 500 1000 1500 2000 2500 0 2000 4000 6000 8000 10000 12000 14000 0 200 400 600 800 1000 1200 1400 1600 1800 0 2000 4000 6000 8000 10000 12000 14000 0 200 400 600 800 1000 1200 1400 1600 1800 2000 0 2000 4000 6000 8000 10000 12000 14000 0 200 400 600 800 1000 1200 1400 1600 0 2000 4000 6000 8000 10000 12000 14000 0 500 1000 1500 2000 2500 3000 0 2000 4000 6000 8000 10000 12000 14000
Figure 6.4: octave 4 の Do のエネルギースペクトルの分布（左上：flute，右上： harp，左中：piano，右中：tuba，下：violin）
横軸を周波数，縦軸をエネルギースペクトルとする Figure 6.4 左上のよう な flute の octave 4 の Do エネルギースペクトル分布図を gnuplot を利用して
課題６−２（Linux）（ＴＡチェック） harp，piano，tuba および violin の octave 4 の Do エネルギースペクトル分布図を gnuplot を利用して描き，印 刷しなさい（Figure 6.4 右上，左中，右中，下参照）．ただし，横軸の範囲は −300 ∼ 15000 としなさい．
6.6
逆離散フーリエ変換
離散データ{f0, f1, f2,· · · , fN−1} に対して (6.2) を満たす a0，a1，a2，· · · ，aN−1，
b1，b2，· · · ，bN−1 を対応させる変換 (6.3) を離散フーリエ変換と呼びました．
逆に，a0，a1，a2，· · · ，aM，b1，b2，· · · ，bM−1 に対して (6.2)，(6.5) や (6.6)
で定まる {f0, f1, f2,· · · , fN−1} を対応させる変換を 逆離散フーリエ変換 と呼 びます（N = 2M ，an=aN−n，bn=−bN−n に注意）． (6.6) に現れる s1，s2，· · · ，sN−1 の最大値を smax とします．(6.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) としています．ここでは，すべての n = 1, 2,· · · , M − 1 ではなく，sn> smax/2 を満たす n だけについてのみ和をとったものを f[1] j （j = 0, 1, 2,· · · , N − 1） としましょう．すなわち， fj[1] = (2a0) 2 + aMcos 2M πj N + ∑ sn>smax/2 2sncos ( 2nπj N − ψn ) とします．同様に m = 2, 3,· · · に対して，sn> smax/2m を満たす n について のみ和をとったものを f[m] j （j = 0, 1, 2,· · · , N − 1）としましょう．すなわち， fj[m] = (2a0) 2 + aMcos 2M πj N + ∑ sn>smax/2m 2sncos ( 2nπj N − ψn ) とします．m が大きくなるにつれて，fj[m] は fj に収束することが期待されま すが，この事実を確認・体験するプログラムが idft.c です． ≪ idft.c の解説≫ サンプルプログラム idft.c は {f0[m], f1[m],· · · , fN[m]−1 } が m が大きくなるにつれて (6.2)，(6.5) や (6.6) の {f0, f1,· · · , fN−1} に収束す ることを確認するために，{f0[m], f1[m],· · · , fN[m]−1 } （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] と (6.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); // (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 個 あることが計算されます．このことに対応して {f0[1], f1[1],· · · , fN[1] } がファイル flu04_FFT_IDFT100009.txt に， {f0[2], f1[2],· · · , fN[2] } がファイル flu04_FFT_IDFT200021.txt に， {f0[3], f1[3],· · · , fN[3] } がファイル flu04_FFT_IDFT300040.txt に， · · · · {f0[9], f1[9],· · · , fN[9] } がファイル flu04_FFT_IDFT900827.txt に 記録されます．ただし，f[m] N = f [m] 0 としています（m = 1, 2,· · · , 9）． idft.c の実行形式ファイル idft は cc idft.c -lm -o idft とすることによって得られます．
課題６−３（Linux & Windows）【flute の octave 4 の Do の近似】 Linux で idft を実行し，ファイル名として flu04_FFT を入力すれば，flute の 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 で音を聴くことによって確認しなさい．
6.7
離散フーリエ変換のまとめ
Windows の Q: ドライブのフォルダ csc の中には音源データ（.wav ファイル） を格納したフォルダ sound_wav_files，sound_wav_files_cats_dogs があり ます．sound_wav_files には楽器の音が，sound_wav_files_cats_dogs には 猫や犬の鳴き声が記録されています．flu02.wav（flute の 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 とする． (4) 作成したテキストファイルをフラッシュメモリーを介して Linux へ移動 する． (5) fft_h_order によって高速フーリエ変換（FFT）を行う．さらに，FFT の 結果を利用してエネルギースペクトル分布図を描き，印刷する． (6) idft を利用して，音源データをフーリエ級数展開で近似する． (7) (6) で得られたファイルをフラッシュメモリーを介して Windows へ移動 する．
(8) フーリエ級数展開の収束の様子を，spwave で音を聴くことによって確認 する．
課題６−５（Windows，Linux & Windows）（ＴＡチェック） 指定され
た音源データに対して，上記の (1) ∼ (8) を行いなさい．ただし，エネルギー スペクトル分布図を描く際の横軸の範囲は −300 ∼ 15000 としなさい．扱うべ き音源データを受講学生ごとにつぎのように指定します： 学籍番号の下 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 としなさい．
なお，体験６−１で求められるレポートにエネルギースペクトル分布図を取 り込みたい場合は，画像を保存する際に
set terminal postscript set output "filename.ps" replot
だけではなく，
set terminal postscript eps color set output "filename.eps"
replot ともしておいたら良いかも知れません．ここに，filename には自分で決める名 称を入れます．
課題６−７（Windows，Linux & Windows）（ＴＡチェック）
sound_wav_files_cats_dogs に格納されている犬の鳴き声のうち，指定され
トル分布図を描く際の横軸の範囲は −200 ∼ 10000 としなさい．扱うべき犬の 鳴き声を受講学生ごとにつぎのように指定します：
学籍番号の下 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 字以内で論じなさい．つぎのようなことがらに対する観 察，解析と見解を述べることが期待されています． ・猫の鳴き声に見られる共通点 ・犬の鳴き声に見られる共通点 ・猫と犬の鳴き声の共通点と相違点 ・本物の猫と人間が真似をした鳴き声の判別可能性 レポートに図を入れることは大変良いことですが，字数としてはカウントさ れません．提出日時は１１月 日１１時０５分 とします．
