src/calf/onepole.h

   1 /* Calf DSP Library
   2  * Basic one-pole one-zero filters based on bilinear transform.
   3  * Copyright (C) 2001-2007 Krzysztof Foltman
   4  *
   5  * This program is free software; you can redistribute it and/or
   6  * modify it under the terms of the GNU Lesser General Public
   7  * License as published by the Free Software Foundation; either
   8  * version 2 of the License, or (at your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful,
  11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * Lesser General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU Lesser General
  16  * Public License along with this program; if not, write to the
  17  * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  18  * Boston, MA 02111-1307, USA.
  19  */
  20 #ifndef __CALF_ONEPOLE_H
  21 #define __CALF_ONEPOLE_H
  22
  23 #include "primitives.h"
  24
  25 namespace dsp {
  26
  27 /**
  28  * one-pole filter, for floating point values
  29  * coefficient calculation is based on bilinear transform, and the code itself is based on my very old OneSignal lib
  30  * lp and hp are *somewhat* tested, allpass is not tested at all
  31  * don't use this for integers because it won't work
  32  */
  33 template<class T = float, class Coeff = float>
  34 class onepole
  35 {
  36 public:
  37     typedef std::complex<double> cfloat;
  38
  39     T x1, y1;
  40     Coeff a0, a1, b1;
  41
  42     onepole()
  43     {
  44         reset();
  45     }
  46
  47     /// Set coefficients for a lowpass filter
  48     void set_lp(float fc, float sr)
  49     {
  50         //   x   x
  51         //  x+1 x-1
  52         Coeff x = tan (M_PI * fc / (2 * sr));
  53         Coeff q = 1/(1+x);
  54         a0 = a1 = x*q;
  55         b1 = (x-1)*q;
  56     }
  57
  58     /// Set coefficients for an allpass filter
  59     void set_ap(float fc, float sr)
  60     {
  61         // x-1  x+1
  62         // x+1  x-1
  63         Coeff x = tan (M_PI * fc / (2 * sr));
  64         Coeff q = 1/(1+x);
  65         b1 = a0 = (x-1)*q;
  66         a1 = 1;
  67     }
  68
  69     /// Set coefficients for an allpass filter, using omega instead of fc and sr
  70     /// omega = (PI / 2) * fc / sr
  71     void set_ap_w(float w)
  72     {
  73         // x-1  x+1
  74         // x+1  x-1
  75         Coeff x = tan (w);
  76         Coeff q = 1/(1+x);
  77         b1 = a0 = (x-1)*q;
  78         a1 = 1;
  79     }
  80
  81     /// Set coefficients for a highpass filter
  82     void set_hp(float fc, float sr)
  83     {
  84         //   x   -x
  85         //  x+1  x-1
  86         Coeff x = tan (M_PI * fc / (2 * sr));
  87         Coeff q = 1/(1+x);
  88         a0 = q;
  89         a1 = -a0;
  90         b1 = (x-1)*q;
  91     }
  92
  93     /// Process one sample
  94     inline T process(T in)
  95     {
  96         T out = in * a0 + x1 * a1 - y1 * b1;
  97         x1 = in;
  98         y1 = out;
  99         return out;
 100     }
 101
 102     /// Process one sample, assuming it's a lowpass filter (optimized special case)
 103     inline T process_lp(T in)
 104     {
 105         T out = (in + x1) * a0 - y1 * b1;
 106         x1 = in;
 107         y1 = out;
 108         return out;
 109     }
 110
 111     /// Process one sample, assuming it's a highpass filter (optimized special case)
 112     inline T process_hp(T in)
 113     {
 114         T out = (in - x1) * a0 - y1 * b1;
 115         x1 = in;
 116         y1 = out;
 117         return out;
 118     }
 119
 120     /// Process one sample, assuming it's an allpass filter (optimized special case)
 121     inline T process_ap(T in)
 122     {
 123         T out = (in - y1) * a0 + x1;
 124         x1 = in;
 125         y1 = out;
 126         return out;
 127     }
 128
 129     /// Process one sample using external state variables
 130     inline T process_ap(T in, float &x1, float &y1)
 131     {
 132         T out = (in - y1) * a0 + x1;
 133         x1 = in;
 134         y1 = out;
 135         return out;
 136     }
 137
 138     /// Process one sample using external state variables, including filter coeff
 139     inline T process_ap(T in, float &x1, float &y1, float a0)
 140     {
 141         T out = (in - y1) * a0 + x1;
 142         x1 = in;
 143         y1 = out;
 144         return out;
 145     }
 146
 147     inline bool empty() const {
 148         return y1 == 0;
 149     }
 150
 151     inline void sanitize()
 152     {
 153         dsp::sanitize(x1);
 154         dsp::sanitize(y1);
 155     }
 156
 157     inline void reset()
 158     {
 159         dsp::zero(x1);
 160         dsp::zero(y1);
 161     }
 162
 163     template<class U>
 164     inline void copy_coeffs(const onepole<U> &src)
 165     {
 166         a0 = src.a0;
 167         a1 = src.a1;
 168         b1 = src.b1;
 169     }
 170
 171     /// Return the filter's gain at frequency freq
 172     /// @param freq   Frequency to look up
 173     /// @param sr     Filter sample rate (used to convert frequency to angular frequency)
 174     float freq_gain(float freq, float sr) const
 175     {
 176         freq *= 2.0 * M_PI / sr;
 177         cfloat z = 1.0 / exp(cfloat(0.0, freq));
 178
 179         return std::abs(h_z(z));
 180     }
 181
 182     /// Return H(z) the filter's gain at frequency freq
 183     /// @param z   Z variable (e^jw)
 184     cfloat h_z(const cfloat &z) const
 185     {
 186         return (cfloat(a0) + double(a1) * z) / (cfloat(1.0) + double(b1) * z);
 187     }
 188 };
 189
 190 };
 191
 192 #endif