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