1 /* This file is part of Clementine.
2 Copyright 2004, Melchior FRANZ <mfranz@kde.org>
3 Copyright 2010, 2014, John Maguire <john.maguire@gmail.com>
4 Copyright 2014, Krzysztof Sobiecki <sobkas@gmail.com>
5
6 Clementine is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 Clementine is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with Clementine. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20 /* Original Author: Melchior FRANZ <mfranz@kde.org> 2004
21 */
22
23 #include <cmath>
24 #include "fht.h"
25
FHT(int n)26 FHT::FHT(int n) : num_((n < 3) ? 0 : 1 << n), exp2_((n < 3) ? -1 : n) {
27 if (n > 3) {
28 buf_vector_.resize(num_);
29 tab_vector_.resize(num_ * 2);
30 makeCasTable();
31 }
32 }
33
~FHT()34 FHT::~FHT() {}
35
sizeExp() const36 int FHT::sizeExp() const { return exp2_; }
size() const37 int FHT::size() const { return num_; }
38
buf_()39 float* FHT::buf_() { return buf_vector_.data(); }
tab_()40 float* FHT::tab_() { return tab_vector_.data(); }
log_()41 int* FHT::log_() { return log_vector_.data(); }
42
makeCasTable(void)43 void FHT::makeCasTable(void) {
44 float* costab = tab_();
45 float* sintab = tab_() + num_ / 2 + 1;
46
47 for (int ul = 0; ul < num_; ul++) {
48 float d = M_PI * ul / (num_ / 2);
49 *costab = *sintab = cos(d);
50
51 costab += 2;
52 sintab += 2;
53 if (sintab > tab_() + num_ * 2) sintab = tab_() + 1;
54 }
55 }
56
scale(float * p,float d)57 void FHT::scale(float* p, float d) {
58 for (int i = 0; i < (num_ / 2); i++) *p++ *= d;
59 }
60
ewma(float * d,float * s,float w)61 void FHT::ewma(float* d, float* s, float w) {
62 for (int i = 0; i < (num_ / 2); i++, d++, s++) *d = *d * w + *s * (1 - w);
63 }
64
logSpectrum(float * out,float * p)65 void FHT::logSpectrum(float* out, float* p) {
66 int n = num_ / 2, i, j, k, *r;
67 if (log_vector_.size() < n) {
68 log_vector_.resize(n);
69 float f = n / log10(static_cast<double>(n));
70 for (i = 0, r = log_(); i < n; i++, r++) {
71 j = static_cast<int>(rint(log10(i + 1.0) * f));
72 *r = j >= n ? n - 1 : j;
73 }
74 }
75 semiLogSpectrum(p);
76 *out++ = *p = *p / 100;
77 for (k = i = 1, r = log_(); i < n; i++) {
78 j = *r++;
79 if (i == j) {
80 *out++ = p[i];
81 } else {
82 float base = p[k - 1];
83 float step = (p[j] - base) / (j - (k - 1));
84 for (float corr = 0; k <= j; k++, corr += step) *out++ = base + corr;
85 }
86 }
87 }
88
semiLogSpectrum(float * p)89 void FHT::semiLogSpectrum(float* p) {
90 power2(p);
91 for (int i = 0; i < (num_ / 2); i++, p++) {
92 float e = 10.0 * log10(sqrt(*p / 2));
93 *p = e < 0 ? 0 : e;
94 }
95 }
96
spectrum(float * p)97 void FHT::spectrum(float* p) {
98 power2(p);
99 for (int i = 0; i < (num_ / 2); i++, p++)
100 *p = static_cast<float>(sqrt(*p / 2));
101 }
102
power(float * p)103 void FHT::power(float* p) {
104 power2(p);
105 for (int i = 0; i < (num_ / 2); i++) *p++ /= 2;
106 }
107
power2(float * p)108 void FHT::power2(float* p) {
109 _transform(p, num_, 0);
110
111 *p = static_cast<float>(2 * pow(*p, 2));
112 p++;
113
114 float* q = p + num_ - 2;
115 for (int i = 1; i < (num_ / 2); i++) {
116 *p = static_cast<float>(pow(*p, 2) + pow(*q, 2));
117 p++;
118 q--;
119 }
120 }
121
transform(float * p)122 void FHT::transform(float* p) {
123 if (num_ == 8)
124 transform8(p);
125 else
126 _transform(p, num_, 0);
127 }
128
transform8(float * p)129 void FHT::transform8(float* p) {
130 float a, b, c, d, e, f, g, h, b_f2, d_h2;
131 float a_c_eg, a_ce_g, ac_e_g, aceg, b_df_h, bdfh;
132
133 a = *p++, b = *p++, c = *p++, d = *p++;
134 e = *p++, f = *p++, g = *p++, h = *p;
135 b_f2 = (b - f) * M_SQRT2;
136 d_h2 = (d - h) * M_SQRT2;
137
138 a_c_eg = a - c - e + g;
139 a_ce_g = a - c + e - g;
140 ac_e_g = a + c - e - g;
141 aceg = a + c + e + g;
142
143 b_df_h = b - d + f - h;
144 bdfh = b + d + f + h;
145
146 *p = a_c_eg - d_h2;
147 *--p = a_ce_g - b_df_h;
148 *--p = ac_e_g - b_f2;
149 *--p = aceg - bdfh;
150 *--p = a_c_eg + d_h2;
151 *--p = a_ce_g + b_df_h;
152 *--p = ac_e_g + b_f2;
153 *--p = aceg + bdfh;
154 }
155
_transform(float * p,int n,int k)156 void FHT::_transform(float* p, int n, int k) {
157 if (n == 8) {
158 transform8(p + k);
159 return;
160 }
161
162 int i, j, ndiv2 = n / 2;
163 float a, *t1, *t2, *t3, *t4, *ptab, *pp;
164
165 for (i = 0, t1 = buf_(), t2 = buf_() + ndiv2, pp = &p[k]; i < ndiv2; i++)
166 *t1++ = *pp++, *t2++ = *pp++;
167
168 std::copy(buf_(), buf_() + n, p + k);
169
170 _transform(p, ndiv2, k);
171 _transform(p, ndiv2, k + ndiv2);
172
173 j = num_ / ndiv2 - 1;
174 t1 = buf_();
175 t2 = t1 + ndiv2;
176 t3 = p + k + ndiv2;
177 ptab = tab_();
178 pp = p + k;
179
180 a = *ptab++ * *t3++;
181 a += *ptab * *pp;
182 ptab += j;
183
184 *t1++ = *pp + a;
185 *t2++ = *pp++ - a;
186
187 for (i = 1, t4 = p + k + n; i < ndiv2; i++, ptab += j) {
188 a = *ptab++ * *t3++;
189 a += *ptab * *--t4;
190
191 *t1++ = *pp + a;
192 *t2++ = *pp++ - a;
193 }
194
195 std::copy(buf_(), buf_() + n, p + k);
196 }
197