1 // This file is part of the FXT library.
2 // Copyright (C) 2010, 2011 Joerg Arndt
3 // License: GNU General Public License version 3 or later,
4 // see the file COPYING.txt in the main directory.
5
6 #include "fft/fft-default.h" // needs fft
7 #include "aux1/arith1.h" // multiply_val()
8 #include "fxttypes.h" // ulong
9
10 #include <cmath> // sqrt()
11
12
13
14 static const int is = +1; // sign of the (forward-) transform
15
inv_sqrt(ulong n)16 static inline double inv_sqrt(ulong n)
17 {
18 return 1.0/sqrt((double)(n));
19 }
20 // -------------------------
21
22 void
harmonic_wavelet(double * fr,double * fi,ulong ldn)23 harmonic_wavelet(double *fr, double *fi, ulong ldn)
24 // Harmonic wavelet transform.
25 // Cf. David E. Newland: "Harmonic Wavelet Analysis",
26 // Proceedings R. Soc. Lond. A,: Mathematical and Physical Sciences,
27 // vol.443, no.1917, pp.203-225, (1993).
28 {
29 ulong n = (1UL<<ldn);
30 if ( n<4 ) return;
31
32 fft(fr, fi, ldn, +is);
33 multiply_val(fr, n, inv_sqrt(n));
34 multiply_val(fi, n, inv_sqrt(n));
35
36 for (ulong ldm=1; ldm<=ldn-1; ++ldm)
37 {
38 ulong m = (1UL<<ldm); // m=2, 4, 8, ..., n/2
39
40 fft(fr+m, fi+m, ldm, +is);
41 multiply_val(fr+m, m, inv_sqrt(m));
42 multiply_val(fi+m, m, inv_sqrt(m));
43 }
44
45 for (ulong k=n/2+1; k<n; ++k) fi[k] = -fi[k];
46 }
47 // -------------------------
48
49
50 void
inverse_harmonic_wavelet(double * fr,double * fi,ulong ldn)51 inverse_harmonic_wavelet(double *fr, double *fi, ulong ldn)
52 // Inverse of harmonic_wavelet()
53 {
54 ulong n = (1UL<<ldn);
55 if ( n<4 ) return;
56
57 for (ulong k=n/2+1; k<n; ++k) fi[k] = -fi[k];
58
59 for (ulong ldm=ldn-1; ldm>=1; --ldm)
60 {
61 ulong m = (1UL<<ldm); // m=n/2, n/4, ..., 4, 2
62 fft(fr+m, fi+m, ldm, -is);
63 multiply_val(fr+m, m, inv_sqrt(m));
64 multiply_val(fi+m, m, inv_sqrt(m));
65 }
66
67 fft(fr, fi, ldn, -is);
68 multiply_val(fr, n, inv_sqrt(n));
69 multiply_val(fi, n, inv_sqrt(n));
70 }
71 // -------------------------
72
73