/*  -*- MaTX -*-
 *
 *  【名前】
 *      sig_demo() - 信号処理ツールボックスのデモンストレーション
 *
 *  【形式】
 *      sig_demo()
 *
 */

Func void sig_demo()
{
	Integer n;
	List demos;

	void demo_spectral_analysis();
	void demo_fourier_expansion();
	void demo_convolution();


	demos = {"<<<  信号処理ツールボックス  >>>",
			 "FFT によるスペクトル解析",
			 "正弦波による矩形波の合成",
			 "畳み込みの(パーシバル)定理",
			 "終  了"};

	n = 1;
	while (1) {
		n = menu(demos, n);
		switch (n) {
		  case 1:
			demo_spectral_analysis();
			break;
		  case 2:
			demo_fourier_expansion();
			break;
		  case 3:
			demo_convolution();
			break;
		  default:
			return;
		}
		n++;
	}
}

Func void demo_spectral_analysis()
{
	Integer win;
	Array Y, x, y, Pyy, f, t;

	win = mgplot_cur_win();

	clear;
	print "\n";
	print "この例題は fft() を用いたスペクトル解析の方法を紹介する。\n";
	print "FFT の通常の目的は，雑音に埋められた信号の周波数成分を見付ける\n";
	print "ことである。\n\n";
	pause;

	print "まず，データを作成する。1000 Hzでサンプルされたデータを考える。\n";
	print "時刻 0 から 0.25 秒まで 1 ミリ秒毎の時間軸を作成することから";
	print "始めよう。\n\n";
	print "\tt = [0:0.001:0.25];\n\n";
	pause;
	t = [0:0.001:0.25];

	print "\n";
	print "次に，40 Hz と 130 Hz の信号の和を作る。\n\n";
	print "\tx = sin(2*PI*40*t) + sin(2*PI*130*t);\n\n";
	pause;
	x = sin(2*PI*40*t) + sin(2*PI*130*t);

	print "\n";
	print "そして，標準偏差 2 の白色雑音を信号 x に加え，信号 y とする。\n\n";
	print "\ty = x + 2*randn(t);\n\n";
	pause;
	y = x + 2*randn(t);

	print "\n";
	print "雑音で汚された信号 y をプロットしよう。\n\n";
	print "\tmgplot_title(win,\"Noisy time domain signal\");\n";
	print "\tmgplot_xlabel(win, \"time [s]\");\n";
	print "\tmgplot(win, t(1:50), y(1:50), {\"noisy signal\"});\n\n";
	pause;

	mgplot_reset(win);
	mgplot_title(win,"Noisy time domain signal");
	mgplot_xlabel(win, "time [s]");
	mgplot(win,t(1:50),y(1:50),{"noisy signal"});

	clear;
	print "\n";
	print "明らかに，この信号から元の周波数成分を同定するのは難しい。\n";
	print "このためスペクトル解析が有効である。\n\n";

	print "雑音で汚された信号 y の離散フーリエ変換は簡単に求まる。\n";
	print "次のように高速フーリエ変換(FFT)すれば良い。\n\n";
	print "\tY = fft(y, 256);\n\n";
	pause;
	Y = fft(y,256);

	print "\n";
	print "それぞれの周波数におけるパワースペクトル密度(エネルギーの尺度)\n";
	print "は次のように求まる。\n\n";
	print "\tPyy = Re(Y .* conj(Y));\n\n";
	pause;
	Pyy = Re(Y .* conj(Y));

	print "\n";
	print "パワースペクトル密度をプロットするため周波数軸を作る。\n\n";
	print "\tf = 1000*[0:127]/256;\n\n";
	pause;
	f = 1000*[0:127]/256;

	print "\n";
	print "パワースペクトル密度のうち始めの 127 点を用いる。残りは\n";
	print "対称的な形となる。さて，パワースペクトル密度をプロットしよう。\n\n";
	print "\tmgplot_title(win, \"Power spectral density\");\n";
	print "\tmgplot_xlabel(win, \"Frequency (Hz)\");\n";
	print "\tmgplot(win, f, Pyy(1:128), {\"power spectral density\"});\n\n";
	pause;

	mgplot_reset(win);
	mgplot_title(win,"Power spectral density");
	mgplot_xlabel(win,"Frequency (Hz)");
	mgplot(win,f,Pyy(1:128),{"power spectral density"});

	print "\n";
	print "ズームインし，200 Hz までをプロットしよう。\n\n";
	print "\tmgplot_title(win,\"Power spectral density\");\n";
	print "\tmgplot_xlabel(win,\"Frequency (Hz)\");\n";
	print "\tmgplot(f(1:50),Pyy(1:50),{\"power spectral density\"});\n\n";
	pause;

	mgplot_reset(win);
	mgplot_title(win,"Power spectral density");
	mgplot_xlabel(win,"Frequency (Hz)");
	mgplot(win,f(1:50),Pyy(1:50),{"power spectral density"});
}

/*  -*- MaTX -*-
 *
 *  【名前】
 *      demo_fourier_expansion() - フーリエ級数展開のデモンストレーション
 *
 *  【形式】
 *      demo_fourier_expansion()
 *
 */

Func void demo_fourier_expansion()
{
	Integer win,k;
	Array t,y,x,yy;

	win = mgplot_cur_win();

	clear;
	print "\n";
	print "矩形波のフーリエ級数展開は奇数次の正弦波の和で表されることを\n";
	print "MaTX を使ってグラフィカルに調べましょう。\n\n";

	print "まず，0 から 10 まで 0.1 きざみで時間ベクトルを作り，全ての点\n";
	print "で sin() を計算する。\n\n";
	print "\tt = [0:0.1:10];\n";
	print "\ty = sin(t);\n";
	print "\tmgplot(win,t,y,{\"sin(t)\"});\n\n";
	print "基本周波数の波形をプロットしましょう。\n\n";
	pause;

	t = [0:0.1:10];
	y = sin(t);
	mgplot(win,t,y,{"sin(t)"});

	print "次に，第 3 周波数の波形を基本周波数の波形に加え，プロットします。";
	print "\n\n";
	print "\ty = sin(t) + sin(3*t)/3;\n";
	print "\tmgplot(win,t,y,{\"sin(t)+sin(3*t)/3\"});\n\n";
	pause;

	y = sin(t) + sin(3*t)/3;
	mgplot(win,t,y,{"sin(t)+sin(3*t)/3"});

	print "さらに，第 1，第 3，第 5，第 7，第 9 周波数の波形を加え，";
	print "プロットします。\n\n";

	print "\ty = sin(t) + sin(3*t)/3 + sin(5*t)/5 + sin(7*t)/7 + sin(9*t)/9;\n";
	print "\tmgplot(win,t,y,{\"\"});\n\n";
	pause;

	y = sin(t) + sin(3*t)/3 + sin(5*t)/5 + sin(7*t)/7 + sin(9*t)/9;
	mgplot(win,t,y,{""});

	print "基本周波数から第 19 周波数まで波形を生成し，それらを\n";
	print "次々に加えた波形を行列に保存します。";
	print "\n\n";

	print "\tt = [0:0.02:3.14];\n";
	print "\ty = Z(10,length(t));\n";
	print "\tx = Z(t);\n\n";

	print "\tfor (k = 1; k <= lenght(t); k = k + (2)) {\n";
	print "\t\tx = x + sin(k*t)/k;\n";
	print "\t\ty((k+1)/2,:) = x;\n";
	print "\t}\n\n";
	pause;

	t = [0:0.02:3.14];
	y = Z(10,length(t));
	x = Z(t);

	for (k = 1; k <= 19; k = k + 2) {
		x = x + sin(k*t)/k; 
		y((k+1)/2,:) = x;
	}

	print "これらの波形をプロットして，矩形波形ができあがる様子を見ま";
	print "しょう。\n";
	print "矩形波の両端が真値に収束しないというギプツ現象に注意しましょう。";
	print "\n\n";

	print "\tmgplot(win,t,y(1:2:9,:),{\"\"});\n";
	print "\tmgplot_title(win,\"The building of a square wave: Gibbs effect\");";
	print "\n\n";
	pause;

	yy = y(1:2:9,:);
	mgplot(win, t, yy, {"1st","5rd","9th","13th","17th"});
	mgplot_title(win, "The building of a square wave: Gibbs effect");
}


/*  -*- MaTX -*-
 *
 *  【名前】
 *      demo_convolution() - 畳み込みのデモンストレーション
 *
 *  【形式】
 *      demo_convolution()
 */

Func void demo_convolution()
{
	Array XY,X,Y,x,y,zfreq,ztime;
	
	clear;
	print "\n";
	print "時間領域における畳み込みは周波数領域における乗算と等価\n";
	print "であるという事実がシステム理論で知られている。MaTX の fft()と\n";
	print "conv() 関数を使って，このことを数値的に調べよう。\n\n";

	print "まず，任意に選ばれた数値からなる2個の数列を作る。\n\n";
	print "\tx = [1 2 3 4];\n";
	print "\ty = [5 4 3 2 1];\n\n";
	pause; 
	x = [1 2 3 4];
	y = [5 4 3 2 1];

	print "時間領域におけるこれらの畳み込みは conv() で簡単に行える。";
	print "\n\n";
	print "\tztime = conv(x,y);\n\n";
	ztime = conv(x,y);
	pause;

	clear;
	print "\n";
	print "周波数領域における乗算を行うため，数列を周波数領域に変換\n";
	print "しなければならない。離散フーリエ変換をするため FFT を使う。\n";
	print "時間領域の畳み込みで計算された z の長さは 8 なので，x と y\n";
	print "の長さが 8 となるよう 0 を追加する。そして，FFT を実行し，\n";
	print "乗算した結果を逆FFTする。\n\n";
	
	print "\tx = [x, Z(1,4)];\n";
	print "\ty = [y, Z(1,3)];\n";
	print "\tX = fft(x);\n";
	print "\tY = fft(y);\n";
	print "\tXY = X * Y;\n";
	print "\tzfreq = Re(ifft(XY));\n\n";
	pause;
	
	x = [x, Z(1,4)];
	y = [y, Z(1,3)];
	X = fft(x);
	Y = fft(y);
	XY = X * Y;
	zfreq = Re(ifft(XY));

	clear;
	print "\n";
	print "FFT の結果と畳み込みの結果を比較する。\n\n";
	print "\tprint ztime;\n";
	print "\tprint round(zfreq);\n\n";
	pause;
	print ztime, "\n", round(zfreq);

	print "\n";
	print "桁落ち誤差の範囲で 2 個の結果が一致する事が分かる。\n";
	print "これを確かめるため周波数領域での結果から丸め誤差による";
	print "誤差を除いた。\n\n";
	pause;

	clear;
	print "\n";
	print "畳み込みの結果と FFT の結果の差を調べよう。\n\n";
	print "\tprint ztime - zfreq;\n\n";
	print ztime - zfreq;
	
	printf("\n\n最大残差は約%16.8E である。\n\n", max(abs(ztime - zfreq)));
}