two pole files with int sizes specified

HP48_notes.txt

@@ -0,0 +1,33 @@
+
+
+To get an HP48 ploting the frequerncy response of a 
+Z transform, first define Z as 'EXP(-(0,1)).
+
+i.e.
+
+> 'EXP(-(0,1))'
+> 'Z'
+> STO
+
+Now define a function to convert a Z transform
+to ABS and then to DBs 
+
+> << ABS LOG10 10 * >>
+> '->ABS'
+> STO
+
+Now place in a test function say the 7/8ths LAG FILTER
+
+> ' 0.125/(1-0.875*Z)'
+> ->ABS
+
+Now store in 'EQ'
+
+Now get the plot menu.
+Set yrange to -40 to + 40
+Set xrange to 0 to 3.142
+
+PLOT func
+
+The Z transform magnitude response will now be displayed in the 
+HP48 graph

Makefile

@@ -0,0 +1,6 @@
+
+
+two_pole:
+	gcc -o tp two_pole_7_8.c -lm
+	./tp > tp.dat
+	gnuplot < tp.gpt

tp.gpt

@@ -0,0 +1,9 @@
+
+
+#plot "tp.dat" using 1:3 title "input", "tp.dat" using 1:5 title "two pole", "tp.dat" using 1:6 title "34LAG", "tp.dat" using 1:7 title "78LAG"
+plot "tp.dat" using 1:3 title "input" with linespoints, "tp.dat" using 1:5 title "two pole" with linespoints,  "tp.dat" using 1:7 title "78LAG" with linespoints
+
+!sleep 10
+!sleep 10
+!sleep 10
+!sleep 90

tp_2.gpt

@@ -0,0 +1,10 @@
+
+set xrange[200:1000]
+
+#plot "tp.dat" using 1:3 title "input", "tp.dat" using 1:5 title "two pole", "tp.dat" using 1:6 title "34LAG", "tp.dat" using 1:7 title "78LAG"
+plot "tp.dat" using 1:3 title "input" with linespoints, "tp.dat" using 1:5 title "two pole" with linespoints,  "tp.dat" using 1:7 title "78LAG" with linespoints
+
+!sleep 10
+!sleep 10
+!sleep 10
+!sleep 90

two_pole_7_8.c

@@ -0,0 +1,108 @@
+
+
+
+#include <stdint.h>
+#include <stdio.h>
+#include <math.h>
+
+// Accuracy here matters.
+// With a max 25k TDS (i.e. 16 bits input signal signed)
+// this means 6 + 16 == 22.
+// This fits within a 32 bit int
+// and leaves ten bits of head room.
+//
+// For two pole the y0, y1, y2 registers must be 32 bit at least.
+// The new rpi is 64 bit so these have had to have been deliberately set
+// as int32_t to simulate the pic18 `long' type.
+// The short int type on the PIC18 is `int'
+//
+#define BIN_FRACS 6
+
+// this is a 64 bit machine (the pi)
+//
+int16_t /* squared version of LAG_7_8 */
+two_pole_7_8 ( int16_t input ) {
+
+	static int32_t y1=0,y2=0;
+	int32_t * res;
+	int32_t y0;
+
+	int32_t x0 = input;
+
+	// x0 times 0.125	
+	// the minus 3 divides by 8	 : DOUBLE POLE
+	x0 <<= (BIN_FRACS-3) ;   // now all calculations are done times BIN_FRACS^2
+
+	// this works really well WHY?
+	//y0 = x0 + y1-(y1>>3) - ((y2>>1) + (y2>>2) + (y2>>6));
+	
+        // should be 	
+	//y0 = x0 + y1+y1-(y1>>2) - (y2>>1 + y2>>2 + y2>>6);
+	//
+	//
+	// OK this works well 
+	y0 = x0 + (14.0 / 8.0) * (double) y1 - (49.0/64.0) * (double) y2; 
+
+
+
+	y0 = x0 +
+	       //	(14.0 / 8.0) * (double) y1
+	       //	this is 7/4
+	       y1 + y1 - (y1>>2)  
+	       //
+	       //
+	       //	- (49.0/64.0) * (double) y2; 
+	       //
+	       -  (    (y2>>1) +  // half
+	               (y2>>2) +  // quarter
+                       (y2>>6) ); // 64th      
+
+
+	//
+	// y0 = x0/64 + (y1 * 7) / 4 + (y2*49)/64;
+
+	// see what happens when y2>>6 is left out. I theory d.c goes unstable
+	// and yes it does! 25AUG2019
+	// y0 = x0 + y1-(y1>>3) - (y2>>1 + y2>>2  /* + y2>>6*/);
+
+	y2 = y1;
+	y1 = y0;
+
+	//res = (short int *)&y0;
+	//res++; 
+	//printf("y0=%X res should be %X res=%x\n",y0, y0>>16, res);
+	//
+
+	// the plus 3 divides by 8	
+	return y0>>(BIN_FRACS+3); // divide back down for scaling and then divide by 8 filter gain // *res;
+}
+
+
+
+int main () {
+
+	int i;
+	int16_t val,res, res34, res78;
+
+	for (i=0;i<1000;i++) {
+
+		// ramp + some sine
+		val = sin ( (double)i*10.0  / (3.142 * 2.0)  ) * 10.0 + 25000 ;
+
+		// fast sine
+		//val = sin ( ((double)i*10.0)  / (3.142 * 2.0)  ) * 10000;
+
+
+		res = two_pole_7_8 ( val );
+
+		res34 = (((res34<<2) - res34)>>2) + (val>>2);
+		res78 = (((res78<<3) - res78)>>3) + (val>>3);
+
+		printf ("%d val %d res %d res34 %d res78 %d\n",i, val,res, res34, res78);
+
+	}
+
+	return 0;
+}
+	
+