/***************************************************************/
/*Copyright (c) by Continocean Tech Inc. 	       	       */
/*12 Mountain Ave., Montville, New Jersey, USA 07045           */
/*Tel: 201-257-1912 Fax: 201-257-9634			       */
/*Email: continocean@webexpert.net			       */
/*All rights reserved. No part of this work covered by the     */
/*copyright hereon may be reproduced or used in any form or by */
/*any means----graphic, electronic, or mechanical, including   */
/*photocopying, recording, taping, or information storage and  */
/*retrieval systems----without written permission of 	       */
/*Continocean Tech inc.					       */
/***************************************************************/

                     ENGSIM for MS VC++

                  Engineering Simulation Library
                    for Microsoft Visual C++

                   Usage of ENGSIM Instructions


IMPORTANT MESSAGE:
    Please use the Microsoft Development Studio's File-OpenWorkspace to
open \engsim\sim project space. Then use File-Open to open file main.cpp
for editing. There are several examples under \engsim\examples which are 
included in main.cpp. You can edit main.cpp and the examples for exercises.


PART I. COMPLEX, MATRIX, COMPLEX MATRIX, INTEGER MATRIX AND STRING MATRIX

Complex
	    declare a complex number
e.g.        Complex a;

matrix	  		
	    declare a matrix
	    matrix a(2,2), h(10,20), p(200);
	    //a(1,1)...... a(2,2)
	    //p(1) ...... p(200)

matrixn			
	    declare a matrix with negative indexing
	    matrixn a(2,2); //a(-2,-2) ...... a(2,2)
	    matrixn p(100); //p(-100) ......p(100)

cmatrix	  	
	    declare a complex matrix
	    cmatrix a(2,2),b(2,2),p(100);

cmatrixn	
	    declare a complex matrix with negative indexing
	    cmatrixn a(2,2),b(2,2),p(100);
	    //a(-2,-2)......a(2,2), p(-100)......p(100)

imatrix	   	
	    declare an integer matrix
	    imatrix a(2,2), b(2,2), p(100);
	    //a(1,1) ...... a(2,2)
	    //p(1) ...... p(100)

imatrixn	
	    declare an integer matrix with negative indexing
	    imatrixn a(2,2),b(2,2),p(100);
	    //a(-2,-2) ... a(2,2)
	    //p(-100) ... p(100)

smatrix	  	
	    declare a string matrix
	    smatrix a(2,2),p(100);
	    //a(1,1) ... a(2,2)
	    //p(1) ... p(100)

=           
	    equal operator for complex, matrix, matrixn, 
	    cmatrix, cmatrixn, imatrix, imatrixn
            b= 1.0+2.0*J;
	    b=a;

+           
	    plus operator for complex, matrix, matrixn, 
	    cmatrix, cmatrixn, imatrix, imatrixn
	    c= a+b;

-           
	    minus operator for complex, matrix, matrixn, 
	    cmatrix, cmatrixn, imatrix, imatrixn
            c= a-b;
            c=-b;

*           multiplication operator for complex, matrix, matrixn, 
	    cmatrix, cmatrixn, imatrix, imatrixn
            c= a*b;

/	    
	    division operator for complex, matrix, matrixn, 
	    cmatrix, cmatrixn
	    c=a/b;

~           
	    conjugate operator for complex, cmatrix, cmatrixn
            c= (~a);

==		    
	    logic equal operator for complex, matrix, matrixn, 
	    cmatrix, cmatrixn, imatrix, imatrixn, smatrix
	    if(a==b){...};

!=		    
	    logic not equal operator for complex, matrix, matrixn, 
	    cmatrix, cmatrixn, imatrix, imatrixn, smatrix
	    while(a!=b){...};

|			
	    or for imatrix, imatrixn
	    c=a|b;

&			
	    and for imatrix, imatrixn
	    c=a&b;

^	    
	    X-or for imatrix, imatrixn
	    c=a^b;

( , )		
	    an element of a matrix, matrixn, cmatrix, cmatrixn
	    imatrix, imatrixn, smatrix
	    a(2,-5)=10.1+2.1*J;
	    p(-5)=3.0+9.9*J;  

J           
	    pure unit imaginary for Complex
            a= 1.0+2.0*J;

abs()		
	    absolute for complex
	    take absolute of each element for matrix, matrixn, cmatrix, cmatrixn
	    Complex ca;
	    double dc;
	    dc=abs(ca);
	    cmatrix a(3,3);
	    matrix b(3,3);
	    b=abs(a);

amplitude()	amplitude of cmatrix and cmatrixn
		cmatrix a(10);
		matrix z(10);
		z=amplitude(a);
		cmatrixn a(10);
		matrixn z(10);
		z=amplitude(a);

arg()		same as phase()

 
corrm(a, b, n)	
			n-by-n correlation matrix of a and b for cmatrix, cmatrixn
			c=corrm(a,b,2);



det()		determinant of a matrix, matrixn, cmatrix, cmatrixn
			imatrix, imatrixn
			double dc;
			dc=det(a);

eigen(a,d,z)
	compute eigenvalues d(1) ... d(n) and eigen vectors 
	  z(1,1)          z(1,2)
	 (......)=z1     (......)=z2 ......
	  z(n,1)          z(n,2)
	a*z1=d(1)*z1     a*z2=d(2)*z2 ......
	matrix a(2,2);
	matrix d(2),z(2,2);	//if a is symmetric
	eigen(a,d,z);
	cmatrix d(2),z(2,2);	//if a is nonsymmetric
	eigen(a,d,z);

exp()		exponential for complex
		    exp(J*2.0*PI);


float2string(fa, sa, len) 
			convert float number fa to string sa(1) ... sa(len)


fprintf(,)	print to a disk file for complex, matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn, smatrix,
		    double, integer, character
		    FILE *fp;
		    fp=fopen("file.out","w");
		    fprintf(fp, c);
		    fclose(fp);
		    fprintf("file.out",c); //append only

fscanf( , )	read from a disk file for complex, matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn, smatrix,
		    double, integer, character
		    FILE *fp;
		    fp=fopen("file.in","r");
		    fscanf(fp, a);
		    fclose(fp);

HH			hermitian for cmatrix, cmatrixn
			c=HH(a);

identity	identity matrix (diagonal elements are all 1's, others all 0's)
		    (for matrix only)
		    a=identity(10);	   //10x10 identity matrix


identityn	identity matrix with negative indexing
		    (for matrix with negative indexing only)
		    a=identityn(10);	   //21x21 identity matrixn


init(, ...)	initialize a matrix, matrixn, cmatrix, cmatrixn
			imatrix, imatrixn, smatrix
			matrix a(2,2);
			init(a, 1.0,2.0,
				4.0,5.0
				);
			matrixn a(2,2);
			a=init(a,1.0,2.0,3.0,4.0,5.0,
				 2.0,2.0,3.0,4.0,5.0,
				 3.0,2.0,3.0,4.0,5.0,
				 4.0,2.0,3.0,4.0,5.0,
				 5.0,2.0,3.0,4.0,5.0
				 );
			cmatrix a(2,2);
			a=init(a,
				1.0+2.0*J, 2.0+3.0*J,
				1.1+2.1*J, 2.1+3.1*J
				);
			cmatrixn a(2,2);
			a=init(a,
				1.0+2.0*J, 2.0+3.0*J, 1.0+2.0*J, 2.0+3.0*J, 1.0+2.0*J,
				1.1+2.1*J, 2.1+3.1*J, 1.0+2.0*J, 2.0+3.0*J, 1.0+2.0*J,
				1.0+2.0*J, 2.0+3.0*J, 1.0+2.0*J, 2.0+3.0*J, 1.0+2.0*J,
				1.1+2.1*J, 2.1+3.1*J, 1.0+2.0*J, 2.0+3.0*J, 1.0+2.0*J,
				1.0+2.0*J, 2.0+3.0*J, 1.0+2.0*J, 2.0+3.0*J, 1.0+2.0*J
				);
			imatrix a(2,2);
			init(a,
		   		1,2,
				4,5
				);
			imatrixn a(2,2);
			init(a,
		   		1,2,3,2,3,
				4,5,6,2,3,
				7,8,9,2,3,
				4,5,6,2,3,
				7,8,9,2,3
				);
			smatrix a(2,2);
			a=init(a,
				'H','i',
				'O','K');

imag()      imaginary of a complex number
            double da;
            da=imag(a);
		    imaginary of Complex matrix
			matrix mc(2,2);
			mc=imag(a);
			matrixn mc(2,2);
			mc=imag(a);




int2string(ia, sa, len)	
			convert integer ia to string sa(1) ... sa(len)

joinr(a, b);
	right join matrices a,b as (a b) for matrix, cmatrix, imatrix
	
joind(a, b);
	down join matrices a,b as 
	   a(1,1)......
	   ............ 
	(  b(1,1)......  )
	   ............
	for matrix, cmatrix, imatrix

log2(x)		base 2 log, x double

log(x,n)	base n log, n=integer, x double

log(x)		
		base e log, x Complex

Max()		find maximum of all elements for matrix
			double a_max=Max(a);

Min()		find minimum of all elements for matrix
			double a_min=Min(a);

Max( , , ,)	find maximum of all elements for matrix
			Max(a, y_max, i_max, j_max); //y_max=a(i_max,j_max)

Min( , , ,)	find minimum of all elements for matrix
			Min(a, y_min, i_min, j_min); //y_min=a(i_min,j_min) 

mean(a)		mean of matrix, cmatrix
		double x;
		matrix a;
		x=mean(a);
		Complex y;
		cmatrix b;
		y=mean(b);

normsq()		real*real+imag*imag for complex	
		    sum a(i,j)*a(i,j) for matrix, matrixn, imatrix, imatrixn
		    sum a(i,j)*(~a(i,j)) for cmatrix, cmatrin
			double dc;
			dc=norm(a);

phase()		phase of a complex number, cmatrix and cmatrixn
		Complex c;
		double x;	
		x=phase(c);
		cmatrix a(10);
		matrix z(10);
		z=phase(a);
		cmatrixn a(10);
		matrixn z(10);
		z=phase(a);


pow( , )	power for complex number
		    Complex c;
		    pow(c,2.0);
		    pow(c,3);
		    pow(c,1.0+J*2.0);


powi(x,n)	x^n
		n=integer

printw()	print to screen for complex, matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn, smatrix,
		    double, integer, character
		    printw(a);
		    printw(2.0+3.0*J);
		    printw(1.0);
		    printw(1);
		    printw('c');
		    printw("test=%g, dog=%d\n",0.5, 100);



real()      real of a complex number
		    double da;
            da=real(a);
			real of Complex matrix
			matrix mc(2,2);
			mc=real(a);
			matrixn mc(2,2);
			mc=real(a);

shiftl( )	shift left by 1 for matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn

shiftr( )	shift right by 1 for matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn

shiftu()	shift up by 1 for matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn

shiftd()	shift down by 1 for matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn

sub(a, m1, n1, m2, n2)	
			Take submatrix beginning at a(m1,n1) and ending at a(m2,n2)
			for matrix, matrixn, cmatrix, imatrix		
			(matrixn becomes a matrix)
			sub(a,m1,n1,m2,n2);

sum(a)		sum for matrix, cmatrix
		double x;
		matrix a;
		x=sum(a);
		Complex y;
		cmatrix b;
		y=sum(b);
		

TT()		transpose for matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn, smatrix
			c=TT(a);



unitary		unitary matrix (all elements are 1's)
		    (for matrix)
		    a=unitary(10);	   //10x10 unitary matrix
		    a=unitary(5,10);	   //5x10 unitary matrix

unitaryn	unitary matrix with negative indexing
		    (for matrix with negative indexing)
		    a=unitaryn(10);	   //21x21 unitary matrixn

var(a)		variance for matrix, cmatrix
		double x;
		matrix a;
		x=var(a);
		cmatrix b;
		x=var(b);


PART II. GRAPHICS (2D and 3D)

In 2D plotting, if y(k)>=1.0e40, the point k is ignored. The feature 
is useful for plotting y1,y2,... with different x(k). Simply combining
x1,x2,... into x and let unwanted y1(k),y1(k),...=1.0e40. The dimension 
of x and y may need to be adjusted.

In 3D plotting, if y(k)>=1.0e40, the point is interrupted and the curve 
is restarted from the next point. The feature is useful for plotting 3D
objects to avoid unwanted lines. Simply insert z(k)=1.0e40 to delete a line
between z(k-1) and z(k+1).

plot(y)		plot array y 
			matrix y(100);
			plot(y);

plot(x, y)	plot array, x=X-axis, y=Y-axis
			matrix x(100), y(110);
			plot(x, y);

plot(x, y1, y2)		
			plot array, x=X-axis, y1,y2,...=Y-axis
			matrix x(100), y1(110), y2(100);
			plot(x, y1, y2);

plot(x, y1, y2, y3)
plot(x, y1, y2, y3, y4);
plot(x, y1, y2, y3, y4, y5);

plot(y, "legend");
plot(x, y, "legend");
plot(x, y1, "legend1", y2, "legend2");
plot(x, y1, "legend1", y2, "legend2", y3, "legend3");
plot(x, y1, "legend1", y2, "legend2", y3, "legend3", y4, "legend4");
plot(x, y1, "legend1", y2, "legend2", y3, "legend3", y4, "legend4", y5, "legend5");


plot3(x, y, z)
			plot 3-D curve using data array, x=X-axis, y=Y-axis, z=Z-axis
			matrix x(100), y(110), z(110);
			plot3(x, y, z);

plot3(x1, y1, z1, x2, y2, z2)		
			plot 3-D curve using data array, x=X-axis, y=Y-axis, z=Z-axis
			matrix x1(100), y1(110), z1(110), x2(100), y2(100), z2(110);
			plot3(x1, y1, z1, x2, y2, z2);

plot3(x, y, z, "legend");
plot3(x1, y1, z1, "legend1", x2, y2, z2, "legend2");


In the above plot(), all x,y,z are matrix (array).

In the following plot(), x is cmatrix.
plot(cmatrix x);
plot(cmatrix x, "legend");
plot(cmatrix x1, cmatrix x2);
plot(cmatrix x1, "legend1", cmatrix x2, "legend2");


printt(x, y, " ", number1, ...)
	print any text or number to a figure with position (x,y)
	x,y are type-int. -10000<x,y<10000 
	printt(2000,3000,"Sweet Energy=%g",eng)
	printt(2000,3000,"Figure 1.")

delete_axis()
	delete axises and legends on the current figure.

set_axis_scale("","")
	for 2D. Set axis scale to logarithmic scale or regular scale.
	set_axis_scale("","log"); //set x-axis to regular and y-axis to logarithmic scale
	set_axis_scale("log","log"); //set x-axis and y-axis to logarithmic scale
	set_axis_scale("",""); //set x-axis and y-axis to regular scale (default)
	


set_rotation_angle(x, y, z);
	for 3D. x,y,z are double-type. 
	If >0, counter-clockwise rotation.
	set_rotation_angle(-PI/24,PI/12,-PI/6); //default
	set_rotation_angle(PI/24,PI/12,-PI/6); 

set_focal_distance(x)
	for 3D.
	default is 100 times the difference between the maximum and minimum X-coordinates.
	set_focal_distance(101.0);

set_picture_size(x);
	scale picture
	x=1.0 default size
	0.0<x<1.0  shrink by x
	1.0<x<2.0  expand by x
	set_picture_size(1.0); //default
	set_picture_size(1.5);
	set_picture_size(0.5);


set_picture_style("style");
	set style to curve (default), bar, scatter
	set_picture_style("curve");
	set_picture_style("bar");
	set_picture_style("scatter");
	for 2-D only

set_picture_position(x, y);
	relocate picture position
	-10000<x,y<10000
	set_picture_position(0, 0);	//default
	set_picture_position(2000, 1000); //move right by 2000, up by 1000
	set_picture_position(-2000, -1000); //move left by 2000, down by 1000

set_axis_title("x_title", "y_title");

set_axis_title("x_title", "y_title", "z_title");

set_text_color("color");
	color=black (default), blue, green, cyan, red, magenta, yellow, white 
	      darkblue, darkgreen, darkcyan, darkred, darkmagenta, darkyellow
	       darkgrey, lightgrey

set_text_position(x, y);	
	x=row, y=column
	x,y>0
	set_text_position(100, 20);
	

set_text_size(x);
	set text font size
	set_text_size(10); //default
	set_text_size(15);




PART III. COMMUNICATIONS ENGINEERING AND SIGNAL PROCESSING SIMULATION


betai(a, b, x)
	Incomplete Beta function I_x(a,b)
	0<x<1, a,b>0
	double x=0.1,a1.0,b=0.1,y;
	y=betai(a, b, x);
	


convol(x, y)
	convolution of arrays x and y for matrix, matrixn, cmatrix and cmatrixn
	matrix x(100),y(100),z(200);
	z=convol(x,y);
	matrixn xn(100),yn(100),zn(200);
	zn=convol(xn,yn);
	cmatrix cx(100),cy(100),cz(200);
	cz=convol(cx,cy);
	cmatrixn cxn(100),cyn(100),czn(200);
	czn=convol(cxn,cyn);

corr(x, y)
	correlation of arrays x and y for matrix (array) and cmatrix (array)
	matrix x(100), y(100);
	matrixn z(100);	//z(-100) ... z(100)
	z=corr(x,y);
	cmatrix x(100), y(100);
	cmatrixn z(100);	//z(-100) ... z(100)
	z=corr(x,y);

corrnorm(x, y)
	normalized correlation of arrays x and y for matrix (array) and cmatrix (array)
	matrix x(100), y(100);
	matrixn z(100);	//z(-100) ... z(100)
	z=corrnorm(x,y);
	cmatrix x(100), y(100);
	cmatrixn z(100);	//z(-100) ... z(100)
	z=corrnorm(x,y);


crandom();
	generate uniform-distributed integer in (0, 2147483647)
	int a;
	a=crandom();

czt(x, m, w, a)
	Chirp-Z Transform
	Complex a=5.0*exp(J*PI/4); //starting point
	Complex w=0.999*exp(J*PI/6/50); //spiral rate
	int m=200;
	matrix x(100);
	cmatrix X(m);
	X=czt(x,m,w,a);


fft(matrix &h)
	compute spectrum of an array h using FFT
	sampling-rate=1
	n=number of taps of h
	if n<pow(2,r), 0's are appended to make =2^r
	h(1)...h(n)=filter taps
	H(1)...H(n)=complex spectrum
	int n=4098;	//2^12
	matrix h(n),freq(n);
	cmatrix H(n);
	H=fft(h);
	matrix sp_amp(n), sp_phase(n);
	sp_amp=amplitude(H);
	sp_phase=phase(H);
	for(i=1;i<=n;i++) H_freq(i)=i-1;
	plot(freq, sp_amp, "Amplitude", sp_phase, "Phase");


fft(matrix &h, matrix &sp_freq, cmatrix & H, int n, float Ts):
	compute spectrum of an array h using FFT
	n=pow(2,r)=number of points in FFT, 
	h(1)...h(n)=filter taps
	H(1)...H(n)=complex spectrum
	Ts=sampling-period (default=1.0)
	int n=4098;	//2^12
	matrix h(n),freq(n);
	cmatrix H(n);
	float Ts=0.01;
	fft(h, freq, H, n, Ts);
	matrix sp_amp(n), sp_phase(n);
	sp_amp=amplitude(H);
	sp_phase=phase(H);
	plot(freq, sp_amp, "Amplitude", sp_phase, "Phase");

fft(x, sp_freq, sp_amp, sp_phase, n, Ts);
	compute spectrum of an array x using FFT
	n=pow(2,r)=number of points in FFT, 
	dimension of sp_freq (frequency), sp_amp (amplitude) and sp_phase (phase) is n/2
	Ts=sampling-period
	int n=4098;	//2^12
	matrix x(n), freq(n/2), sp_amp(n/2), sp_phase(n/2);
	float Ts=0.01;
	fft(x, freq, sp_amp, sp_phase, n, Ts);
	plot(freq, sp_amp, "Amplitude", sp_phase, "Phase");

f_test(a1, a2)
	Compute F test for testing the
	difference of variances for  two arrays 
	a1(1)...a1(n1), a2(1)...a2(n2) 
	A very small returned value (<0.05) implies 
	that the variances are significantly different
	matrix a1(100),a2(100);
	double y;
	y=f_test(a1, a2);

f_stat(a1, a2)
	Compute F statistic for testing the
	difference of variances for  two arrays 
	a1(1)...a1(n1), a2(1)...a2(n2) 
	matrix a1(100),a2(100);
	double y;
	y=f_stat(a1, a2);

gamma(xx)
	Compute gamma function
	double x=1.5,y;
	y=gamma(x);

gammln(xx)
	Compute ln(gamma) function
	double x=1.5,y;
	y=gammln(x);

ifft(cmatrix & H, matrix &h_time, matrix &h, int n, float fs)
	Inverse FFT
	n=pow(2,r)=number of points in IFFT, 
	h(1)...h(n)=filter taps
	H(1)...H(n)=complex spectrum
	fs=unit frequency (default=1.0)
	h_time(1)...h_time(n)=k*Ts  (Ts=1/fs/n=sampling period)
	int n=4098;	//2^12
	matrix h(n),h_time(n);
	cmatrix H(n);
	float fs=0.01;
	ifft(H, h_time, h, n, fs);

	

integral(func, a, b, step)	    
        integrate function func from a to b with step-size step
	double function(double x){
	       double y;
	       ......
	       return y;
	}
	double a=-1.0,b=2.0,step=0.01;
	integral(func,a,b,step);


lpf(fmT, n, h);
	generate low pass filter coefficient h(1) ...... h(n) with cut-off frequency fmT
	int n=20;
	float fmT=1000;
	matrix h(n);
	lpf(fmT, n, h);

median(a)
	Compute the median of an array a(1), ..., a(n)
	double x;
	matrix a(100);
	x=median(a);
	

min(func, dfunc, n, xm, fm, mu)
	find minimum point xm and minimum value fm for function func 
	(first order derivative dfunc), n=number of variables
	mu=stepsize (using steepest descending algorithm)
	double function(matrix &x){...}
	matrix dfunction(matrix &x){...}
	int n=2;
	matrix xm(n);
	double fm,mu=0.1;
	min(function, dfunction, n, xm, fm, mu);

min(func, dfunc, xi, n, xm, fm, mu)
	same as the above except that xi=initial point
	......
	xi=init(xi,1.0,1.0);
	min(func, dfunc, xi, n, xm, fm, mu);

norm(m, v)
	generate Gaussian (Normal) distributed random variable with mean m 
	and variance v*v
	float a;
	a=norm(1.0, 2.0); //Gaussian distributed random variable with mean 1.0 and variance 4.0

pcm(x,step,n)
	n-bit uniform pulse-code-modulation (PCM) of x with step size step
	......
	x=b10...01 if step/2       <= x < step/2+step
	x=b10...00 if 0            <= x < step/2
	x=b01...11 if -step/2      <= x < 0
	x=b01...10 if -step/2-step <= x < -step/2
	......
	double x=10.0, step=0.625;  //step=10/2^4
	int n=5;
	binary b(n);
	b=pcm(x,step,n);
	printw("%g=",x);
	printw(b);

pcmd(b)
	decode uniform PCM
	double x;
	int n=5;
	binary b(n);
	......
	x=pcmd(b);
	printw("%g=",x);
	printw(b);


pcm(x, step, x_max, mu, n)
	n-bit mu-law pulse-code-modulation (PCM) of x with step size step
	double x=5.0, step=0.0625, x_max=10.0, mu=255; //step=1.0/2^4
	int n=5;
	binary b(n);
	b=pcm(x,step,x_max,mu,n);
	printw("%g=",x);
	printw(b);

pcmd(b,mu)
	decode mu-law PCM
	double x, mu=255;
	int n=5;
	binary b(n);
	......
	x=pcmd(b,mu);
	printw("%g=",x);
	printw(b);



pdf(a, x, y, n=10)
	Given an array of data a(1), ..., a(m), 
	estimate the probability density function
	of the data a(i). a(i) is divided into n (default=10)
	steps from Max(a) to Min(a) and returned as x(1), ..., x(n)
	the probability density is returned as y(1), ..., y(n)
	matrix a(100000),x(10),y(10);
	pdf(a,x,y);	//default n=10
	pdf(a,x,y,10);	
	plot(x,y,"Estimated PDF");

 
phasecontinue(p)
	make the phase array p to be continuous by subtracting or adding 2PI.
	matrix p(100);
	phasecontinue(p);
	plot(p);	//the plot  will be continuous

phasefollow(p1, p2)
	make phase array p1 follow phase array p2 by subtracting or adding 2PI.
	matrix p1(100),p2(100),x(100);
	phasefollow(p1,p2);
	......
	plot(x,p1,p2);	//

qpiqpskmod(a2, a1)
	PI/4-QPSK modulation. 
	binary sequence (a2 a1)to QPSK mapping: 00=PI/4, 01=3*PI/4, 11=5PI/4, 10=7PI/4.
	cmatrix x(100);
	imatrix a(200);
	a(2)=0;
	a(1)=1;
	x(1)=qpiqpskmod(a(2), a(1));

qpiqpskdem(x, a2, a1);
	PI/4-QPSK demodulation. 
	binary sequence (a2 a1)to QPSK mapping: 00=PI/4, 01=3*PI/4, 11=5PI/4, 10=7PI/4.
	cmatrix x(100);
	imatrix a(200);
	x(1)=exp(J*PI/4.0);
	qpiqpskdem(x(1), a(2), a(1));

qpidqpskmod(a2, a1, x1);
	PI/4-DQPSK modulation. 
	binary sequence (a2 a1)to QPSK mapping: 00=PI/4, 01=3*PI/4, 11=5PI/4, 10=7PI/4.
	cmatrixn x(100);
	imatrixn a(200);
	a(0)=0;
	a(-1)=1;
	x(-1)=exp(J*PI/4.0);
	x(0)=qpidqpskmod(a(0), a(-1),x(-1));

qpidqpskdem(y0, y1, a2, a1);
	PI/4-DQPSK demodulation. 
	binary sequence (a2 a1)to QPSK mapping: 00=PI/4, 01=3*PI/4, 11=5PI/4, 10=7PI/4.
	cmatrix y(100);
	imatrix a(200);
	y(-1)=exp(J*PI/4.0);
	y(0)=exp(J*3.0*PI/4.0);
	qpidqpskdem(y(-1), y(0), a(0), a(-1));


ray(m, v)
	generate Rayleigh (m=0) or Rician (m!=0) distributed variable
	float a;
	a=ray(1.0, 2.0);

ray_fade(fm, T, k);
	generate Jakes Rayleigh fading channel coefficient
	Doppler frequency is fm Hz, T is symbol period, 
	k (>0) is the kth time.
	(setting k=1 and k=1000 respectively gives two-ray channel or diversity channel)
	cmatrixn c1(100), c2(100);
	float fm=75.0, T=41.0e-6;
	c1(0)=ray_fade(fm, T, 1);
	c2(0)=ray_fade(fm, T, 1000); //c1 and c2 are two independent rays.

rcos_sr(float b, float T, float step, int k):
	Square-root raised-cosine pulse-shaping function.
	b=rolloff parameter
	spectrum is
              T                                       0         <= |f| <= (1-b)/T/2
	X(f)= T cos{(PI*T/b/2)*[|f|-(1-b)/T/2]}     (1-b)/T/2 <= |f| <= (1+b)/T/2
              0                                       otherwise
	step*k=time point
	matrixn h(16);
	float b=0.3, T=40.0e-6, step=T/4;  //each symbol has 4 samples, 4 times oversampling
	for(int k=-16;k<=16;k++) h(k)=rcos_sr(b, T, step, k);
	Ref: J. Proakis, Digital Communications, 2ed edition, McGraw-Hill, 1989.

rcos_ps(matrixn &h, int n, float b, float T, int m):
rcos_ps(matrix &h, int n, float b, float T, int m):
	compute filter taps of square-root raised-cosine pulse-shaping
	matrixn h: h(-n),..., h(0),...h(n)=filter taps, 2n+1=number of taps
	matrix  h: h(1),...h(n)=filter taps, n=number of taps
	b=rolloff parameter
	T=symbol period
	T/m=sampling-period, default m=4
	matrixn h(16);
	float b=0.3, T=40.0e-6;  
	//each symbol has 4 samples, 4 times oversampling
	rcos_ps(h, 16, b, T): //default m=4, matrixn h has 2*16+1 taps
	int m=2; 
	matrix g(17);
	rcos_ps(g, 17, b, T, m): //m=2, matrix g has 2*8+1=17 taps
	
	Ref: J. Proakis, Digital Communications, 2ed edition, McGraw-Hill, 1989.
	
	
select(k, a)
	Select the kth smallest number in an array a(1),...,a(n).
	Returned array is reordered in such a way that the kth element is the kth smallest
	with all the smaller elements moved to the left-hand side of it and all the larger elements 
	moved to the right-hand side of it.
	int n=100;
	matrix a(n),x(n);
	x=select(10,a);

sign(x)
	return the sign (1.0 or -1.0) of x
	double x,xs;
	xs=sign(x);

sinc(x);
	sinc function sin(x)/x
	double y;
	y=sinc(1.5*PI);

sa(x);
	same as sinc function sin(x)/x

spectrum(x, sp_freq, sp_amp, sp_phase, n, Ts);
	compute spectrum of an array x using FFT
	n=pow(2,r)=number of points in FFT, 
	dimension of sp_freq (frequency), sp_amp (amplitude) and sp_phase (phase) is n/2
	Ts=sampling-period
	int n=4098;	//2^12
	matrix x(n), sp_freq(n/2), sp_amp(n/2), sp_phase(n/2);
	float Ts=0.01;
	fft(x, sp_freq, sp_amp, sp_phase, n, Ts);
	plot(freq, sp_amp, "Amplitude", sp_phase, "Phase");


spectrum(func, sp_freq, sp_amp, sp_phase, n, T);
	compute spectrum of function func
	T=function period, from 0 to T
	n=pow(2,r)=number of points in FFT, 
	dimension of sp_freq (frequency), sp_amp (amplitude) and sp_phase (phase) 
	is n/2
	sampling-period=T/n
	float function(float x)
	{...... }
	int n=4098;	//2^12
	matrix sp_freq(n/2), sp_amp(n/2), sp_phase(n/2);
	float T=100.0;
	spectrum(function, sp_freq, sp_amp, sp_phase, n, T);
	plot(freq, sp_amp, "Amplitude", sp_phase, "Phase");

tail(func, a)
	integrate the tail of the function func from a to +infinity
	double function(double x){
	       double y;
	       ......
	       return y;
	}
	double a=1.0;
	tail(func,a);


t_stat(a1, a2)
	Compute Student-t statistic for testing the
	difference of means for two arrays 
	a1(1)...a1(n1), a2(1)...a2(n2) 
	with the SAME variance
	matrix a1(100),a2(100);
	double y;
	y=t_stat(a1, a2);

t_test(data1, data2);
	Compute Student-t test for testing the
	difference of mean for  two arrays 
	a1(1)...a1(n1), a2(1)...a2(n2) 
	with the SAME variance.
	A very small returned value (<0.05) implies 
	that the means are significantly different
	matrix a1(100),a2(100);
	double y;
	y=t_test(a1, a2);



t_statv(a1, a2)
	Compute Student-t statistic for testing the
	difference of means for two arrays 
	a1(1)...a1(n1), a2(1)...a2(n2) 
	with DIFFERENT variance
	matrix a1(100),a2(100);
	double y;
	y=t_statv(a1, a2);

t_testv(data1, data2);
	Compute Student-t test for testing the
	difference of mean for  two arrays 
	a1(1)...a1(n1), a2(1)...a2(n2) 
	with DIFFERENT variances.
	A very small returned value (<0.05) implies 
	that the means are significantly different
	matrix a1(100),a2(100);
	double y;
	y=t_testv(a1, a2);
	
ztransform(matrix &h,matrix &f,matrix &a,matrix &p,int m,int n,float T):
	compute spectrum of an FIR filter using Z-transform
	h(m)=taps of an FIR filter=h(1)+h(2)z^{-1}+...+h(m)z^{-(m-1)}
	f(n)=array of frequencies
	a(n)=array of amplitudes of H(z)=H(e^(2*PI*f*T))
	p(n)=array of phases of H(z)=H(e^(2*PI*f*T))
	m=number of taps of h
	n=number of discrete frequencies
	T=sampling period of h, 1/T=freq-domain-period
	int n=4096,m=33;	
	matrix h(n), sp_amp(n/2), sp_freq(n/2), sp_phase(n/2);
	float T=float(41.0e-6), b=(float)0.3;
	rcos_ps(h, 33, b, T);	//square-root raised-cosine pulse-shaping
	ztransform(h, sp_freq, sp_amp, sp_phase, m, n, Ts);
	set_axis_title("Frequency (Hz)"," ");
	set_caption("Spectrum of an FIR filter using Z transform");
	plot(sp_freq, sp_amp, "Amp");