tsiry-sandratraina.com
A modern Music Player Daemon based on Rockbox open source high quality audio player
libadwaita
audio
rust
zig
deno
mpris
rockbox
mpd
1/*
2** Algorithm: complex multiplication
3**
4** Performance: Bad accuracy, great speed.
5**
6** The simplest, way of generating trig values. Note, this method is
7** not stable, and errors accumulate rapidly.
8**
9** Computation: 2 *, 1 + per value.
10*/
11
12#include "m_pd.h"
13
14#include "m_fixed.h"
15
16
17#define TRIG_FAST
18
19#ifdef TRIG_FAST
20char mtrig_algorithm[] = "Simple";
21#define TRIG_VARS \
22 t_fixed t_c,t_s;
23#define TRIG_INIT(k,c,s) \
24 { \
25 t_c = fcostab[k]; \
26 t_s = fsintab[k]; \
27 c = itofix(1); \
28 s = 0; \
29 }
30#define TRIG_NEXT(k,c,s) \
31 { \
32 t_fixed t = c; \
33 c = mult(t,t_c) - mult(s,t_s); \
34 s = mult(t,t_s) + mult(s,t_c); \
35 }
36#define TRIG_23(k,c1,s1,c2,s2,c3,s3) \
37 { \
38 c2 = mult(c1,c1) - mult(s1,s1); \
39 s2 = (mult(c1,s1)<<2); \
40 c3 = mult(c2,c1) - mult(s2,s1); \
41 s3 = mult(c2,s1) + mult(s2,c1); \
42 }
43#endif
44#define TRIG_RESET(k,c,s)
45
46/*
47** Algorithm: O. Buneman's trig generator.
48**
49** Performance: Good accuracy, mediocre speed.
50**
51** Manipulates a log(n) table to stably create n evenly spaced trig
52** values. The newly generated values lay halfway between two
53** known values, and are calculated by appropriatelly scaling the
54** average of the known trig values appropriatelly according to the
55** angle between them. This process is inherently stable; and is
56** about as accurate as most trig library functions. The errors
57** caused by this code segment are primarily due to the less stable
58** method to calculate values for 2t and s 3t. Note: I believe this
59** algorithm is patented(!), see readme file for more details.
60**
61** Computation: 1 +, 1 *, much integer math, per trig value
62**
63*/
64
65#ifdef TRIG_GOOD
66#define TRIG_VARS \
67 int t_lam=0; \
68 double coswrk[TRIG_TABLE_SIZE],sinwrk[TRIG_TABLE_SIZE];
69#define TRIG_INIT(k,c,s) \
70 { \
71 int i; \
72 for (i=0 ; i<=k ; i++) \
73 {coswrk[i]=fcostab[i];sinwrk[i]=fsintab[i];} \
74 t_lam = 0; \
75 c = 1; \
76 s = 0; \
77 }
78#define TRIG_NEXT(k,c,s) \
79 { \
80 int i,j; \
81 (t_lam)++; \
82 for (i=0 ; !((1<<i)&t_lam) ; i++); \
83 i = k-i; \
84 s = sinwrk[i]; \
85 c = coswrk[i]; \
86 if (i>1) \
87 { \
88 for (j=k-i+2 ; (1<<j)&t_lam ; j++); \
89 j = k - j; \
90 sinwrk[i] = halsec[i] * (sinwrk[i-1] + sinwrk[j]); \
91 coswrk[i] = halsec[i] * (coswrk[i-1] + coswrk[j]); \
92 } \
93 }
94#endif
95
96
97#define TRIG_TAB_SIZE 22
98
99#ifndef ROCKBOX
100static long long halsec[TRIG_TAB_SIZE]= {1,2,3};
101#endif
102
103#define FFTmult(x,y) mult(x,y)
104
105
106
107
108#include "d_imayer_tables.h"
109
110#define SQRT2 ftofix(1.414213562373095048801688724209698)
111
112
113static void imayer_fht(t_fixed *fz, int n)
114{
115 int k,k1,k2,k3,k4,kx;
116 t_fixed *fi,*fn,*gi;
117 TRIG_VARS;
118
119 for (k1=1,k2=0;k1<n;k1++)
120 {
121 t_fixed aa;
122 for (k=n>>1; (!((k2^=k)&k)); k>>=1);
123 if (k1>k2)
124 {
125 aa=fz[k1];fz[k1]=fz[k2];fz[k2]=aa;
126 }
127 }
128 for ( k=0 ; (1<<k)<n ; k++ );
129 k &= 1;
130 if (k==0)
131 {
132 for (fi=fz,fn=fz+n;fi<fn;fi+=4)
133 {
134 t_fixed f0,f1,f2,f3;
135 f1 = fi[0 ]-fi[1 ];
136 f0 = fi[0 ]+fi[1 ];
137 f3 = fi[2 ]-fi[3 ];
138 f2 = fi[2 ]+fi[3 ];
139 fi[2 ] = (f0-f2);
140 fi[0 ] = (f0+f2);
141 fi[3 ] = (f1-f3);
142 fi[1 ] = (f1+f3);
143 }
144 }
145 else
146 {
147 for (fi=fz,fn=fz+n,gi=fi+1;fi<fn;fi+=8,gi+=8)
148 {
149 t_fixed bs1,bc1,bs2,bc2,bs3,bc3,bs4,bc4,
150 bg0,bf0,bf1,bg1,bf2,bg2,bf3,bg3;
151 bc1 = fi[0 ] - gi[0 ];
152 bs1 = fi[0 ] + gi[0 ];
153 bc2 = fi[2 ] - gi[2 ];
154 bs2 = fi[2 ] + gi[2 ];
155 bc3 = fi[4 ] - gi[4 ];
156 bs3 = fi[4 ] + gi[4 ];
157 bc4 = fi[6 ] - gi[6 ];
158 bs4 = fi[6 ] + gi[6 ];
159 bf1 = (bs1 - bs2);
160 bf0 = (bs1 + bs2);
161 bg1 = (bc1 - bc2);
162 bg0 = (bc1 + bc2);
163 bf3 = (bs3 - bs4);
164 bf2 = (bs3 + bs4);
165 bg3 = FFTmult(SQRT2,bc4);
166 bg2 = FFTmult(SQRT2,bc3);
167 fi[4 ] = bf0 - bf2;
168 fi[0 ] = bf0 + bf2;
169 fi[6 ] = bf1 - bf3;
170 fi[2 ] = bf1 + bf3;
171 gi[4 ] = bg0 - bg2;
172 gi[0 ] = bg0 + bg2;
173 gi[6 ] = bg1 - bg3;
174 gi[2 ] = bg1 + bg3;
175 }
176 }
177 if (n<16) return;
178
179 do
180 {
181 t_fixed s1,c1;
182 int ii;
183 k += 2;
184 k1 = 1 << k;
185 k2 = k1 << 1;
186 k4 = k2 << 1;
187 k3 = k2 + k1;
188 kx = k1 >> 1;
189 fi = fz;
190 gi = fi + kx;
191 fn = fz + n;
192 do
193 {
194 t_fixed g0,f0,f1,g1,f2,g2,f3,g3;
195 f1 = fi[0 ] - fi[k1];
196 f0 = fi[0 ] + fi[k1];
197 f3 = fi[k2] - fi[k3];
198 f2 = fi[k2] + fi[k3];
199 fi[k2] = f0 - f2;
200 fi[0 ] = f0 + f2;
201 fi[k3] = f1 - f3;
202 fi[k1] = f1 + f3;
203 g1 = gi[0 ] - gi[k1];
204 g0 = gi[0 ] + gi[k1];
205 g3 = FFTmult(SQRT2, gi[k3]);
206 g2 = FFTmult(SQRT2, gi[k2]);
207 gi[k2] = g0 - g2;
208 gi[0 ] = g0 + g2;
209 gi[k3] = g1 - g3;
210 gi[k1] = g1 + g3;
211 gi += k4;
212 fi += k4;
213 } while (fi<fn);
214 TRIG_INIT(k,c1,s1);
215 for (ii=1;ii<kx;ii++)
216 {
217 t_fixed c2,s2;
218 TRIG_NEXT(k,c1,s1);
219 c2 = FFTmult(c1,c1) - FFTmult(s1,s1);
220 s2 = 2*FFTmult(c1,s1);
221 fn = fz + n;
222 fi = fz +ii;
223 gi = fz +k1-ii;
224 do
225 {
226 t_fixed a,b,g0,f0,f1,g1,f2,g2,f3,g3;
227 b = FFTmult(s2,fi[k1]) - FFTmult(c2,gi[k1]);
228 a = FFTmult(c2,fi[k1]) + FFTmult(s2,gi[k1]);
229 f1 = fi[0 ] - a;
230 f0 = fi[0 ] + a;
231 g1 = gi[0 ] - b;
232 g0 = gi[0 ] + b;
233 b = FFTmult(s2,fi[k3]) - FFTmult(c2,gi[k3]);
234 a = FFTmult(c2,fi[k3]) + FFTmult(s2,gi[k3]);
235 f3 = fi[k2] - a;
236 f2 = fi[k2] + a;
237 g3 = gi[k2] - b;
238 g2 = gi[k2] + b;
239 b = FFTmult(s1,f2) - FFTmult(c1,g3);
240 a = FFTmult(c1,f2) + FFTmult(s1,g3);
241 fi[k2] = f0 - a;
242 fi[0 ] = f0 + a;
243 gi[k3] = g1 - b;
244 gi[k1] = g1 + b;
245 b = FFTmult(c1,g2) - FFTmult(s1,f3);
246 a = FFTmult(s1,g2) + FFTmult(c1,f3);
247 gi[k2] = g0 - a;
248 gi[0 ] = g0 + a;
249 fi[k3] = f1 - b;
250 fi[k1] = f1 + b;
251 gi += k4;
252 fi += k4;
253 } while (fi<fn);
254 }
255 TRIG_RESET(k,c1,s1);
256 } while (k4<n);
257}
258
259
260void imayer_fft(int n, t_fixed *real, t_fixed *imag)
261{
262 t_fixed a,b,c,d;
263 t_fixed q,r,s,t;
264 int i,j,k;
265 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
266 a = real[i]; b = real[j]; q=a+b; r=a-b;
267 c = imag[i]; d = imag[j]; s=c+d; t=c-d;
268 real[i] = (q+t)>>1; real[j] = (q-t)>>1;
269 imag[i] = (s-r)>>1; imag[j] = (s+r)>>1;
270 }
271 imayer_fht(real,n);
272 imayer_fht(imag,n);
273}
274
275void imayer_ifft(int n, t_fixed *real, t_fixed *imag)
276{
277 t_fixed a,b,c,d;
278 t_fixed q,r,s,t;
279 int i,j,k;
280 imayer_fht(real,n);
281 imayer_fht(imag,n);
282 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
283 a = real[i]; b = real[j]; q=a+b; r=a-b;
284 c = imag[i]; d = imag[j]; s=c+d; t=c-d;
285 imag[i] = (s+r)>>1; imag[j] = (s-r)>>1;
286 real[i] = (q-t)>>1; real[j] = (q+t)>>1;
287 }
288}
289
290void imayer_realfft(int n, t_fixed *real)
291{
292#ifdef ROCKBOX
293 t_fixed a,b;
294#else
295 t_fixed a,b,c,d;
296#endif
297 int i,j,k;
298 imayer_fht(real,n);
299 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
300 a = real[i];
301 b = real[j];
302 real[j] = (a-b)>>1;
303 real[i] = (a+b)>>1;
304 }
305}
306
307void imayer_realifft(int n, t_fixed *real)
308{
309#ifdef ROCKBOX
310 t_fixed a,b;
311#else
312 t_fixed a,b,c,d;
313#endif
314 int i,j,k;
315 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
316 a = real[i];
317 b = real[j];
318 real[j] = (a-b);
319 real[i] = (a+b);
320 }
321 imayer_fht(real,n);
322}
323
324
325
326#ifdef TEST
327
328static double dfcostab[TRIG_TAB_SIZE]=
329 {
330 .00000000000000000000000000000000000000000000000000,
331 .70710678118654752440084436210484903928483593768847,
332 .92387953251128675612818318939678828682241662586364,
333 .98078528040323044912618223613423903697393373089333,
334 .99518472667219688624483695310947992157547486872985,
335 .99879545620517239271477160475910069444320361470461,
336 .99969881869620422011576564966617219685006108125772,
337 .99992470183914454092164649119638322435060646880221,
338 .99998117528260114265699043772856771617391725094433,
339 .99999529380957617151158012570011989955298763362218,
340 .99999882345170190992902571017152601904826792288976,
341 .99999970586288221916022821773876567711626389934930,
342 .99999992646571785114473148070738785694820115568892,
343 .99999998161642929380834691540290971450507605124278,
344 .99999999540410731289097193313960614895889430318945,
345 .99999999885102682756267330779455410840053741619428,
346 .99999999971275670684941397221864177608908945791828,
347 .99999999992818917670977509588385049026048033310951
348 };
349static double dfsintab[TRIG_TAB_SIZE]=
350 {
351 1.0000000000000000000000000000000000000000000000000,
352 .70710678118654752440084436210484903928483593768846,
353 .38268343236508977172845998403039886676134456248561,
354 .19509032201612826784828486847702224092769161775195,
355 .09801714032956060199419556388864184586113667316749,
356 .04906767432741801425495497694268265831474536302574,
357 .02454122852291228803173452945928292506546611923944,
358 .01227153828571992607940826195100321214037231959176,
359 .00613588464915447535964023459037258091705788631738,
360 .00306795676296597627014536549091984251894461021344,
361 .00153398018628476561230369715026407907995486457522,
362 .00076699031874270452693856835794857664314091945205,
363 .00038349518757139558907246168118138126339502603495,
364 .00019174759731070330743990956198900093346887403385,
365 .00009587379909597734587051721097647635118706561284,
366 .00004793689960306688454900399049465887274686668768,
367 .00002396844980841821872918657716502182009476147489,
368 .00001198422490506970642152156159698898480473197753
369 };
370
371static double dhalsec[TRIG_TAB_SIZE]=
372 {
373 0,
374 0,
375 .54119610014619698439972320536638942006107206337801,
376 .50979557910415916894193980398784391368261849190893,
377 .50241928618815570551167011928012092247859337193963,
378 .50060299823519630134550410676638239611758632599591,
379 .50015063602065098821477101271097658495974913010340,
380 .50003765191554772296778139077905492847503165398345,
381 .50000941253588775676512870469186533538523133757983,
382 .50000235310628608051401267171204408939326297376426,
383 .50000058827484117879868526730916804925780637276181,
384 .50000014706860214875463798283871198206179118093251,
385 .50000003676714377807315864400643020315103490883972,
386 .50000000919178552207366560348853455333939112569380,
387 .50000000229794635411562887767906868558991922348920,
388 .50000000057448658687873302235147272458812263401372,
389 .50000000014362164661654736863252589967935073278768,
390 .50000000003590541164769084922906986545517021050714
391 };
392
393
394#include <stdio.h>
395
396
397int writetables()
398{
399 int i;
400
401 printf("/* Tables for fixed point lookup with %d bit precision*/\n\n",fix1);
402
403 printf("static int fsintab[TRIG_TAB_SIZE]= {\n");
404
405 for (i=0;i<TRIG_TAB_SIZE-1;i++)
406 printf("%ld, \n",ftofix(dfsintab[i]));
407 printf("%ld }; \n",ftofix(dfsintab[i]));
408
409
410 printf("\n\nstatic int fcostab[TRIG_TAB_SIZE]= {\n");
411
412 for (i=0;i<TRIG_TAB_SIZE-1;i++)
413 printf("%ld, \n",ftofix(dfcostab[i]));
414 printf("%ld }; \n",ftofix(dfcostab[i]));
415
416}
417
418
419#define OUTPUT \
420 fprintf(stderr,"Integer - Float\n");\
421 for (i=0;i<NP;i++)\
422 fprintf(stderr,"%f %f --> %f %f\n",fixtof(r[i]),fixtof(im[i]),\
423 fr[i],fim[i]);\
424 fprintf(stderr,"\n\n");
425
426
427
428int main()
429{
430 int i;
431 t_fixed r[256];
432 t_fixed im[256];
433 float fr[256];
434 float fim[256];
435
436
437#if 1
438 writetables();
439 exit(0);
440#endif
441
442
443#if 0
444 t_fixed c1,s1,c2,s2,c3,s3;
445 int k;
446 int i;
447
448 TRIG_VARS;
449
450 for (k=2;k<10;k+=2) {
451 TRIG_INIT(k,c1,s1);
452 for (i=0;i<8;i++) {
453 TRIG_NEXT(k,c1,s1);
454 TRIG_23(k,c1,s1,c2,s2,c3,s3);
455 printf("TRIG value k=%d,%d val1 = %f %f\n",k,i,fixtof(c1),fixtof(s1));
456 }
457 }
458#endif
459
460
461
462#if 1
463
464 #define NP 16
465
466 for (i=0;i<NP;i++) {
467 fr[i] = 0.;
468 r[i] = 0.;
469 fim[i] = 0.;
470 im[i] = 0;
471 }
472
473#if 0
474 for (i=0;i<NP;i++) {
475 if (i&1) {
476 fr[i] = 0.1*i;
477 r[i] = ftofix(0.1*i);
478 }
479 else {
480 fr[i] = 0.;
481 r[i] = 0.;
482 }
483 }
484#endif
485#if 0
486 for (i=0;i<NP;i++) {
487 fr[i] = 0.1;
488 r[i] = ftofix(0.1);
489 }
490#endif
491
492 r[1] = ftofix(0.1);
493 fr[1] = 0.1;
494
495
496
497 /* TEST RUN */
498
499 OUTPUT
500
501 imayer_fft(NP,r,im);
502 mayer_fft(NP,fr,fim);
503// imayer_fht(r,NP);
504// mayer_fht(fr,NP);
505
506#if 1
507 for (i=0;i<NP;i++) {
508 r[i] = mult(ftofix(0.01),r[i]);
509 fr[i] = 0.01*fr[i];
510 }
511#endif
512
513 OUTPUT
514
515
516 imayer_fft(NP,r,im);
517 mayer_fft(NP,fr,fim);
518
519 OUTPUT
520
521
522#endif
523
524
525}
526
527#endif
528