Project

General

Profile

Revision 1602

Added some test files for fixed point math.

Something definitely broke in cosine when I changed to smaller tables
(which shouldn't have sacrificed any accuracy). Need to investigate...

View differences:

trunk/code/projects/fp_math/fp_math.c
9 9
#define FP_TWO_PI 411775
10 10

  
11 11
#define TABLE_LENGTH 64
12
#define TABLE_STEP 1621//.160757
13
#define EXP_MAGICONSTANT
12
#define TABLE_STEP 1635//.160757
13
#define EXP_MAGICONSTANT 726817 //ln(2^16) * 2^16
14 14

  
15
//#define DEBUG
16

  
17
//TODO This can probably be removed at some point.
15 18
static int32_t linspace[TABLE_LENGTH] PROGMEM = { 
16
	0, 1621, 3242, 4863, 6485, 8106, 9727, 11348, 
17
	12969, 14590, 16212, 17833, 19454, 21075, 22696, 24317, 
18
	25939, 27560, 29181, 30802, 32423, 34044, 35666, 37287, 
19
	38908, 40529, 42150, 43771, 45393, 47014, 48635, 50256, 
20
	51877, 53498, 55119, 56741, 58362, 59983, 61604, 63225, 
21
	64846, 66468, 68089, 69710, 71331, 72952, 74573, 76195, 
22
	77816, 79437, 81058, 82679, 84300, 85922, 87543, 89164, 
23
	90785, 92406, 94027, 95648, 97270, 98891, 100512, 102133
19
	0, 1634, 3268, 4902, 6536, 8170, 9804, 11438, 
20
	13072, 14706, 16340, 17974, 19608, 21242, 22876, 24510, 
21
	26144, 27778, 29412, 31047, 32681, 34315, 35949, 37583, 
22
	39217, 40851, 42485, 44119, 45753, 47387, 49021, 50655, 
23
	52289, 53923, 55557, 57191, 58825, 60459, 62093, 63727, 
24
	65361, 66995, 68629, 70263, 71897, 73531, 75165, 76799, 
25
	78433, 80067, 81701, 83335, 84969, 86603, 88237, 89871, 
26
	91506, 93140, 94774, 96408, 98042, 99676, 101310, 102944, 
24 27
};
25 28

  
26 29
static int32_t fp_cos_table[TABLE_LENGTH] PROGMEM = {
27
	65536, 65516, 65456, 65356, 65215, 65035, 64815, 64556, 
28
	64257, 63919, 63541, 63125, 62670, 62176, 61645, 61076, 
29
	60470, 59826, 59146, 58430, 57678, 56890, 56068, 55212, 
30
	54322, 53398, 52442, 51454, 50434, 49384, 48303, 47193, 
31
	46053, 44886, 43691, 42470, 41222, 39949, 38652, 37331, 
32
	35988, 34622, 33235, 31828, 30401, 28956, 27492, 26013, 
33
	24517, 23006, 21481, 19943, 18393, 16831, 15260, 13679, 
34
	12089, 10492, 8889, 7280, 5667, 4050, 2431, 811
30
	65536, 65516, 65455, 65353, 65210, 65027, 64804, 64540, 
31
	64237, 63893, 63509, 63087, 62624, 62123, 61584, 61006, 
32
	60390, 59736, 59046, 58319, 57555, 56756, 55921, 55052, 
33
	54148, 53211, 52241, 51238, 50203, 49138, 48041, 46915, 
34
	45760, 44576, 43364, 42126, 40861, 39571, 38256, 36918, 
35
	35556, 34173, 32768, 31343, 29898, 28435, 26954, 25456, 
36
	23943, 22415, 20872, 19317, 17750, 16171, 14583, 12986, 
37
	11380, 9768, 8149, 6525, 4898, 3267, 1634, 0, 
35 38
};
36 39

  
37 40
//FIXME Lazy implementations of tangent and sine.
......
48 51
	int32_t xsquare = fp_mult(x,x);
49 52
	int32_t result, i;
50 53
	
51
	if(theta > EXP_MAGICONSTANT) { 
54
	if(x > EXP_MAGICONSTANT) { 
52 55
		usb_puts("Output value is out of range.\n\r");
53 56
		return -1;
54 57
	}
55 58
	
56 59
	//This is again, massive amounts of lazyness.
57 60
	//The approximation fails outside this range (approximately).
58
	if(theta < -17)
61
	if(x < (-17l << 16))
59 62
		return 0;
60 63
	
61 64
	result = xsquare;
62
	for(i= (22 << 16); i > (2 << 16); i-= (4 << 16)) { 
65
	for(i= (22l << 16); i > (2l << 16); i-= (4l << 16)) { 
63 66
		result += i;
64 67
		result = fp_div(xsquare, result);
65 68
	}
66 69
	
67
	result += (2 << 16) - x;
68
	return (1 << 16) + fp_div(fp_mult((2 << 16), x), result);
70
	result += (2l << 16) - x;
71
	return (1l << 16) + fp_div(fp_mult((2l << 16), x), result);
69 72

  
70 73
}
71 74

  
......
76 79
	int32_t x_n, x_np1, x_np2;
77 80
	int32_t y_n, y_np1, y_np2;
78 81
	int32_t dd_n, dd_np1, second_dd, result;
79
	
82

  
83
#ifdef DEBUG
84
	char buf[128];
85
#endif
86

  
80 87
	//Move theta into [0, pi/2] w/ appropriate sign.
81 88
	theta = ABS(theta) % FP_TWO_PI; 
82 89

  
......
90 97

  
91 98
	//Find the nearest table values. FIXME
92 99
	n = theta / TABLE_STEP;
93
	while( n < TABLE_LENGTH - 1 && (x_np1 = pgm_read_dword(&linspace[n+1]))) 
100
	while( n < TABLE_LENGTH - 1 
101
		&& (x_np1 = pgm_read_dword(&linspace[n+1]) < theta)) 
94 102
		n++;
95 103
	
96 104
	//theta is between x_{n} and x_{n+1}
105

  
97 106
	if(n == TABLE_LENGTH - 1) { 
98 107
		//Perform linear interpolation, since we're close to zero anyway.
99 108
		x_n = pgm_read_dword(&linspace[TABLE_LENGTH - 1]);
100 109
		y_n = pgm_read_dword(&fp_cos_table[TABLE_LENGTH - 1]);
101 110

  
102
		result = fp_mult(fp_div(FP_PI_OVER_TWO - x_n, 0 - y_n), theta - x_n) + y_n;
111
		result = fp_mult(
112
			fp_div(FP_PI_OVER_TWO - x_n, 0 - y_n), theta - x_n) + y_n;
113
		
103 114
		return negative ? -result : result;
104 115
	}
105 116

  
106 117
	if(n == TABLE_LENGTH) { 
107 118
		//We didn't find a value! Oh no!
108
		usb_puts("fp_math: We couldn't find surrounding table values! \n\r
109
				  This should never happen!!!\n\r\n\r");
119
		usb_puts("fp_math: Fatal! We couldn't find surrounding table values! \n\r");
120
				  
110 121
		return 0;
111 122
	}
112
	
123

  
113 124
	//Address the general case. Quadratic interpolation.
114 125
	//Load in the necessary values.
115 126
	x_n = pgm_read_dword(&linspace[n]);
127
	x_np1 = pgm_read_dword(&linspace[n + 1]);
116 128
	x_np2  = pgm_read_dword(&linspace[n + 2]);
117 129
	
118 130
	y_n = pgm_read_dword(&fp_cos_table[n]);
119 131
	y_np1 = pgm_read_dword(&fp_cos_table[n + 1]);
120 132
	y_np2 = pgm_read_dword(&fp_cos_table[n + 2]);
121 133

  
134
#ifdef DEBUG
135
	sprintf(buf, "x_n = %ld, theta = %ld, x_{n+1} = %ld", x_n, theta, x_np1);
136
	usb_puts(buf);
137
#endif
138

  
122 139
	//Calculate first divided differences.
123 140
	dd_n = fp_div(y_np1 - y_n, x_np1 - x_n);
124 141
	dd_np1 = fp_div(y_np2 - y_np1, x_np2 - x_np1);
trunk/code/projects/fp_math/gentables.m
1 1

  
2
space = int32(round(linspace(0,pi,128) * 2^16))
3
costable = int32(round(cos(linspace(0,pi,128)) * 2^16));
4
sintable = int32(round(sin(linspace(0,pi,128)) * 2^16));
2
TABLE_LENGTH = 64;
5 3

  
4
space = int32(round(linspace(0, pi/2, TABLE_LENGTH) * 2^16));
5
costable = int32(round(cos(linspace(0, pi/2, TABLE_LENGTH)) * 2^16));
6

  
6 7
space_file = fopen('space_table.txt', 'w');
7 8
cos_file = fopen('cos_table.txt', 'w');
8
sin_file = fopen('sin_table.txt', 'w');
9 9

  
10 10
for i=1:length(space)
11 11
	fprintf(space_file, '%d, ', space(i));
......
20 20
		fprintf(cos_file, '\n');
21 21
	end
22 22
end
23

  
24
for i=1:length(sintable)
25
	fprintf(sin_file, '%d, ', sintable(i));
26
	if ~mod(i,8)
27
		fprintf(sin_file, '\n');
28
	end
29
end
30

  
trunk/code/projects/fp_math/test/serial.c
1
#include <stdio.h>
2
#include <stdlib.h>
3
#include <string.h>
4
#include <unistd.h>
5
#include <fcntl.h>
6
#include <errno.h>
7
#include <termios.h>
8
#include <unistd.h>
9
#include "serial.h"
10

  
11
/* open and configure the port */
12
int openport(char* portname){
13
  struct termios oldtio, newtio;
14
  char *error = 0;
15
  int fd;
16

  
17
  /* open read/write, non-controlling, non-blocking */
18
  fd = open(portname, O_RDWR | O_NOCTTY | O_NONBLOCK);
19
  if (fd == -1){
20
    fprintf(stderr, "Unable to open %s: ", portname);
21
    perror(error);
22
    exit(-1);
23
  }
24

  
25
  tcgetattr(fd, &oldtio);          //save current port settings
26
  bzero(&newtio, sizeof(newtio));  //clear struct for new settings
27

  
28
  newtio.c_cflag = B115200 | CS8 | CLOCAL | CREAD;
29
  newtio.c_iflag = IGNPAR | ICRNL;
30

  
31
  newtio.c_oflag = 0;
32
  newtio.c_lflag = 0;
33

  
34
  newtio.c_cc[VINTR]    = 0;     /* Ctrl-c */
35
  newtio.c_cc[VQUIT]    = 0;     /* Ctrl-\ */
36
  newtio.c_cc[VERASE]   = 0;     /* del */
37
  newtio.c_cc[VKILL]    = 0;     /* @ */
38
  newtio.c_cc[VEOF]     = 4;     /* Ctrl-d */
39
  newtio.c_cc[VTIME]    = 0;     /* inter-character timer unused */
40
  newtio.c_cc[VMIN]     = 1;     /* blocking read until 1 character arrives 
41
*/
42
  newtio.c_cc[VSWTC]    = 0;     /* '\0' */
43
  newtio.c_cc[VSTART]   = 0;     /* Ctrl-q */
44
  newtio.c_cc[VSTOP]    = 0;     /* Ctrl-s */
45
  newtio.c_cc[VSUSP]    = 0;     /* Ctrl-z */
46
  newtio.c_cc[VEOL]     = 0;     /* '\0' */
47
  newtio.c_cc[VREPRINT] = 0;     /* Ctrl-r */
48
  newtio.c_cc[VDISCARD] = 0;     /* Ctrl-u */
49
  newtio.c_cc[VWERASE]  = 0;     /* Ctrl-w */
50
  newtio.c_cc[VLNEXT]   = 0;     /* Ctrl-v */
51
  newtio.c_cc[VEOL2]    = 0;     /* '\0' */
52

  
53
  /* flush the port, and apply the new settings now */
54
  tcflush(fd, TCIFLUSH);
55
  tcsetattr(fd, TCSANOW, &newtio);
56

  
57
  return (fd);
58
}
59

  
60
int putbyte (int fd, char c)
61
{
62
   int result;
63
   result= write(fd, &c, 1);
64
//   Sleep(25);
65
   return result;
66
}
67

  
68
unsigned char getbyte (int fd)
69
{
70
  char c;
71
  ssize_t retval;
72

  
73
  retval = read(fd, &c, 1); 
74

  
75
  if (retval <= 0)
76
    return 255;
77

  
78
  return c;
79
}
80

  
81
unsigned char getbyteblock (int fd)
82
{
83
  char c;
84
  ssize_t retval;
85

  
86
  while ((retval = read(fd, &c, 1)) <= 0){}
87
//    usleep(5000);
88

  
89
  return c;
90
}
91

  
92
void clearbuf(int fd)
93
{
94
  tcflush(fd, TCIFLUSH);
95
}
96

  
trunk/code/projects/fp_math/test/serial.h
1
#ifndef __SERIAL_H
2

  
3
int openport(char *portname);
4
int putbyte (int fd, char c);
5
unsigned char getbyte (int fd);
6
unsigned char getbyteblock (int fd);
7
void clearbuf(int fd);
8

  
9
#define __SERIAL_H
10
#endif
trunk/code/projects/fp_math/test/server_test.c
1

  
2
#include <math.h>
3
#include <stdio.h>
4
#include <stdlib.h>
5
#include "serial.h"
6

  
7
#define X_IN 0
8
#define Y_OUT 1
9

  
10
#define FP_COS 2
11
#define FP_SIN 3
12
#define FP_TAN 4
13

  
14
#define EXIT 5
15

  
16
typedef union fp_t { 
17
	char d[4];
18
	int v;
19
} fp_t;
20

  
21
int main(int argc, char** argv) { 
22
	
23
	int fd, i, n = 0;
24
	float input, output, error, avg_error, max_error;
25
	fp_t fp_input, fp_output;
26
	
27
	
28
	fd = openport("/dev/ttyUSB0");
29
	
30
	printf("Beginning test loop.");
31

  
32
	while(1) {
33

  
34
		switch(getbyte(fd)) { 
35
		
36
		case X_IN:
37
			for(i=0; i < 4; i++)
38
				fp_input.d[i] = getbyte(fd);
39
			
40
			input = ((float)fp_input.v) / 65536.;
41

  
42
			break;
43

  
44
		case Y_OUT:
45
			for(i=0; i < 4; i++)
46
				fp_output.d[i] = getbyte(fd);
47
			
48
			output = ((float)fp_output.v) / 65536.;
49

  
50
			break;
51

  
52
		case FP_COS: 
53

  
54
			n++;
55
			printf("%f = cos(%f)\n", output, input);
56
			error = abs(cos(input) - output);
57
			if(error > max_error) max_error = error;
58
			avg_error = (avg_error * (n - 1) + error) / n;
59
			
60
			break;
61

  
62
		case FP_SIN:
63
			
64
			break;
65

  
66
		case FP_TAN:
67
			
68
			break;
69

  
70
		case EXIT:
71
		
72
			printf("Average error: %f\n", avg_error);
73
			printf("Max error: %f\n", max_error);
74
			exit(EXIT_SUCCESS);
75

  
76
			break;
77

  
78
		default: break;
79
			printf("Uh oh - defauled\n");
80
			
81
			break;
82

  
83
		}
84
	}
85
}
86

  
87

  
88

  
trunk/code/projects/fp_math/test/Makefile
1

  
2
all:  serial.o
3
	gcc server_test.c serial.o -lm -o server_test
4

  
5
clean: 
6
	rm *.o server_test
7

  
trunk/code/projects/fp_math/space_table.txt
1
0, 1634, 3268, 4902, 6536, 8170, 9804, 11438, 
2
13072, 14706, 16340, 17974, 19608, 21242, 22876, 24510, 
3
26144, 27778, 29412, 31047, 32681, 34315, 35949, 37583, 
4
39217, 40851, 42485, 44119, 45753, 47387, 49021, 50655, 
5
52289, 53923, 55557, 57191, 58825, 60459, 62093, 63727, 
6
65361, 66995, 68629, 70263, 71897, 73531, 75165, 76799, 
7
78433, 80067, 81701, 83335, 84969, 86603, 88237, 89871, 
8
91506, 93140, 94774, 96408, 98042, 99676, 101310, 102944, 
trunk/code/projects/fp_math/main.c
1 1
#include <serial.h>
2 2
#include <stdint.h>
3
#include <stdlib.h>
3 4
#include "fp_math.h"
4 5

  
6
#define X_IN 0
7
#define Y_OUT 1
8

  
9
#define FP_COS 2
10
#define FP_SIN 3
11
#define FP_TAN 4
12

  
13
#define EXIT 5
14

  
5 15
int main(int argc, char** argv) { 
6
	int32_t i,j;
16
	int32_t i, x, y;
7 17
	char buf[128];
18

  
19
	usb_init();
8 20
	
9
	usb_init();
10
	usb_puts("USB Initialized\n\r");
21
	for(i=0; i < 10000; i++) { 
22

  
23
		//Random number / sign
24
		x = random();
25
		x = ((rand() >> 2) & 1) ? x : -x;
11 26
		
12
	for(i=10000; i < 11000; i+=5) { 
13
		sprintf(buf, "fp_cos(%ld) = %ld\n\r", i,  fp_cos(i));
14
		usb_puts(buf);
27
		usb_putc(X_IN);
28
		for(i=0; i < 32; i += 8) 
29
			usb_putc((x >> i)  & 0xFF);
30

  
31
		y = fp_cos(x);	
32
		
33
		usb_putc(Y_OUT);
34
		for(i=0; i < 32; i += 8) 
35
			usb_putc((y >> i)  & 0xFF);
36
	
37
		usb_putc(FP_COS);
38
		
39
		//sprintf(buf, "%ld = cos(%ld)\n\r", y, x);
40
		//usb_puts(buf);
15 41
	}
42
	
43
	usb_putc(EXIT);
44

  
45
	while(1);
16 46
}
17 47

  
trunk/code/projects/fp_math/cos_table.txt
1
65536, 65516, 65455, 65353, 65210, 65027, 64804, 64540, 
2
64237, 63893, 63509, 63087, 62624, 62123, 61584, 61006, 
3
60390, 59736, 59046, 58319, 57555, 56756, 55921, 55052, 
4
54148, 53211, 52241, 51238, 50203, 49138, 48041, 46915, 
5
45760, 44576, 43364, 42126, 40861, 39571, 38256, 36918, 
6
35556, 34173, 32768, 31343, 29898, 28435, 26954, 25456, 
7
23943, 22415, 20872, 19317, 17750, 16171, 14583, 12986, 
8
11380, 9768, 8149, 6525, 4898, 3267, 1634, 0, 
trunk/code/projects/fp_math/Makefile
12 12

  
13 13
# com1 = serial port. Use lpt1 to connect to parallel port.
14 14
AVRDUDE_PORT = $(shell if uname -s |grep -i w32 >/dev/null; then echo 'COM4:'; else echo '/dev/ttyUSB0'; fi)
15
AVRDUDE_PORT = /dev/tty.usbserial-A4001gUV
15
AVRDUDE_PORT = /dev/ttyUSB0
16 16

  
17 17

  
18 18
else

Also available in: Unified diff