nivedita velagaleti / 3DAudio

ese 519

Dependents:   PROJECT_3D_AUDIO

fhtdit2.cc@0:2076b4d80327, 2015-04-07 (annotated)

Committer:
niv17
Date:
Tue Apr 07 21:09:22 2015 +0000
Revision:
0:2076b4d80327
sonic initial april 7

Who changed what in which revision?

UserRevisionLine numberNew contents of line
niv17 0:2076b4d80327 1 // This file is part of the FXT library.
niv17 0:2076b4d80327 2 // Copyright (C) 2010, 2012 Joerg Arndt
niv17 0:2076b4d80327 3 // License: GNU General Public License version 3 or later,
niv17 0:2076b4d80327 4 // see the file COPYING.txt in the main directory.
niv17 0:2076b4d80327 5
niv17 0:2076b4d80327 6 #include "fxttypes.h"
niv17 0:2076b4d80327 7 #include "complextype.h"
niv17 0:2076b4d80327 8 #include "sincos.h"
niv17 0:2076b4d80327 9 #include "revbinpermute.h"
niv17 0:2076b4d80327 10 #include "sumdiff.h"
niv17 0:2076b4d80327 11
niv17 0:2076b4d80327 12 #include <cmath> // M_PI
niv17 0:2076b4d80327 13
niv17 0:2076b4d80327 14
niv17 0:2076b4d80327 15 void
niv17 0:2076b4d80327 16 fht_depth_first_dit2(double *f, ulong ldn)
niv17 0:2076b4d80327 17 // Radix-2 decimation in time (DIT) FHT.
niv17 0:2076b4d80327 18 // Depth-first version.
niv17 0:2076b4d80327 19 // Compared to usual fht this one
niv17 0:2076b4d80327 20 // - does more trig computations
niv17 0:2076b4d80327 21 // - is (far) better memory local
niv17 0:2076b4d80327 22 {
niv17 0:2076b4d80327 23 const ulong n = 1UL<<ldn;
niv17 0:2076b4d80327 24 revbin_permute(f, n);
niv17 0:2076b4d80327 25
niv17 0:2076b4d80327 26 for (ulong ldm=1; ldm<=ldn; ++ldm)
niv17 0:2076b4d80327 27 {
niv17 0:2076b4d80327 28 const ulong m = (1UL<<ldm);
niv17 0:2076b4d80327 29 const ulong mh = (m>>1);
niv17 0:2076b4d80327 30 const ulong m4 = (mh>>1);
niv17 0:2076b4d80327 31 const double phi0 = M_PI/(double)mh;
niv17 0:2076b4d80327 32
niv17 0:2076b4d80327 33 for (ulong r=0; r<n; r+=m)
niv17 0:2076b4d80327 34 {
niv17 0:2076b4d80327 35 { // j == 0:
niv17 0:2076b4d80327 36 ulong t1 = r;
niv17 0:2076b4d80327 37 ulong t2 = t1 + mh;
niv17 0:2076b4d80327 38 sumdiff(f[t1], f[t2]);
niv17 0:2076b4d80327 39 }
niv17 0:2076b4d80327 40
niv17 0:2076b4d80327 41 if ( m4 )
niv17 0:2076b4d80327 42 {
niv17 0:2076b4d80327 43 ulong t1 = r + m4;
niv17 0:2076b4d80327 44 ulong t2 = t1 + mh;
niv17 0:2076b4d80327 45 sumdiff(f[t1], f[t2]);
niv17 0:2076b4d80327 46 }
niv17 0:2076b4d80327 47
niv17 0:2076b4d80327 48 for (ulong j=1, k=mh-1; j<k; ++j, --k)
niv17 0:2076b4d80327 49 {
niv17 0:2076b4d80327 50 double s, c;
niv17 0:2076b4d80327 51 SinCos(phi0*(double)j, &s, &c);
niv17 0:2076b4d80327 52
niv17 0:2076b4d80327 53 ulong tj = r + mh + j;
niv17 0:2076b4d80327 54 ulong tk = r + mh + k;
niv17 0:2076b4d80327 55 double fj = f[tj];
niv17 0:2076b4d80327 56 double fk = f[tk];
niv17 0:2076b4d80327 57 f[tj] = fj * c + fk * s;
niv17 0:2076b4d80327 58 f[tk] = fj * s - fk * c;
niv17 0:2076b4d80327 59
niv17 0:2076b4d80327 60 ulong t1 = r + j;
niv17 0:2076b4d80327 61 ulong t2 = tj; // == t1 + mh;
niv17 0:2076b4d80327 62 sumdiff(f[t1], f[t2]);
niv17 0:2076b4d80327 63
niv17 0:2076b4d80327 64 t1 = r + k;
niv17 0:2076b4d80327 65 t2 = tk; // == t1 + mh;
niv17 0:2076b4d80327 66 sumdiff(f[t1], f[t2]);
niv17 0:2076b4d80327 67 }
niv17 0:2076b4d80327 68 }
niv17 0:2076b4d80327 69 }
niv17 0:2076b4d80327 70 }
niv17 0:2076b4d80327 71 // -------------------------
niv17 0:2076b4d80327 72
niv17 0:2076b4d80327 73
niv17 0:2076b4d80327 74 void
niv17 0:2076b4d80327 75 fht_dit2(double *f, ulong ldn)
niv17 0:2076b4d80327 76 // Radix-2 decimation in time (DIT) FHT.
niv17 0:2076b4d80327 77 {
niv17 0:2076b4d80327 78 const ulong n = 1UL<<ldn;
niv17 0:2076b4d80327 79 revbin_permute(f, n);
niv17 0:2076b4d80327 80
niv17 0:2076b4d80327 81 for (ulong ldm=1; ldm<=ldn; ++ldm)
niv17 0:2076b4d80327 82 {
niv17 0:2076b4d80327 83 const ulong m = (1UL<<ldm);
niv17 0:2076b4d80327 84 const ulong mh = (m>>1);
niv17 0:2076b4d80327 85 const ulong m4 = (mh>>1);
niv17 0:2076b4d80327 86 const double phi0 = M_PI/(double)mh;
niv17 0:2076b4d80327 87
niv17 0:2076b4d80327 88 for (ulong r=0; r<n; r+=m)
niv17 0:2076b4d80327 89 {
niv17 0:2076b4d80327 90 { // j == 0:
niv17 0:2076b4d80327 91 ulong t1 = r;
niv17 0:2076b4d80327 92 ulong t2 = t1 + mh;
niv17 0:2076b4d80327 93 sumdiff(f[t1], f[t2]);
niv17 0:2076b4d80327 94 }
niv17 0:2076b4d80327 95
niv17 0:2076b4d80327 96 if ( m4 )
niv17 0:2076b4d80327 97 {
niv17 0:2076b4d80327 98 ulong t1 = r + m4;
niv17 0:2076b4d80327 99 ulong t2 = t1 + mh;
niv17 0:2076b4d80327 100 sumdiff(f[t1], f[t2]);
niv17 0:2076b4d80327 101 }
niv17 0:2076b4d80327 102 }
niv17 0:2076b4d80327 103
niv17 0:2076b4d80327 104 for (ulong j=1, k=mh-1; j<k; ++j, --k)
niv17 0:2076b4d80327 105 {
niv17 0:2076b4d80327 106 double s, c;
niv17 0:2076b4d80327 107 SinCos(phi0*(double)j, &s, &c);
niv17 0:2076b4d80327 108
niv17 0:2076b4d80327 109 for (ulong r=0; r<n; r+=m)
niv17 0:2076b4d80327 110 {
niv17 0:2076b4d80327 111 ulong tj = r + mh + j;
niv17 0:2076b4d80327 112 ulong tk = r + mh + k;
niv17 0:2076b4d80327 113 double fj = f[tj];
niv17 0:2076b4d80327 114 double fk = f[tk];
niv17 0:2076b4d80327 115 f[tj] = fj * c + fk * s;
niv17 0:2076b4d80327 116 f[tk] = fj * s - fk * c;
niv17 0:2076b4d80327 117
niv17 0:2076b4d80327 118 ulong t1 = r + j;
niv17 0:2076b4d80327 119 ulong t2 = tj; // == t1 + mh;
niv17 0:2076b4d80327 120 sumdiff(f[t1], f[t2]);
niv17 0:2076b4d80327 121
niv17 0:2076b4d80327 122 t1 = r + k;
niv17 0:2076b4d80327 123 t2 = tk; // == t1 + mh;
niv17 0:2076b4d80327 124 sumdiff(f[t1], f[t2]);
niv17 0:2076b4d80327 125 }
niv17 0:2076b4d80327 126 }
niv17 0:2076b4d80327 127 }
niv17 0:2076b4d80327 128 }
niv17 0:2076b4d80327 129 // -------------------------
niv17 0:2076b4d80327 130