From 211544b0b8dd7e6a396c4a99167bf2cca9f3fcfc Mon Sep 17 00:00:00 2001
From: "Robin. P. Clark" <robin@energytechnologycontrol.com>
Date: Mon, 16 Sep 2019 13:01:05 +0100
Subject: [PATCH] Z plane analysis in C use

make fa
---
 Makefile        |   7 ++
 fa.gpt          |   8 +++
 freq_analysis.c | 171 ++++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 186 insertions(+)
 create mode 100644 fa.gpt
 create mode 100644 freq_analysis.c

diff --git a/Makefile b/Makefile
index ab4a4cf..7393a4b 100644
--- a/Makefile
+++ b/Makefile
@@ -5,3 +5,10 @@ two_pole:
 	./tp > tp.dat
 	vi tp.dat
 	gnuplot < tp.gpt
+
+
+fa: freq_analysis.c 
+	gcc freq_analysis.c -o fa -lm
+	./fa > fa.dat
+	vi fa.dat
+	gnuplot < fa.gpt
diff --git a/fa.gpt b/fa.gpt
new file mode 100644
index 0000000..590bdc9
--- /dev/null
+++ b/fa.gpt
@@ -0,0 +1,8 @@
+
+
+
+plot "fa.dat" using 1:8
+
+!sleep 10
+!sleep 10
+!sleep 10
diff --git a/freq_analysis.c b/freq_analysis.c
new file mode 100644
index 0000000..9210845
--- /dev/null
+++ b/freq_analysis.c
@@ -0,0 +1,171 @@
+#include <stdio.h>
+#include <math.h>
+
+double xz[10]; //  = {0.0}; // Z transform numerator or zeros Z^0 + Z^1 + Z^2 etc. 
+double yp[10]; // = {0.0}; // Z transform denominator or poles Z^0 + Z^1 + Z^2 etc.
+
+
+typedef struct complex_number {
+	double real;
+	double imaginary;
+} COMPLEX;
+
+#define NYQUIST 1000
+#define PI 3.14159265358979323844 
+
+double arg ( COMPLEX x) {
+	return atan2(x.imaginary, x.real);
+}
+
+double mag ( COMPLEX x ) {
+	return sqrt ( x.real * x.real + x.imaginary * x.imaginary );
+}
+
+// add x1 and x2 result in x1
+void complex_add ( COMPLEX * x1, COMPLEX  x2 ) {
+	x1->real = x1->real + x2.real;
+	x1->imaginary = x1->imaginary + x2.imaginary;
+}
+
+void complex_print ( COMPLEX x ) {
+	printf("x.r = %f x.i = %f ",x.real,x.imaginary);
+}
+
+void complex_scalar_mul ( double m, COMPLEX * x) {
+	x->real = m * x->real;
+	x->imaginary = m * x->imaginary;
+}
+
+void complex_pow ( double p, COMPLEX * x ) {
+	double r  = arg(*x) + PI * 2.0;
+	double m  = mag(*x);
+	
+	r *= p;
+	while (r>360.0) r -= 360.0; // r %= ( PI * 2.0);
+	m = pow(m,p);
+
+	x->real = cos(r) * m;
+        x->imaginary = sin(r) * m;	
+}
+
+COMPLEX complex_mul ( COMPLEX a, COMPLEX b ) {
+
+	COMPLEX ans;
+
+	ans.real = a.real * b.real - a.imaginary * b.imaginary;
+	ans.imaginary = a.real * b.imaginary + b.real * a.imaginary;
+
+	return ans;
+}
+
+COMPLEX complex_div ( COMPLEX n, COMPLEX d ) {
+	COMPLEX ans;
+	// find overall divisor
+	//
+	double oad;
+
+	oad  = d.real * d.real + d.imaginary * d.imaginary;
+	oad  = 1.0/oad;
+	//printf("\noad == %f numerator ",oad);
+	//complex_print (n);
+	//printf("denominator ",oad);
+	//complex_print (d);
+
+
+	d.imaginary = - d.imaginary; // conjugate
+	ans = complex_mul ( n, d );
+	complex_scalar_mul ( oad, &ans);
+
+	//printf("ans ");
+	//complex_print (ans);
+	//printf("\n");
+	return ans;
+}
+
+main () {
+
+	COMPLEX x,x_num,x_den, x_ans_den, x_ans_num, ccc;
+	int i,j,k; // counters
+	double r; // angle in radians
+
+	for (k=0;k<10;k++) {  // zero all the Z parameters
+		xz[k] = 0.0;
+		yp[k] = 0.0;
+	}
+
+	// simple lag filter
+	// xz[0] = 0.125;
+	// yp[0] = 1.0;
+	// yp[1] = -0.875;
+	//
+	//
+	//simple 7/8th lag filter squared
+	xz[0] = (0.125*0.125);
+	yp[0] = 1.0;
+	yp[1] = -(2.0*0.875);
+	yp[2] = (0.875*0.875); 
+
+	for (i=0;i<NYQUIST;i++) { 
+
+		r = (double)i/NYQUIST * PI;
+
+		x.real = cos(r);
+		x.imaginary = sin(r);
+
+		x_ans_num.real = x_ans_num.imaginary = 0.0;
+		x_ans_den.real = x_ans_den.imaginary = 0.0;
+                
+		// now take all the z transforms to their powers and multiple by coeffecients 
+		// to get the denominator and numerator for the frequency under inspection
+		for (j=0;j<10;j++) {
+
+			x_num = x;
+			complex_pow((double)-j,&x_num);
+			complex_scalar_mul(xz[j],&x_num);
+			complex_add ( &x_ans_num, x_num);
+
+			x_den = x;
+			complex_pow((double)-j,&x_den);
+			complex_scalar_mul(yp[j],&x_den);
+			complex_add ( &x_ans_den, x_den);
+
+		}
+
+		// OK now perform the division to obtain the frequency response 
+		//
+		ccc = complex_div ( x_ans_num, x_ans_den) ;
+
+		//ccc = x;
+		//complex_pow(2.0,&ccc);
+		printf(" %d arg = %f mag %f maglog10 %f\n", i, arg(ccc), mag(ccc), log(mag(ccc))*10.0 );
+	}
+
+	//----------------------------------------------------------------------------------------------------------
+	// Test pow and scalar multiply
+	/*
+	x.imaginary = 0.0;
+	x.real = 2.0;
+	complex_pow(3.0,&x);
+	printf(" x.r %f  x.i %f\n", x.real, x.imaginary);
+	complex_scalar_mul(2.1,&x);
+	printf(" x.r %f  x.i %f\n", x.real, x.imaginary);
+	printf(" after scalar * 2.1 x.r %f  x.i %f\n", x.real, x.imaginary);
+
+	x.imaginary = 1.0;
+	x.real = 0.0;
+	complex_pow(3.0,&x);
+	printf(" x.r %f  x.i %f\n", x.real, x.imaginary);
+	complex_scalar_mul(-2.1,&x);
+	printf(" after scalar * -2.1 x.r %f  x.i %f\n", x.real, x.imaginary);
+	*/
+	//----------------------------------------------------------------------------------------------------------
+	//
+	x_ans_num.real =  1.0; x_ans_num.imaginary = 10.0;
+	x_ans_den.real = 2.0; x_ans_den.imaginary = 80.0;
+	x = complex_div ( x_ans_num, x_ans_den );
+	printf(" after divide  x.r %f  x.i %f\n", x.real, x.imaginary);
+
+	
+
+}
+