File indexing completed on 2024-05-12 04:50:07
0001 /* 0002 FHT - Fast Hartley Transform Class 0003 0004 SPDX-FileCopyrightText: 2004 Melchior FRANZ <mfranz@kde.org> 0005 0006 SPDX-License-Identifier: GPL-2.0-or-later 0007 */ 0008 0009 #include "fht.h" 0010 0011 #include <math.h> 0012 #include <string.h> 0013 0014 FHT::FHT(int n) 0015 : m_buf(nullptr) 0016 , m_tab(nullptr) 0017 , m_log(nullptr) 0018 { 0019 if (n < 3) { 0020 m_num = 0; 0021 m_exp2 = -1; 0022 return; 0023 } 0024 m_exp2 = n; 0025 m_num = 1 << n; 0026 if (n > 3) { 0027 m_buf = new float[m_num]; 0028 m_tab = new float[m_num * 2]; 0029 makeCasTable(); 0030 } 0031 } 0032 0033 FHT::~FHT() 0034 { 0035 delete[] m_buf; 0036 delete[] m_tab; 0037 delete[] m_log; 0038 } 0039 0040 void FHT::makeCasTable(void) 0041 { 0042 float d, *costab, *sintab; 0043 int ul, ndiv2 = m_num / 2; 0044 0045 for (costab = m_tab, sintab = m_tab + m_num / 2 + 1, ul = 0; ul < m_num; ul++) { 0046 d = M_PI * ul / ndiv2; 0047 *costab = *sintab = cos(d); 0048 0049 costab += 2, sintab += 2; 0050 if (sintab > m_tab + m_num * 2) 0051 sintab = m_tab + 1; 0052 } 0053 } 0054 0055 float *FHT::copy(float *d, float *s) 0056 { 0057 return (float *)memcpy(d, s, m_num * sizeof(float)); 0058 } 0059 0060 float *FHT::clear(float *d) 0061 { 0062 return (float *)memset(d, 0, m_num * sizeof(float)); 0063 } 0064 0065 void FHT::scale(float *p, float d) 0066 { 0067 for (int i = 0; i < (m_num / 2); i++) 0068 *p++ *= d; 0069 } 0070 0071 void FHT::ewma(float *d, float *s, float w) 0072 { 0073 for (int i = 0; i < (m_num / 2); i++, d++, s++) 0074 *d = *d * w + *s * (1 - w); 0075 } 0076 0077 void FHT::logSpectrum(float *out, float *p) 0078 { 0079 int n = m_num / 2, i, j, k, *r; 0080 if (!m_log) { 0081 m_log = new int[n]; 0082 float f = n / log10((double)n); 0083 for (i = 0, r = m_log; i < n; i++, r++) { 0084 j = int(rint(log10(i + 1.0) * f)); 0085 *r = j >= n ? n - 1 : j; 0086 } 0087 } 0088 semiLogSpectrum(p); 0089 *out++ = *p = *p / 100; 0090 for (k = i = 1, r = m_log; i < n; ++i) { 0091 j = *r++; 0092 if (i == j) 0093 *out++ = p[i]; 0094 else { 0095 float base = p[k - 1]; 0096 float step = (p[j] - base) / (j - (k - 1)); 0097 for (float corr = 0; k <= j; k++, corr += step) 0098 *out++ = base + corr; 0099 } 0100 } 0101 } 0102 0103 void FHT::semiLogSpectrum(float *p) 0104 { 0105 float e; 0106 power2(p); 0107 for (int i = 0; i < (m_num / 2); i++, p++) { 0108 e = 10.0 * log10(sqrt(*p * .5)); 0109 *p = e < 0 ? 0 : e; 0110 } 0111 } 0112 0113 void FHT::spectrum(float *p) 0114 { 0115 power2(p); 0116 for (int i = 0; i < (m_num / 2); i++, p++) 0117 *p = (float)sqrt(*p * .5); 0118 } 0119 0120 void FHT::power(float *p) 0121 { 0122 power2(p); 0123 for (int i = 0; i < (m_num / 2); i++) 0124 *p++ *= .5; 0125 } 0126 0127 void FHT::power2(float *p) 0128 { 0129 int i; 0130 float *q; 0131 _transform(p, m_num, 0); 0132 0133 *p = (*p * *p), *p += *p, p++; 0134 0135 for (i = 1, q = p + m_num - 2; i < (m_num / 2); i++, --q) 0136 *p = (*p * *p) + (*q * *q), p++; 0137 } 0138 0139 void FHT::transform(float *p) 0140 { 0141 if (m_num == 8) 0142 transform8(p); 0143 else 0144 _transform(p, m_num, 0); 0145 } 0146 0147 void FHT::transform8(float *p) 0148 { 0149 float a, b, c, d, e, f, g, h, b_f2, d_h2; 0150 float a_c_eg, a_ce_g, ac_e_g, aceg, b_df_h, bdfh; 0151 0152 a = *p++, b = *p++, c = *p++, d = *p++; 0153 e = *p++, f = *p++, g = *p++, h = *p; 0154 b_f2 = (b - f) * M_SQRT2; 0155 d_h2 = (d - h) * M_SQRT2; 0156 0157 a_c_eg = a - c - e + g; 0158 a_ce_g = a - c + e - g; 0159 ac_e_g = a + c - e - g; 0160 aceg = a + c + e + g; 0161 0162 b_df_h = b - d + f - h; 0163 bdfh = b + d + f + h; 0164 0165 *p = a_c_eg - d_h2; 0166 *--p = a_ce_g - b_df_h; 0167 *--p = ac_e_g - b_f2; 0168 *--p = aceg - bdfh; 0169 *--p = a_c_eg + d_h2; 0170 *--p = a_ce_g + b_df_h; 0171 *--p = ac_e_g + b_f2; 0172 *--p = aceg + bdfh; 0173 } 0174 0175 void FHT::_transform(float *p, int n, int k) 0176 { 0177 if (n == 8) { 0178 transform8(p + k); 0179 return; 0180 } 0181 0182 int i, j, ndiv2 = n / 2; 0183 float a, *t1, *t2, *t3, *t4, *ptab, *pp; 0184 0185 for (i = 0, t1 = m_buf, t2 = m_buf + ndiv2, pp = &p[k]; i < ndiv2; i++) 0186 *t1++ = *pp++, *t2++ = *pp++; 0187 0188 memcpy(p + k, m_buf, sizeof(float) * n); 0189 0190 _transform(p, ndiv2, k); 0191 _transform(p, ndiv2, k + ndiv2); 0192 0193 j = m_num / ndiv2 - 1; 0194 t1 = m_buf; 0195 t2 = t1 + ndiv2; 0196 t3 = p + k + ndiv2; 0197 ptab = m_tab; 0198 pp = p + k; 0199 0200 a = *ptab++ * *t3++; 0201 a += *ptab * *pp; 0202 ptab += j; 0203 0204 *t1++ = *pp + a; 0205 *t2++ = *pp++ - a; 0206 0207 for (i = 1, t4 = p + k + n; i < ndiv2; i++, ptab += j) { 0208 a = *ptab++ * *t3++; 0209 a += *ptab * *--t4; 0210 0211 *t1++ = *pp + a; 0212 *t2++ = *pp++ - a; 0213 } 0214 memcpy(p + k, m_buf, sizeof(float) * n); 0215 }