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