1 //Spencer Jackson
2 //waves.c
3
4 #include<waves.h>
5 #include<stdlib.h>
6 #include<time.h>
7 #include<math.h>
8
9 #define _GNU_SOURCE
10
init_waves(WAVESOURCE * self)11 void init_waves(WAVESOURCE* self)
12 {
13 uint16_t i =0;
14 uint8_t j;
15 int8_t k=0;
16 double phase = 0;
17 self->half_phase = PI;
18 self->saw_step = 2*PI/TABLE_LENGTH;
19 self->phase_coeff = TABLE_LENGTH/(2*PI);
20 self->phase_offset = TABLE_LENGTH/2;
21
22 //saw
23 for(i=0; i<TABLE_LENGTH; i++)
24 {
25 self->saw_table[i] = 0;
26 k=1;
27 for(j=1; j<=MAX_N_HARMONICS; j++)
28 {
29 self->saw_table[i] += k*sin(j*phase)/(j);
30 k = -k;
31 }
32 phase += self->saw_step;
33 self->saw_table[i] *=.56694;
34 }
35
36 //tri
37 for(i=0; i<TABLE_LENGTH; i++)
38 {
39 self->tri_table[i] = 0;
40 k=1;
41 for(j=1; j<=MAX_N_HARMONICS; j+=2)
42 {
43 self->tri_table[i] += k*sin(j*phase)/(j*j);
44 k = -k;
45 }
46 phase += self->saw_step;
47 self->tri_table[i] *= .82922;
48 }
49
50 //white and random
51 srand ((uint16_t) time (NULL));
52 self->V = 2*rand() / (float)RAND_MAX - 1;
53 self->V2 = self->V*self->V;
54
55 //pointers
56 self->sin_func = &mySin;
57 self->saw_func = &saw;
58 self->sqr_func = □
59 self->tri_func = ▵
60 self->white_func = &white;
61 self->rand_func = &randomsnh;
62
63 self->wave_func[FUNC_SIN] = self->sin_func;
64 self->wave_func[FUNC_SAW] = self->saw_func;
65 self->wave_func[FUNC_SQR] = self->sqr_func;
66 self->wave_func[FUNC_TRI] = self->tri_func;
67 self->wave_func[FUNC_WHITE] = self->white_func;
68 self->wave_func[FUNC_RAND] = self->rand_func;
69
70
71 self->func_max = PI;
72 self->func_min = -PI;
73 self->func_domain = self->func_max - self->func_min;
74 }
75
init_hysteresis(HYSTERESIS * self)76 void init_hysteresis(HYSTERESIS *self)
77 {
78 self->prev_val = 2*rand() / (float)RAND_MAX - 1;
79 self->prev_phase = 0;
80 }
81
82 //8th order series approximation of pow(2,x) suggested by benjamin guihaire
myPow2(double x)83 double myPow2(double x)
84 {
85 double p;
86 if(x>=0) p = (1+x *0.00270760617406228636491106297445); //LN2/256 = 0.00270760617406228636491106297445
87 else p = 1/(1-x*0.00270760617406228636491106297445);
88 p *=p;
89 p *=p;
90 p *=p;
91 p *=p;
92 p *=p;
93 p *=p;
94 p *=p;
95 return p*p;
96 }
97
98 //based on an algorithm by Nicolas Capens
mySin(WAVESOURCE * self,HYSTERESIS * mem,double x)99 double mySin(WAVESOURCE *self, HYSTERESIS *mem, double x)
100 {
101 double y = 1.27323954474*x - 0.40528473456*x*(x>0?x:-x);
102 return 0.225*(y*(y>0?y:-y) - y) + y;
103 }
104
saw(WAVESOURCE * self,HYSTERESIS * mem,double phase)105 double saw(WAVESOURCE* self, HYSTERESIS* mem, double phase)
106 {
107 phase = phase*self->phase_coeff + self->phase_offset;
108 uint16_t hi,lo;
109 lo = (uint16_t) phase;
110 hi = lo +1;
111 return self->saw_table[lo] + (phase - lo)*(self->saw_table[hi] - self->saw_table[lo]);
112 }
113
square(WAVESOURCE * self,HYSTERESIS * mem,double phase)114 double square(WAVESOURCE* self, HYSTERESIS* mem, double phase)
115 {
116 if(phase > 0)
117 {
118 return saw(self,mem,phase) - saw(self,mem,phase - self->half_phase);
119 }
120 else
121 {
122 return saw(self,mem,phase) - saw(self,mem,phase + self->half_phase);
123 }
124 }
125
triangle(WAVESOURCE * self,HYSTERESIS * mem,double phase)126 double triangle(WAVESOURCE* self, HYSTERESIS *mem, double phase)
127 {
128 phase = phase*self->phase_coeff + self->phase_offset;
129 uint16_t hi,lo;
130 lo = (uint16_t) phase;
131 hi = lo +1;
132 return self->tri_table[lo] + (phase - lo)*(self->tri_table[hi] - self->tri_table[lo]);
133 }
134
135 //normal distribution approximation calculated by a modified Marsaglia polar method
white(WAVESOURCE * self,HYSTERESIS * mem,double phase)136 double white(WAVESOURCE* self, HYSTERESIS *mem, double phase)
137 {
138 float U = 2.0* rand() / (float)RAND_MAX - 1;//rand E(-1,1)
139 float S = U*U + self->V2;//map 2 random vars to unit circle
140
141 if(S>=1)//repull RV if outside unit circle
142 {
143 U = 2.0* rand() / (float)RAND_MAX - 1;
144 S = U*U + self->V2;
145 if(S>=1)
146 {
147 U = 2.0* rand() / (float)RAND_MAX - 1;
148 S = U*U + self->V2;
149 if(S>=1)
150 {
151 //guarantee an exit, value will be unchanged
152 U=0;
153 }
154 }
155 }
156
157 if(U)
158 {
159 self->V = U;//store V for next round
160 self->V2 = U*U;
161 return U*sqrt(-2*log(S)/S);
162 }
163 else
164 {
165 return self->V;
166 }
167
168
169 }
170
171 //currently a uniform distribution
randomsnh(WAVESOURCE * self,HYSTERESIS * mem,double phase)172 double randomsnh(WAVESOURCE* self, HYSTERESIS* mem, double phase)
173 {
174 if(phase < mem->prev_phase)
175 {
176 mem->prev_val = 2*rand() / (float)RAND_MAX - 1;
177 }
178 mem->prev_phase = phase;
179 return mem->prev_val;
180 }
181