/*  -*- MaTX -*-
 *
 *  NAME
 *	    sig_demo() - Demonstrations of signal toolbox
 *
 *  SYNOPSIS
 *      sig_demo()
 *
 *  DESCRIPTION
 *      sig_demo() brings up a menu of the available demonstrations.
 *
 *  NOTE
 *      gnuplot is required to display graphs
 *      gnuplot is availabel by ftp from ftp.matx.org.
 *
 */

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

	void demo_fourie_expansion();
	void demo_convolution();
	void demo_spectral_analysis();

	demos = {"Demo of Signal Toolbox",
			 "Spectral analysis using FFTs",
			 "Fourie series expansion",
			 "Convolution and FFT",
			 "E N D"};

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

/*  -*- MaTX -*-
 *
 *  NAME
 *      demo_spectral_analysis() - Demonstration of spectral analysis
 *
 *  SYNOPSIS
 *      demo_spectral_analysis()
 *
 */

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

	win = mgplot_cur_win();

	clear;
	printf("\n");
	printf("Spectral analysis with FFT\n");
	printf("\n");
	pause;

	printf("Make a time axis from t=0 to t=0.25\n\n");
	printf("\tt = [0:0.001:0.25];\n");
	pause;
	t = [0:0.001:0.25];

	printf("\n");
	printf("Make a signal containing 40 Hz and 130 Hz:\n\n");
	printf("\tx = sin(2*PI*40*t) + sin(2*PI*130*t);\n");
	pause;
	x = sin(2*PI*40*t) + sin(2*PI*130*t);

	printf("\n");
	printf("Contaminate the signal with white noise of standard deviation 2\n\n");
	printf("\ty = x + 2*randn(t);\n");
	pause;
	y = x + 2*randn(t);

	printf("\n");
	printf("Plot noisy signal y\n\n");
	printf("\tmgplot_title(win,\"Noisy time domain signal\");\n");
	printf("\tmgplot(win,y(1:50));\n\n");
	pause;

	mgplot_reset(win);
	mgplot_title(win,"Noisy time domain signal");
	mgplot(win,y(1:50));

	printf("Perform fast-Fourier transform (FFT) :\n\n");
	printf("\tY = fft(y,256);\n");
	pause;
	Y = fft(y,256);

	printf("\n");
	printf("Power spectral density:\n\n");
	printf("\tPyy = Re(Y .* conj(Y));\n");
	pause;
	Pyy = Re(Y .* conj(Y));

	printf("\n");
	printf("Make a frequency axis:\n\n");
	printf("\tf = 1000*[0:127]/256;\n");
	pause;
	f = 1000*[0:127]/256;

	printf("\n");
	printf("Plot the power spectral density:\n\n");
	printf("\tmgplot_title(win,\"Power spectral density\");\n");
	printf("\tmgplot_xlabel(win,\"Frequency (Hz)\");\n");
	printf("\tmgplot(win,f,Pyy(1:128));\n");
	pause;

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

	printf("\n");
	printf("Plot only up to 200 Hz:\n\n");
	printf("\tmgplot_title(win,\"Power spectral density\");\n");
	printf("\tmgplot_xlabel(win,\"Frequency (Hz)\");\n");
	printf("\tmgplot(f(1:50),Pyy(1:50));\n");
	pause;

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

/*  -*- MaTX -*-
 *
 *  NAME
 *      demo_fourie_expansion() - Demonstration of Fourie series expansion
 *
 *  SYNOPSIS
 *      demo_fourie_expansion()
 *
 */
Func void demo_fourie_expansion()
{
	Integer win,k;
	Array t,tt,y,x;

	win = mgplot_cur_win();

	clear;
	print "\n";
	print "Make a time vector rangin from 0 to 10 in step 0.1,\n";
	print "and taking the sine of all the elements:\n\n";
	print "\tt = [0:0.1:10];\n";
	print "\ty = sin(t);\n";
	print "\tmgplot(win,t,y,{\"\"});\n\n";
	print "Plot the fundamental.\n\n";
	pause;

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

	print "Add the 3rd harmonic to the fundamental\n\n";
	print "\ty = sin(t) + sin(3*t)/3;\n";
	print "\tmgplot(win,t,y,{\"\"});\n\n";
	pause;

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

	print "Use the 1st, 3rd, 5th, 7th, and 9th harmonics:\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 "Go from the fundamental to the 19th harmonic\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 "Plot these waves successively\n\n";

	print "\tmgplot(win,y(1:2:9,:),{\"\"});\n";
	print "\tmgplot_title(win,\"Gibbs' effect\");";
	print "\n\n";
	pause;

	tt = ONE(5,1)*Matrix(t);

	mgplot(win, tt, y(1:2:9,:), {"1st","5rd","9th","13th","17th"});
	mgplot_title(win, "Gibbs effect");

	mgplot_reset(win);
}


/*  -*- MaTX -*-
 *
 *  NAME
 *      demo_convolution() - Demonstration of convolution
 *
 *  SYNOPSIS
 *      demo_convolution()
 */

Func void demo_convolution()
{
	Array XY,X,Y,x,y,zfreq,ztime;
	
	clear;
	print "\n";
	
	print "Make two sequences of numbers:\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 "Convolve them in the time domain\n";
	print "\tztime = conv(x,y);\n\n";
	ztime = conv(x,y);
	pause;

	print "Perform FFT to the two sequences\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));

	print "Comparie the FFT result with the convolution result:\n\n";
	print "\tprint ztime;\n";
	print "\tprint round(zfreq);\n\n";
	print ztime, "\n", round(zfreq);

	print "\n";
	print "The two are the same.\n";
	print "\n";
	pause;
}