libs/qm-dsp/dsp/rateconversion/Decimator.cpp

   1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
   2
   3 /*
   4     QM DSP Library
   5
   6     Centre for Digital Music, Queen Mary, University of London.
   7     This file 2005-2006 Christian Landone.
   8
   9     This program is free software; you can redistribute it and/or
  10     modify it under the terms of the GNU General Public License as
  11     published by the Free Software Foundation; either version 2 of the
  12     License, or (at your option) any later version.  See the file
  13     COPYING included with this distribution for more information.
  14 */
  15
  16 #include "Decimator.h"
  17
  18 #include <iostream>
  19
  20 //////////////////////////////////////////////////////////////////////
  21 // Construction/Destruction
  22 //////////////////////////////////////////////////////////////////////
  23
  24 Decimator::Decimator( unsigned int inLength, unsigned int decFactor )
  25 {
  26
  27     m_inputLength = 0;
  28     m_outputLength = 0;
  29     m_decFactor = 1;
  30
  31     initialise( inLength, decFactor );
  32 }
  33
  34 Decimator::~Decimator()
  35 {
  36     deInitialise();
  37 }
  38
  39 void Decimator::initialise( unsigned int inLength, unsigned int decFactor)
  40 {
  41     m_inputLength = inLength;
  42     m_decFactor = decFactor;
  43     m_outputLength = m_inputLength / m_decFactor;
  44
  45     decBuffer = new double[ m_inputLength ];
  46
  47     // If adding new factors here, add them to
  48     // getHighestSupportedFactor in the header as well
  49
  50     if(m_decFactor == 8)
  51     {
  52         //////////////////////////////////////////////////
  53         b[0] = 0.060111378492136;
  54         b[1] = -0.257323420830598;
  55         b[2] = 0.420583503165928;
  56         b[3] = -0.222750785197418;
  57         b[4] = -0.222750785197418;
  58         b[5] = 0.420583503165928;
  59         b[6] = -0.257323420830598;
  60         b[7] = 0.060111378492136;
  61
  62         a[0] = 1;
  63         a[1] = -5.667654878577432;
  64         a[2] = 14.062452278088417;
  65         a[3] = -19.737303840697738;
  66         a[4] = 16.889698874608641;
  67         a[5] = -8.796600612325928;
  68         a[6] = 2.577553446979888;
  69         a[7] = -0.326903916815751;
  70         //////////////////////////////////////////////////
  71     }
  72     else if( m_decFactor == 4 )
  73     {
  74         //////////////////////////////////////////////////
  75         b[ 0 ] = 0.10133306904918619;
  76         b[ 1 ] = -0.2447523353702363;
  77         b[ 2 ] = 0.33622528590120965;
  78         b[ 3 ] = -0.13936581560633518;
  79         b[ 4 ] = -0.13936581560633382;
  80         b[ 5 ] = 0.3362252859012087;
  81         b[ 6 ] = -0.2447523353702358;
  82         b[ 7 ] = 0.10133306904918594;
  83
  84         a[ 0 ] = 1;
  85         a[ 1 ] = -3.9035590278139427;
  86         a[ 2 ] = 7.5299379980621133;
  87         a[ 3 ] = -8.6890803793177511;
  88         a[ 4 ] = 6.4578667096099176;
  89         a[ 5 ] = -3.0242979431223631;
  90         a[ 6 ] = 0.83043385136748382;
  91         a[ 7 ] = -0.094420800837809335;
  92         //////////////////////////////////////////////////
  93     }
  94     else if( m_decFactor == 2 )
  95     {
  96         //////////////////////////////////////////////////
  97         b[ 0 ] = 0.20898944260075727;
  98         b[ 1 ] = 0.40011234879814367;
  99         b[ 2 ] = 0.819741973072733;
 100         b[ 3 ] = 1.0087419911682323;
 101         b[ 4 ] = 1.0087419911682325;
 102         b[ 5 ] = 0.81974197307273156;
 103         b[ 6 ] = 0.40011234879814295;
 104         b[ 7 ] = 0.20898944260075661;
 105
 106         a[ 0 ] = 1;
 107         a[ 1 ] = 0.0077331184208358217;
 108         a[ 2 ] = 1.9853971155964376;
 109         a[ 3 ] = 0.19296739275341004;
 110         a[ 4 ] = 1.2330748872852182;
 111         a[ 5 ] = 0.18705341389316466;
 112         a[ 6 ] = 0.23659265908013868;
 113         a[ 7 ] = 0.032352924250533946;
 114     }
 115     else
 116     {
 117         if ( m_decFactor != 1 ) {
 118             std::cerr << "WARNING: Decimator::initialise: unsupported decimation factor " << m_decFactor << ", no antialiasing filter will be used" << std::endl;
 119         }
 120
 121         //////////////////////////////////////////////////
 122         b[ 0 ] = 1;
 123         b[ 1 ] = 0;
 124         b[ 2 ] = 0;
 125         b[ 3 ] = 0;
 126         b[ 4 ] = 0;
 127         b[ 5 ] = 0;
 128         b[ 6 ] = 0;
 129         b[ 7 ] = 0;
 130
 131         a[ 0 ] = 1;
 132         a[ 1 ] = 0;
 133         a[ 2 ] = 0;
 134         a[ 3 ] = 0;
 135         a[ 4 ] = 0;
 136         a[ 5 ] = 0;
 137         a[ 6 ] = 0;
 138         a[ 7 ] = 0;
 139     }
 140
 141     resetFilter();
 142 }
 143
 144 void Decimator::deInitialise()
 145 {
 146     delete [] decBuffer;
 147 }
 148
 149 void Decimator::resetFilter()
 150 {
 151     Input = Output = 0;
 152
 153     o1=o2=o3=o4=o5=o6=o7=0;
 154 }
 155
 156 void Decimator::doAntiAlias(const double *src, double *dst, unsigned int length)
 157 {
 158
 159     for( unsigned int i = 0; i < length; i++ )
 160     {
 161         Input = (double)src[ i ];
 162
 163         Output = Input * b[ 0 ] + o1;
 164
 165         o1 = Input * b[ 1 ] - Output * a[ 1 ] + o2;
 166         o2 = Input * b[ 2 ] - Output * a[ 2 ] + o3;
 167         o3 = Input * b[ 3 ] - Output * a[ 3 ] + o4;
 168         o4 = Input * b[ 4 ] - Output * a[ 4 ] + o5;
 169         o5 = Input * b[ 5 ] - Output * a[ 5 ] + o6;
 170         o6 = Input * b[ 6 ] - Output * a[ 6 ] + o7;
 171         o7 = Input * b[ 7 ] - Output * a[ 7 ] ;
 172
 173         dst[ i ] = Output;
 174     }
 175
 176 }
 177
 178 void Decimator::doAntiAlias(const float *src, double *dst, unsigned int length)
 179 {
 180
 181     for( unsigned int i = 0; i < length; i++ )
 182     {
 183         Input = (double)src[ i ];
 184
 185         Output = Input * b[ 0 ] + o1;
 186
 187         o1 = Input * b[ 1 ] - Output * a[ 1 ] + o2;
 188         o2 = Input * b[ 2 ] - Output * a[ 2 ] + o3;
 189         o3 = Input * b[ 3 ] - Output * a[ 3 ] + o4;
 190         o4 = Input * b[ 4 ] - Output * a[ 4 ] + o5;
 191         o5 = Input * b[ 5 ] - Output * a[ 5 ] + o6;
 192         o6 = Input * b[ 6 ] - Output * a[ 6 ] + o7;
 193         o7 = Input * b[ 7 ] - Output * a[ 7 ] ;
 194
 195         dst[ i ] = Output;
 196     }
 197
 198 }
 199
 200 void Decimator::process(const double *src, double *dst)
 201 {
 202     if( m_decFactor != 1 )
 203     {
 204         doAntiAlias( src, decBuffer, m_inputLength );
 205     }
 206     unsigned idx = 0;
 207
 208     for( unsigned int i = 0; i < m_outputLength; i++ )
 209     {
 210         dst[ idx++ ] = decBuffer[ m_decFactor * i ];
 211     }
 212 }
 213
 214 void Decimator::process(const float *src, float *dst)
 215 {
 216     if( m_decFactor != 1 )
 217     {
 218         doAntiAlias( src, decBuffer, m_inputLength );
 219     }
 220     unsigned idx = 0;
 221
 222     for( unsigned int i = 0; i < m_outputLength; i++ )
 223     {
 224         dst[ idx++ ] = decBuffer[ m_decFactor * i ];
 225     }
 226 }