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...
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