libs/Music/Filter.cpp

   1 // Copyright 2007 "Gilles Degottex"
   2
   3 // This file is part of "Music"
   4
   5 // "Music" is free software; you can redistribute it and/or modify
   6 // it under the terms of the GNU Lesser General Public License as published by
   7 // the Free Software Foundation; either version 2.1 of the License, or
   8 // (at your option) any later version.
   9 //
  10 // "Music" 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
  13 // GNU Lesser General Public License for more details.
  14 //
  15 // You should have received a copy of the GNU Lesser General Public License
  16 // along with this program; if not, write to the Free Software
  17 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  18
  19 #include "Filter.h"
  20
  21 // #include <iostream>
  22 using namespace std;
  23 #include <CppAddons/CAMath.h>
  24 using namespace Math;
  25 #include "Music.h"
  26
  27 // b = fir1(32,[0.00001 0.23]);
  28 vector<double> Music::fir1_lowpass(int n, double cutoff)
  29 {
  30         vector<double> b(n);
  31
  32         int d = (n-1)/2;
  33
  34         for(size_t i=0; i<b.size(); i++)
  35                 b[i] = cutoff*sinc(cutoff*(i-d));
  36
  37         return b;
  38 }
  39
  40 vector<double> Music::fir1_highpass(int n, double cutoff)
  41 {
  42         vector<double> b(n);
  43
  44         int d = (n-1)/2;
  45
  46         double s=0.0;
  47         for(size_t i=0; i<b.size(); i++)
  48         {
  49                 b[i] = cutoff*sinc(cutoff*(i-d));
  50                 s += b[i];
  51         }
  52
  53         for(size_t i=0; i<b.size(); i++)
  54                 b[i] = -b[i];
  55         b[d] = s+b[d];
  56
  57         return b;
  58 }
  59
  60 vector<double> Music::fir1_bandpass(int n, double low_cutoff, double high_cutoff)
  61 {
  62         vector<double> b(n, 0.0);
  63
  64         if(low_cutoff>high_cutoff)
  65                 return b;
  66
  67         vector<double> lowf = fir1_lowpass(n, low_cutoff);
  68         vector<double> highf = fir1_highpass(n, high_cutoff);
  69
  70         return conv(lowf, highf);
  71 }
  72
  73 Music::FIRRTFilter::FIRRTFilter(std::vector<double>& imp_res)
  74 {
  75         assert(imp_res.size()>0);
  76
  77         m_imp_res = imp_res;
  78 //                      m_to_filter.reserve(imp_res.size());
  79 }
  80
  81 double Music::FIRRTFilter::operator()(double v)
  82 {
  83         double value = 0.0;
  84
  85         m_to_filter.push_front(v);
  86
  87         if(m_to_filter.size()>=m_imp_res.size())
  88         {
  89                                 // convolve
  90                 for(size_t i=0; i<m_imp_res.size(); i++)
  91                         value += m_imp_res[i]*m_to_filter[i];
  92
  93 //              if(decim++%4==0)
  94 //                      m_queue.push_front(value);
  95
  96                 // drop unused data
  97                 m_to_filter.pop_back();
  98         }
  99
 100         return value;
 101 }
 102
 103 Music::RectangularHighPassRTFilter::RectangularHighPassRTFilter(int N)
 104 {
 105         reset(N);
 106 }
 107 void Music::RectangularHighPassRTFilter::reset(int N)
 108 {
 109         if(N<1) N=1;
 110
 111         m_N = N;
 112         m_sum = 0.0;
 113         m_summed_values.clear();
 114 }
 115 double Music::RectangularHighPassRTFilter::operator()(double v)
 116 {
 117         m_summed_values.push_front(v);
 118         m_sum += v;
 119         while(int(m_summed_values.size())>m_N)
 120         {
 121                 m_sum -= m_summed_values.back();
 122
 123                 m_summed_values.pop_back();
 124         }
 125
 126         return m_summed_values[m_summed_values.size()/2] - m_sum/m_summed_values.size();
 127 }
 128
 129 double Music::RectangularLowPassRTFilter::operator()(double v)
 130 {
 131         return v;
 132 }
 133
 134 double Music::RectangularBandPassRTFilter::operator()(double v)
 135 {
 136         return v;
 137 }
 138
 139 /*
 140 LP and HP filter
 141 Type : biquad, tweaked butterworthReferences : Posted by Patrice TarrabiaCode : www.musicdsp.org
 142 r  = rez amount, from sqrt(2) to ~ 0.1
 143 f  = cutoff frequency
 144 (from ~0 Hz to SampleRate/2 - though many
 145 synths seem to filter only  up to SampleRate/4)
 146
 147 The filter algo:
 148 out(n) = a1 * in + a2 * in(n-1) + a3 * in(n-2) - b1*out(n-1) - b2*out(n-2)
 149
 150 Lowpass:
 151       c = 1.0 / tan(pi * f / sample_rate);
 152
 153       a1 = 1.0 / ( 1.0 + r * c + c * c);
 154       a2 = 2* a1;
 155       a3 = a1;
 156       b1 = 2.0 * ( 1.0 - c*c) * a1;
 157       b2 = ( 1.0 - r * c + c * c) * a1;
 158
 159 Hipass:
 160       c = tan(pi * f / sample_rate);
 161
 162       a1 = 1.0 / ( 1.0 + r * c + c * c);
 163       a2 = -2*a1;
 164       a3 = a1;
 165       b1 = 2.0 * ( c*c - 1.0) * a1;
 166       b2 = ( 1.0 - r * c + c * c) * a1;
 167 */