1 #if !defined HAVE_HARTLEYSHIFT_H__
2 #define HAVE_HARTLEYSHIFT_H__
3 // This file is part of the FXT library.
4 // Copyright (C) 2010, 2012, 2017 Joerg Arndt
5 // License: GNU General Public License version 3 or later,
6 // see the file COPYING.txt in the main directory.
7
8 #include "aux0/sincos.h" // SinCos()
9 #include "aux0/sumdiff.h" // sumdiff()
10
11 #include "fxttypes.h"
12
13 #include <cmath> // sin()
14
15 // for trig recursions:
16 #define Tdouble long double
17 #define Sin sinl
18 // use long double, else loss of precision with fht_loc_*()
19
20 template <typename Type>
hartley_shift_05_v1(Type * f,ulong n)21 inline void hartley_shift_05_v1(Type *f, ulong n)
22 // Hartley transform analogue to fourier_shift(f, n, 0.5)
23 // n := length of array
24 // Used for negacyclic convolution and recursive FHTs.
25 // Self-inverse.
26 {
27 const double phi0 = M_PI/n;
28 double phi = 0.0;
29 for (ulong i=1, j=n-1; i<j; ++i, --j)
30 {
31 phi += phi0;
32 double s, c;
33 SinCos(phi, &s, &c);
34
35 Type fi = f[i], fj = f[j];
36 f[i] = fi * c + fj * s;
37 f[j] = fi * s - fj * c;
38 }
39 }
40 // -------------------------
41
42 template <typename Type>
hartley_shift_05_v1rec(Type * f,ulong n)43 inline void hartley_shift_05_v1rec(Type *f, ulong n)
44 // Same as hartley_shift_05_v1() but with trigonometric recursion.
45 {
46 const double phi0 = M_PI/n;
47 Tdouble s = 0.0, c = 1.0;
48 Tdouble be = Sin(phi0), al = Sin(0.5*phi0); al *= (2.0*al);
49 for (ulong i=1, j=n-1; i<j; ++i, --j)
50 {
51 { Tdouble tt = c; c -= (al*tt+be*s); s -= (al*s-be*tt); }
52
53 Type fi = f[i], fj = f[j];
54 f[i] = fi * (double)c + fj * (double)s; // jjcast
55 f[j] = fi * (double)s - fj * (double)c; // jjcast
56 }
57 }
58 // -------------------------
59
60 template <typename Type>
hartley_shift_05_v2(Type * f,ulong n)61 inline void hartley_shift_05_v2(Type *f, ulong n)
62 // Optimized version, n must be a power of 2.
63 {
64 const double phi0 = M_PI/n;
65 double phi = 0.0;
66 const ulong nh = n/2;
67 if ( n>=4 )
68 {
69 ulong im=nh/2, jm=3*im;
70 Type fi = f[im], fj = f[jm];
71 double cs = M_SQRT1_2;
72 f[im] = (fi + fj) * cs;
73 f[jm] = (fi - fj) * cs;
74
75 if ( n>=8 )
76 {
77 for (ulong i=1, j=n-1, k=nh-1, l=nh+1; i<k; ++i, --j, --k, ++l)
78 {
79 phi += phi0;
80 double s, c;
81 SinCos(phi, &s, &c);
82
83 fi = f[i];
84 fj = f[j];
85 f[i] = fi * c + fj * s;
86 f[j] = fi * s - fj * c;
87
88 // l == i + nh; k == j - nh; swap2(c, s);
89 fi = f[k];
90 fj = f[l];
91 f[k] = fi * s + fj * c;
92 f[l] = fi * c - fj * s;
93 }
94 }
95 }
96 }
97 // -------------------------
98
99 template <typename Type>
hartley_shift_05_v2rec(Type * f,ulong n)100 inline void hartley_shift_05_v2rec(Type *f, ulong n)
101 // Optimized version, n must be a power of 2.
102 // Same as hartley_shift_05_v2() but with trigonometric recursion.
103 {
104 const ulong nh = n/2;
105 if ( n>=4 )
106 {
107 ulong im=nh/2, jm=3*im;
108 Type fi = f[im], fj = f[jm];
109 double cs = M_SQRT1_2;
110 f[im] = (fi + fj) * cs;
111 f[jm] = (fi - fj) * cs;
112
113 if ( n>=8 )
114 {
115 const Tdouble phi0 = M_PI/(double)n;
116 Tdouble be = Sin(phi0), al = Sin(0.5*phi0); al *= (2.0*al);
117 Tdouble s = 0.0, c = 1.0;
118 for (ulong i=1, j=n-1, k=nh-1, l=nh+1; i<k; ++i, --j, --k, ++l)
119 {
120 { Tdouble tt = c; c -= (al*tt+be*s); s -= (al*s-be*tt); }
121
122 fi = f[i];
123 fj = f[j];
124 f[i] = fi * (double)c + fj * (double)s; // jjcast
125 f[j] = fi * (double)s - fj * (double)c; // jjcast
126
127 // l = i + nh; k = j - nh;
128 fi = f[k];
129 fj = f[l];
130 f[k] = fi * (double)s + fj * (double)c; // jjcast
131 f[l] = fi * (double)c - fj * (double)s; // jjcast
132 }
133 }
134 }
135 }
136 // -------------------------
137
138
139
140
141 //: inlines give default implementation:
142 template <typename Type>
hartley_shift_05(Type * a,ulong n)143 inline void hartley_shift_05(Type *a, ulong n)
144 { hartley_shift_05_v2rec(a, n); }
145
146
147 #undef Tdouble
148 #undef Sin
149
150 #endif // !defined HAVE_HARTLEYSHIFT_H__
151