まず，テキスト336ページのFig. 13-4に示されているようなテキストファイルを，テキストエディタで作成して下さい。
その際，ファイル名が，確かにhw.datになっていることを確認してください。
よくある間違い（「先生，ちゃんとテキストとおりに書いているのに，プログラムが動きません！」の原因）は，ファイル名が，hw.dat.txtなどになってしまっていることです。MacOSやWindowsなどのGUIのOSは，ファイル名の拡張子（ . の後ろの部分）を表示しなかったり，表示したり，勝手にひっつけたり，とっぱらったり，など，（コンピュータの初心者用に）とにかく余計なことをしくさる（笑）ので，注意が必要です。
（コンピュータに我々が振り回されるのではなく，我々のほうがコンピュータを完璧にコントロールしてやりましょう。）
ファイル名が，確かにhw.datになっていることを間違いなく確認するためには，ターミナルで，ls コマンドを使って，確認するのが確実です。
（GUIよりもCUIのほうが優れている点のひとつは，このような点です。）
そして，テキスト337ページの List13-2 の動作を確認してソースプログラムを送信してください。
repo14-1にアップロードしてください。
(重要) 出席登録のための実習問題は，授業時間中に必ず提出して下さい。
音響信号処理に関連したプログラムを作ってみましょう。
下図は，周期が２πとなる矩形波を表しています。
このような矩形波の，フーリェ級数展開は以下の式で表現されます。(ホント?と思った人は，フーリェ係数を実際に紙と鉛筆で計算してみましょう。)
では本当に上式によって，上図のような矩形波が表せるのでしょうか? 実際に計算して確かめてみましょう。
フーリェ級数展開は，たくさんの(理論的には無限個の)項数を足し合わせないといけませんが，このような作業こそ，プログラミング言語の得意とするところなのです。
それでは，
◯ 上式を使って，引数 x に対する f(x) の値を計算し，それを戻り値として返す関数を作成してみましょう。
◯ その関数を使って，0 <= x <= ２π の値域における f(x) の値を計算し，その結果をファイルへ書き込むようなプログラムに仕上げて下さい。
その際，サンプリング間隔は，0.01として下さい。(つまり，xの値を0.01ずつ増やしていって，f(x)を上の式で計算する。)
また，101次までの成分を使って計算して下さい。(つまり，sin(x) の項から，sin(101x) の項までを足し合わせる。)
計算結果は，1行ごとに，
xの値 (スペース) f(x)の値 (改行)
と組にして，ファイルへ書き込むようにして下さい。
----例：以下のようなテキストファイルができあがるようにして下さい----
0.000000 0.000000
0.010000 0.612974
0.020000 1.033412
0.030000 1.178326
0.040000 1.109464
0.050000 0.974797
0.060000 0.904130
0.070000 0.934903
0.080000 1.014467
0.090000 1.064231
0.100000 1.047940
.
.
.
このプログラムを実行してできあがったテキストファイルを，Excelなどのグラフ作成ソフトに読み込み，計算した f(x) をグラフにして見てみましょう。
“ギブスの現象”が確認できますか?
ソースプログラムを，repo14-2にアップロードしてください。
上記のプログラムは，フーリェ級数を計算するプログラムでしたが，
それではフーリェ変換を計算するプログラムのほうも，作ってみましょう。
離散時間信号のフーリェ変換については，以下のとおり。
上記を参考にして，
音響信号をファイルから読み込み，
その音響信号のフーリェ変換を計算し，
その計算結果をファイルに書き出す
プログラムを作ってみてください。以下の手順に従うとよいでしょう。
◯ テキストファイル：1kHzSine_1kTap_44100Hz.txt に保存されている音響信号を読み込むコードを，以下のようにして（List13-2をマネして）main( )関数内に書いてみましょう。
short tap;
double fs, T, Nf, time, data, x[15000];
char fname[64];
FILE *fp;
printf("時間波形のファイル名を入力して下さい：");
scanf("%s", fname);
printf("%sから読み込みます\n", fname);
printf("標本化周波数を入力して下さい：");
scanf("%lf", &fs);
printf("標本化周波数は%f Hzです\n", fs);
T = 1./fs; /* 標本化周期[s]*/
Nf = fs/2.; /* ナイキスト周波数[Hz] */
if( (fp = fopen(fname, "r")) == NULL ){
printf("\aファイルをオープンできません。\n");
}
else{
tap = 0;
while( fscanf(fp, "%lf%lf", &time, &data) == 2 ){
x[tap] = data;
tap++;
}
printf("時間波形の長さは%hd tapsです\n", tap);
fclose(fp);
}
◯ 離散時間信号のフーリェ変換 C(ω) の実数部分を計算し，それを戻り値として返す関数を，以下のようにして作成してみましょう。
double dftReal(double x[ ], short tap, double omega, double T)
/* 離散時間信号のフーリエ変換の実部を計算する関数
x[ ]：離散時間信号，tap：離散時間信号の長さ[taps]，omega：角周波数[1/s]，T：標本化周期[s] */
{
.
.
.
}
◯ 同様にして，離散時間信号のフーリェ変換 C(ω) の虚数部分を計算し，それを戻り値として返す関数を作成してみましょう。
◯ それらの関数を使って，0 <= frequency <= fs / 2 [Hz] の値域における フーリェ変換 C(ω) の実数部分と虚数部分 の値を計算し，それらの値からエネルギスペクトル（ 2*|C(ω)|2 ）を計算して，その結果をファイルへ書き込むようなプログラムに仕上げて下さい。
その際，周波数間隔は，1 [Hz]として下さい。(つまり，frequencyの値を1ずつ増やしていって，エネルギスペクトルを計算する。)
計算結果は，1行ごとに，
frequencyの値 (スペース) エネルギスペクトルの値 (改行)
と組にして，ファイルへ書き込むようにして下さい。
----例：以下のようなテキストファイルができあがるようにして下さい
（以下は1kHzSine_1kTap_44100Hz.txtを読み込んで，標本化周期 T = 1/44100 [s] として，エネルギスペクトルを計算した例）----
0.000000 224.953707
1.000000 223.957228
2.000000 220.987465
3.000000 216.103058
4.000000 209.400450
5.000000 201.011994
6.000000 191.103337
7.000000 179.870161
8.000000 167.534318
9.000000 154.339465
10.000000 140.546251
.
.
.
例えば，以下のようなコードをmain( )関数内に書いてみるとよいでしょう。
double frequency, omega, cReal, cImag, energy;
printf("フーリエ変換のファイル名を入力して下さい：");
scanf("%s", fname);
printf("%sに保存します\n", fname);
if( (fp = fopen(fname, "w")) == NULL ){
printf("\aファイルをオープンできません。\n");
}
else{
frequency = 0.;
while(frequency<=Nf){ /* frequency：周波数を，0 Hzから，Nf：ナイキスト周波数[Hz]まで */
omega = 2.*PAI * frequency; /* 角周波数 */
cReal = dftReal(x, tap, omega, T); /* フーリエ変換の実部 */
cImag = dftImag(x, tap, omega, T); /* フーリエ変換の虚部 */
energy = 2.*(cReal*cReal + cImag*cImag); /* エネルギスペクトル */
fprintf(fp, "%f %f\n", frequency, energy);
frequency += 1.; /* freqを1 Hz増やす */
}
fclose(fp);
}
さて，プログラムが完成したら，それを使って，実際に音響信号の周波数分析をしてみましょう。
◯ まず最初に，1 [kHz] の正弦波を周波数分析してみましょう。
さきのテキストファイル：1kHzSine_1kTap_44100Hz.txt は，実は1 [kHz]の正弦波を，標本化周波数44100 [Hz]でサンプルした数値を，1000 [taps]の長さのぶん，書き込んだテキストファイルになっています。
それではこのプログラムを実行してできあがったテキストファイルを，グラフ作成ソフトに読み込み，計算した エネルギスペクトル をグラフにして見てみましょう。以下のようなグラフが描けるでしょうか？
◯ 次に，「ア」の音声の周波数分析をしてみましょう。
音響信号のテキストファイル：a_15kTap_44100Hz.txt を読み込むテキストファイルとして使用してください。
これは，音声「ア」を，標本化周波数44100 [Hz]でサンプルした数値を，15000 [taps]の長さのぶん，書き込んだテキストファイルになっています。
実際に音を聞いてみたいひとは，WAVEサウンドファイル：a_15kTap_44100Hz.wav を再生してみてください。
同じようにこのプログラムを実行してできあがったテキストファイルを，グラフ作成ソフトに読み込み，計算した エネルギスペクトル をグラフにして見てみましょう。以下のようなグラフが描けるでしょうか？
周波数分析がC言語で，できるようになりましたか？
|
最近は便利なソフト（しかし高価なもの，10数万円もする！）がたくさんあって，周波数分析（をはじめとして他のいろいろな科学技術計算）をするソフトをC言語で一から作らなくてもよい世の中になってきていますが， 皆さんが将来，たとえそのようなソフトが買えない貧乏・どん底状態になったとしても， (1) そこらへんに捨てられてあるパソコンを拾ってきて， (2) お昼ごはんを一回がまんして，例えば「自分でできるLinuxサーバ構築 CD-ROM付き」（仮名，古本で約170円！），みたいな本を買って， (3) Linux（UNIXの一種）をインストールし， (4) C言語で科学技術計算をするソフトを自分で作って， (5) いろいろな技術を独力で開発・駆使できるような， “本物の技術者” になってほしいと，私は思っています！！！ |
ソースプログラムを，repo14-3にアップロードしてください。
（上記のような分析結果グラフではなくて，ソースプログラムを，提出してください。）
(重要) 今日の宿題その1とその2の締切は，2月19日(金)までとします。じっくりと考えて，いいものを作って下さいね。(これらの宿題の解答例は，3/1になったらこのweb上に掲示します。)