File indexing completed on 2024-05-12 05:46:32
0001 // FHT - Fast Hartley Transform Class 0002 // 0003 // Copyright (C) 2004 Melchior FRANZ - mfranz@kde.org 0004 // 0005 // This program is free software; you can redistribute it and/or 0006 // modify it under the terms of the GNU General Public License as 0007 // published by the Free Software Foundation; either version 2 of the 0008 // License, or (at your option) any later version. 0009 // 0010 // This program is distributed in the hope that it will be useful, but 0011 // WITHOUT ANY WARRANTY; without even the implied warranty of 0012 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0013 // General Public License for more details. 0014 // 0015 // You should have received a copy of the GNU General Public License 0016 // along with this program; if not, write to the Free Software 0017 // Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA 0018 // 0019 // $Id: fht.cpp 871912 2008-10-16 02:27:10Z mitchell $ 0020 0021 0022 #include "fht.h" 0023 0024 #include <math.h> 0025 #include <string.h> 0026 0027 FHT::FHT(int n) : 0028 m_buf(nullptr), 0029 m_tab(nullptr), 0030 m_log(nullptr) 0031 { 0032 if (n < 3) { 0033 m_num = 0; 0034 m_exp2 = -1; 0035 return; 0036 } 0037 m_exp2 = n; 0038 m_num = 1 << n; 0039 if (n > 3) { 0040 m_buf = new float[m_num]; 0041 m_tab = new float[m_num * 2]; 0042 makeCasTable(); 0043 } 0044 } 0045 0046 0047 FHT::~FHT() 0048 { 0049 delete[] m_buf; 0050 delete[] m_tab; 0051 delete[] m_log; 0052 } 0053 0054 0055 void FHT::makeCasTable(void) 0056 { 0057 float d, *costab, *sintab; 0058 int ul, ndiv2 = m_num / 2; 0059 0060 for (costab = m_tab, sintab = m_tab + m_num / 2 + 1, ul = 0; ul < m_num; ul++) { 0061 d = M_PI * ul / ndiv2; 0062 *costab = *sintab = cos(d); 0063 0064 costab += 2, sintab += 2; 0065 if (sintab > m_tab + m_num * 2) 0066 sintab = m_tab + 1; 0067 } 0068 } 0069 0070 0071 float* FHT::copy(float *d, float *s) 0072 { 0073 return (float *)memcpy(d, s, m_num * sizeof(float)); 0074 } 0075 0076 0077 float* FHT::clear(float *d) 0078 { 0079 return (float *)memset(d, 0, m_num * sizeof(float)); 0080 } 0081 0082 0083 void FHT::scale(float *p, float d) 0084 { 0085 for (int i = 0; i < (m_num / 2); i++) 0086 *p++ *= d; 0087 } 0088 0089 0090 void FHT::ewma(float *d, float *s, float w) 0091 { 0092 for (int i = 0; i < (m_num / 2); i++, d++, s++) 0093 *d = *d * w + *s * (1 - w); 0094 } 0095 0096 0097 void FHT::logSpectrum(float *out, float *p) 0098 { 0099 int n = m_num / 2, i, j, k, *r; 0100 if (!m_log) { 0101 m_log = new int[n]; 0102 float f = n / log10((double)n); 0103 for (i = 0, r = m_log; i < n; i++, r++) { 0104 j = int(rint(log10(i + 1.0) * f)); 0105 *r = j >= n ? n - 1 : j; 0106 } 0107 } 0108 semiLogSpectrum(p); 0109 *out++ = *p = *p / 100; 0110 for (k = i = 1, r = m_log; i < n; ++i) { 0111 j = *r++; 0112 if (i == j) 0113 *out++ = p[i]; 0114 else { 0115 float base = p[k - 1]; 0116 float step = (p[j] - base) / (j - (k - 1)); 0117 for (float corr = 0; k <= j; k++, corr += step) 0118 *out++ = base + corr; 0119 } 0120 } 0121 } 0122 0123 0124 void FHT::semiLogSpectrum(float *p) 0125 { 0126 float e; 0127 power2(p); 0128 for (int i = 0; i < (m_num / 2); i++, p++) { 0129 e = 10.0 * log10(sqrt(*p * .5)); 0130 *p = e < 0 ? 0 : e; 0131 } 0132 } 0133 0134 0135 void FHT::spectrum(float *p) 0136 { 0137 power2(p); 0138 for (int i = 0; i < (m_num / 2); i++, p++) 0139 *p = (float)sqrt(*p * .5); 0140 } 0141 0142 0143 void FHT::power(float *p) 0144 { 0145 power2(p); 0146 for (int i = 0; i < (m_num / 2); i++) 0147 *p++ *= .5; 0148 } 0149 0150 0151 void FHT::power2(float *p) 0152 { 0153 int i; 0154 float *q; 0155 _transform(p, m_num, 0); 0156 0157 *p = (*p * *p), *p += *p, p++; 0158 0159 for (i = 1, q = p + m_num - 2; i < (m_num / 2); i++, --q) 0160 *p = (*p * *p) + (*q * *q), p++; 0161 } 0162 0163 0164 void FHT::transform(float *p) 0165 { 0166 if (m_num == 8) 0167 transform8(p); 0168 else 0169 _transform(p, m_num, 0); 0170 } 0171 0172 0173 void FHT::transform8(float *p) 0174 { 0175 float a, b, c, d, e, f, g, h, b_f2, d_h2; 0176 float a_c_eg, a_ce_g, ac_e_g, aceg, b_df_h, bdfh; 0177 0178 a = *p++, b = *p++, c = *p++, d = *p++; 0179 e = *p++, f = *p++, g = *p++, h = *p; 0180 b_f2 = (b - f) * M_SQRT2; 0181 d_h2 = (d - h) * M_SQRT2; 0182 0183 a_c_eg = a - c - e + g; 0184 a_ce_g = a - c + e - g; 0185 ac_e_g = a + c - e - g; 0186 aceg = a + c + e + g; 0187 0188 b_df_h = b - d + f - h; 0189 bdfh = b + d + f + h; 0190 0191 *p = a_c_eg - d_h2; 0192 *--p = a_ce_g - b_df_h; 0193 *--p = ac_e_g - b_f2; 0194 *--p = aceg - bdfh; 0195 *--p = a_c_eg + d_h2; 0196 *--p = a_ce_g + b_df_h; 0197 *--p = ac_e_g + b_f2; 0198 *--p = aceg + bdfh; 0199 } 0200 0201 0202 void FHT::_transform(float *p, int n, int k) 0203 { 0204 if (n == 8) { 0205 transform8(p + k); 0206 return; 0207 } 0208 0209 int i, j, ndiv2 = n / 2; 0210 float a, *t1, *t2, *t3, *t4, *ptab, *pp; 0211 0212 for (i = 0, t1 = m_buf, t2 = m_buf + ndiv2, pp = &p[k]; i < ndiv2; i++) 0213 *t1++ = *pp++, *t2++ = *pp++; 0214 0215 memcpy(p + k, m_buf, sizeof(float) * n); 0216 0217 _transform(p, ndiv2, k); 0218 _transform(p, ndiv2, k + ndiv2); 0219 0220 j = m_num / ndiv2 - 1; 0221 t1 = m_buf; 0222 t2 = t1 + ndiv2; 0223 t3 = p + k + ndiv2; 0224 ptab = m_tab; 0225 pp = p + k; 0226 0227 a = *ptab++ * *t3++; 0228 a += *ptab * *pp; 0229 ptab += j; 0230 0231 *t1++ = *pp + a; 0232 *t2++ = *pp++ - a; 0233 0234 for (i = 1, t4 = p + k + n; i < ndiv2; i++, ptab += j) { 0235 a = *ptab++ * *t3++; 0236 a += *ptab * *--t4; 0237 0238 *t1++ = *pp + a; 0239 *t2++ = *pp++ - a; 0240 } 0241 memcpy(p + k, m_buf, sizeof(float) * n); 0242 }