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 
Decimator(int inLength,int decFactor)24 Decimator::Decimator( int inLength, int decFactor )
25 {
26     m_inputLength = 0;
27     m_outputLength = 0;
28     m_decFactor = 1;
29 
30     initialise( inLength, decFactor );
31 }
32 
~Decimator()33 Decimator::~Decimator()
34 {
35     deInitialise();
36 }
37 
initialise(int inLength,int decFactor)38 void Decimator::initialise( int inLength, int decFactor)
39 {
40     m_inputLength = inLength;
41     m_decFactor = decFactor;
42     m_outputLength = m_inputLength / m_decFactor;
43 
44     decBuffer = new double[ m_inputLength ];
45 
46     // If adding new factors here, add them to
47     // getHighestSupportedFactor in the header as well
48 
49     if(m_decFactor == 8) {
50 
51         //////////////////////////////////////////////////
52         b[0] = 0.060111378492136;
53         b[1] = -0.257323420830598;
54         b[2] = 0.420583503165928;
55         b[3] = -0.222750785197418;
56         b[4] = -0.222750785197418;
57         b[5] = 0.420583503165928;
58         b[6] = -0.257323420830598;
59         b[7] = 0.060111378492136;
60 
61         a[0] = 1;
62         a[1] = -5.667654878577432;
63         a[2] = 14.062452278088417;
64         a[3] = -19.737303840697738;
65         a[4] = 16.889698874608641;
66         a[5] = -8.796600612325928;
67         a[6] = 2.577553446979888;
68         a[7] = -0.326903916815751;
69         //////////////////////////////////////////////////
70 
71     } else if( m_decFactor == 4 ) {
72 
73         //////////////////////////////////////////////////
74         b[ 0 ] = 0.10133306904918619;
75         b[ 1 ] = -0.2447523353702363;
76         b[ 2 ] = 0.33622528590120965;
77         b[ 3 ] = -0.13936581560633518;
78         b[ 4 ] = -0.13936581560633382;
79         b[ 5 ] = 0.3362252859012087;
80         b[ 6 ] = -0.2447523353702358;
81         b[ 7 ] = 0.10133306904918594;
82 
83         a[ 0 ] = 1;
84         a[ 1 ] = -3.9035590278139427;
85         a[ 2 ] = 7.5299379980621133;
86         a[ 3 ] = -8.6890803793177511;
87         a[ 4 ] = 6.4578667096099176;
88         a[ 5 ] = -3.0242979431223631;
89         a[ 6 ] = 0.83043385136748382;
90         a[ 7 ] = -0.094420800837809335;
91         //////////////////////////////////////////////////
92 
93     } else if( m_decFactor == 2 ) {
94 
95         //////////////////////////////////////////////////
96         b[ 0 ] = 0.20898944260075727;
97         b[ 1 ] = 0.40011234879814367;
98         b[ 2 ] = 0.819741973072733;
99         b[ 3 ] = 1.0087419911682323;
100         b[ 4 ] = 1.0087419911682325;
101         b[ 5 ] = 0.81974197307273156;
102         b[ 6 ] = 0.40011234879814295;
103         b[ 7 ] = 0.20898944260075661;
104 
105         a[ 0 ] = 1;
106         a[ 1 ] = 0.0077331184208358217;
107         a[ 2 ] = 1.9853971155964376;
108         a[ 3 ] = 0.19296739275341004;
109         a[ 4 ] = 1.2330748872852182;
110         a[ 5 ] = 0.18705341389316466;
111         a[ 6 ] = 0.23659265908013868;
112         a[ 7 ] = 0.032352924250533946;
113 
114     } else {
115 
116         if ( m_decFactor != 1 ) {
117             std::cerr << "WARNING: Decimator::initialise: unsupported decimation factor " << m_decFactor << ", no antialiasing filter will be used" << std::endl;
118         }
119 
120         //////////////////////////////////////////////////
121         b[ 0 ] = 1;
122         b[ 1 ] = 0;
123         b[ 2 ] = 0;
124         b[ 3 ] = 0;
125         b[ 4 ] = 0;
126         b[ 5 ] = 0;
127         b[ 6 ] = 0;
128         b[ 7 ] = 0;
129 
130         a[ 0 ] = 1;
131         a[ 1 ] = 0;
132         a[ 2 ] = 0;
133         a[ 3 ] = 0;
134         a[ 4 ] = 0;
135         a[ 5 ] = 0;
136         a[ 6 ] = 0;
137         a[ 7 ] = 0;
138     }
139 
140     resetFilter();
141 }
142 
deInitialise()143 void Decimator::deInitialise()
144 {
145     delete [] decBuffer;
146 }
147 
resetFilter()148 void Decimator::resetFilter()
149 {
150     Input = Output = 0;
151 
152     o1=o2=o3=o4=o5=o6=o7=0;
153 }
154 
doAntiAlias(const double * src,double * dst,int length)155 void Decimator::doAntiAlias(const double *src, double *dst, int length)
156 {
157     for (int i = 0; i < length; i++ ) {
158 
159         Input = (double)src[ i ];
160 
161         Output = Input * b[ 0 ] + o1;
162 
163         o1 = Input * b[ 1 ] - Output * a[ 1 ] + o2;
164         o2 = Input * b[ 2 ] - Output * a[ 2 ] + o3;
165         o3 = Input * b[ 3 ] - Output * a[ 3 ] + o4;
166         o4 = Input * b[ 4 ] - Output * a[ 4 ] + o5;
167         o5 = Input * b[ 5 ] - Output * a[ 5 ] + o6;
168         o6 = Input * b[ 6 ] - Output * a[ 6 ] + o7;
169         o7 = Input * b[ 7 ] - Output * a[ 7 ] ;
170 
171         dst[ i ] = Output;
172     }
173 }
174 
doAntiAlias(const float * src,double * dst,int length)175 void Decimator::doAntiAlias(const float *src, double *dst, int length)
176 {
177     for (int i = 0; i < length; i++ ) {
178 
179         Input = (double)src[ i ];
180 
181         Output = Input * b[ 0 ] + o1;
182 
183         o1 = Input * b[ 1 ] - Output * a[ 1 ] + o2;
184         o2 = Input * b[ 2 ] - Output * a[ 2 ] + o3;
185         o3 = Input * b[ 3 ] - Output * a[ 3 ] + o4;
186         o4 = Input * b[ 4 ] - Output * a[ 4 ] + o5;
187         o5 = Input * b[ 5 ] - Output * a[ 5 ] + o6;
188         o6 = Input * b[ 6 ] - Output * a[ 6 ] + o7;
189         o7 = Input * b[ 7 ] - Output * a[ 7 ] ;
190 
191         dst[ i ] = Output;
192     }
193 }
194 
process(const double * src,double * dst)195 void Decimator::process(const double *src, double *dst)
196 {
197     if (m_decFactor == 1) {
198         for (int i = 0; i < m_outputLength; i++ ) {
199             dst[i] = src[i];
200         }
201         return;
202     }
203 
204     doAntiAlias( src, decBuffer, m_inputLength );
205 
206     int idx = 0;
207 
208     for (int i = 0; i < m_outputLength; i++ ) {
209         dst[ idx++ ] = decBuffer[ m_decFactor * i ];
210     }
211 }
212 
process(const float * src,float * dst)213 void Decimator::process(const float *src, float *dst)
214 {
215     if (m_decFactor == 1) {
216         for (int i = 0; i < m_outputLength; i++ ) {
217             dst[i] = src[i];
218         }
219         return;
220     }
221 
222     doAntiAlias( src, decBuffer, m_inputLength );
223 
224     int idx = 0;
225 
226     for (int i = 0; i < m_outputLength; i++ ) {
227         dst[ idx++ ] = decBuffer[ m_decFactor * i ];
228     }
229 }
230