src/fht.h

   1 // FHT - Fast Hartley Transform Class
   2 //
   3 // Copyright (C) 2004  Melchior FRANZ - mfranz@kde.org
   4 //
   5 // This program is free software; you can redistribute it and/or
   6 // modify it under the terms of the GNU General Public License as
   7 // published by the Free Software Foundation; either version 2 of the
   8 // License, or (at your option) any later version.
   9 //
  10 // This program is distributed in the hope that it will be useful, but
  11 // WITHOUT ANY WARRANTY; without even the implied warranty of
  12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 // General Public License for more details.
  14 //
  15 // You should have received a copy of the GNU General Public License
  16 // along with this program; if not, write to the Free Software
  17 // Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
  18 //
  19 // $Id$
  20
  21 #ifndef FHT_H
  22 #define FHT_H
  23
  24 /**
  25  * Implementation of the Hartley Transform after Bracewell's discrete
  26  * algorithm. The algorithm is subject to US patent No. 4,646,256 (1987)
  27  * but was put into public domain by the Board of Trustees of Stanford
  28  * University in 1994 and is now freely available[1].
  29  *
  30  * [1] Computer in Physics, Vol. 9, No. 4, Jul/Aug 1995 pp 373-379
  31  */
  32 class FHT
  33 {
  34         int     m_exp2;
  35         int     m_num;
  36         float   *m_buf;
  37         float   *m_tab;
  38         int     *m_log;
  39
  40         /**
  41          * Create a table of "cas" (cosine and sine) values.
  42          * Has only to be done in the constructor and saves from
  43          * calculating the same values over and over while transforming.
  44          */
  45         void    makeCasTable();
  46
  47         /**
  48          * Recursive in-place Hartley transform. For internal use only!
  49          */
  50         void    _transform(float *, int, int);
  51
  52    public:
  53         /**
  54         * Prepare transform for data sets with @f$2^n@f$ numbers, whereby @f$n@f$
  55         * should be at least 3. Values of more than 3 need a trigonometry table.
  56         * @see makeCasTable()
  57         */
  58         FHT(int);
  59
  60         ~FHT();
  61         inline int sizeExp() const { return m_exp2; }
  62         inline int size() const { return m_num; }
  63         float   *copy(float *, float *);
  64         float   *clear(float *);
  65         void    scale(float *, float);
  66
  67         /**
  68          * Exponentially Weighted Moving Average (EWMA) filter.
  69          * @param d is the filtered data.
  70          * @param s is fresh input.
  71          * @param w is the weighting factor.
  72          */
  73         void    ewma(float *d, float *s, float w);
  74
  75         /**
  76          * Logarithmic audio spectrum. Maps semi-logarithmic spectrum
  77          * to logarithmic frequency scale, interpolates missing values.
  78          * A logarithmic index map is calculated at the first run only.
  79          * @param p is the input array.
  80          * @param out is the spectrum.
  81          */
  82         void    logSpectrum(float *out, float *p);
  83
  84         /**
  85          * Semi-logarithmic audio spectrum.
  86          */
  87         void    semiLogSpectrum(float *);
  88
  89         /**
  90          * Fourier spectrum.
  91          */
  92         void    spectrum(float *);
  93
  94         /**
  95          * Calculates a mathematically correct FFT power spectrum.
  96          * If further scaling is applied later, use power2 instead
  97          * and factor the 0.5 in the final scaling factor.
  98          * @see FHT::power2()
  99          */
 100         void    power(float *);
 101
 102         /**
 103          * Calculates an FFT power spectrum with doubled values as a
 104          * result. The values need to be multiplied by 0.5 to be exact.
 105          * Note that you only get @f$2^{n-1}@f$ power values for a data set
 106          * of @f$2^n@f$ input values. This is the fastest transform.
 107          * @see FHT::power()
 108          */
 109         void    power2(float *);
 110
 111         /**
 112          * Discrete Hartley transform of data sets with 8 values.
 113          */
 114         void    transform8(float *);
 115
 116         void    transform(float *);
 117 };
 118
 119 #endif