1 /*---------------------------------------------------------------------------*
2 * filter.c *
3 * *
4 * Copyright 2007, 2008 Nuance Communciations, Inc. *
5 * *
6 * Licensed under the Apache License, Version 2.0 (the 'License'); *
7 * you may not use this file except in compliance with the License. *
8 * *
9 * You may obtain a copy of the License at *
10 * http://www.apache.org/licenses/LICENSE-2.0 *
11 * *
12 * Unless required by applicable law or agreed to in writing, software *
13 * distributed under the License is distributed on an 'AS IS' BASIS, *
14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
15 * See the License for the specific language governing permissions and *
16 * limitations under the License. *
17 * *
18 *---------------------------------------------------------------------------*/
19
20 #define _ROUNDOFF
21
22 /****************************************************************************
23 * FILENAME
24 * pcm44pcm11.c
25 *
26 * DESCRIPTION
27 *
28 * Apply FIR filter to 44 kHz raw 16-bit PCM linear audio then
29 * downsample to 11 kHz.
30 ****************************************************************************/
31
32 #include <stdio.h>
33 #include <stdlib.h>
34 #include <string.h>
35 #include <ctype.h>
36 #include <time.h>
37
38 #include "filter.h"
39
40 /****************************************************************************
41 * FIR FILTER *
42 ****************************************************************************/
43
44 /* b = firls(120, [0 5000 6000 44000/2]/44000*2, [1 1 0 0]); */
45 /* bRounded = round(b*2^15); */
46
47 const typeCoeff ps16FilterCoeff_up1_down4[] = {
48
49 1, 9, 13, 10, -1, -16, -24, -18, 2, 25,
50 38, 28, -2, -38, -57, -42, 3, 55, 81, 60,
51 -3, -76, -113, -84, 4, 104, 153, 114, -5, -138,
52 -205, -152, 5, 183, 271, 202, -6, -241, -360, -269,
53 6, 321, 482, 364, -6, -438, -667, -511, 7, 632,
54 986, 778, -7, -1030, -1704, -1450, 7, 2451, 5204, 7366,
55 8185, 7366, 5204, 2451, 7, -1450, -1704, -1030, -7, 778,
56 986, 632, 7, -511, -667, -438, -6, 364, 482, 321,
57 6, -269, -360, -241, -6, 202, 271, 183, 5, -152,
58 -205, -138, -5, 114, 153, 104, 4, -84, -113, -76,
59 -3, 60, 81, 55, 3, -42, -57, -38, -2, 28,
60 38, 25, 2, -18, -24, -16, -1, 10, 13, 9,
61 1,
62 };
63 unsigned int filter_length = sizeof(ps16FilterCoeff_up1_down4)/sizeof(typeCoeff);
64
65 /****************************************************************************
66 * FIR_struct *FIR_construct(unsigned int nTaps, FIR_type *pCoeffs)
67 *
68 * DESCRIPTION
69 * allocate and initialize FIR structure
70 *
71 * INPUT
72 * nTaps - length of FIR filter
73 * pCoeffs - pointer to FIR filter coefficients
74 * scale - fixed point scale factor
75 *
76 * OUTPUT
77 * returns pointer to FIR structure; NULL if error
78 ****************************************************************************/
FIR_construct(unsigned int nTaps,const typeCoeff * pCoeffs,int scale,int factor_up,int factor_down)79 FIR_struct* FIR_construct(unsigned int nTaps, const typeCoeff *pCoeffs, int scale, int factor_up, int factor_down)
80 {
81 FIR_struct *pFIR;
82
83 if (nTaps == 0)
84 return NULL;
85
86 pFIR = malloc(sizeof(FIR_struct));
87 if (pFIR == NULL)
88 return NULL;
89
90 // alloocate space for delay line use calloc to avoid uninitialized memory
91 // that causes an audible "pop" at the beginning of audio. SteveR
92 pFIR->z = calloc(nTaps * sizeof(typeSample), 1);
93 if (pFIR->z == NULL)
94 {
95 free(pFIR);
96 return NULL;
97 }
98
99 pFIR->factor_up = factor_up;
100 pFIR->factor_down = factor_down;
101
102 pFIR->state = 0;
103 pFIR->h = pCoeffs;
104 pFIR->nTaps = nTaps;
105 pFIR->scale = scale;
106 pFIR->round = (1 << (scale-1));
107
108 return pFIR;
109 }
110
111 /****************************************************************************
112 * int FIR_deconstruct(FIR_struct *pFIR)
113 *
114 * DESCRIPTION
115 * deallocate FIR structure
116 *
117 * INPUT
118 * pFIR - pointer to FIR structure
119 *
120 * OUTPUT
121 * returns 0 for success; 1 for failure
122 ****************************************************************************/
123
FIR_deconstruct(FIR_struct * pFIR)124 int FIR_deconstruct(FIR_struct *pFIR)
125 {
126 if (pFIR == NULL)
127 return 1;
128
129 if (pFIR->z == NULL)
130 return 1;
131
132 free(pFIR->z);
133 free(pFIR);
134
135 return 0;
136 }
137
138 /****************************************************************************
139 * void FIR_reset(FIR_struct *pFIR)
140 *
141 * DESCRIPTION
142 *
143 * reset FIR state
144 *
145 * INPUT
146 * pFIR - pointer to FIR structure
147 *
148 * OUTPUT
149 * pFIR->state initialized
150 ****************************************************************************/
151
FIR_reset(FIR_struct * pFIR)152 void FIR_reset(FIR_struct *pFIR)
153 {
154 pFIR->state = 0;
155
156 memset(pFIR->z, pFIR->nTaps, sizeof(typeSample));
157 }
158
159 /*****************************************************************************
160 * FIR_type FIR_downsample(unsigned int nInput, typeSample *pInput,
161 * typeSample *pOutput, FIR_struct *pFIR)
162 *
163 * DESCRIPTION
164 *
165 * Apply FIR filter to input data. If nInput > 1, this will also
166 * decimate by a factor of nInput. That is, the filter will only be
167 * evaluated every nInput samples, not at each of the nInput samples.
168 *
169 * Breakup filter computation into 2 parts to avoid doing a wraparound
170 * check inside the loop.
171 *
172 * Example:
173 *
174 * pFIR->nTaps = 8
175 * pFIR->state = 2
176 * nInput = 1
177 * *pInput = s20
178 *
179 * (a) Store new sample(s) in delay buffer z[]
180 *
181 * Since pFIR->state == 2, store new sample s20 at location z[2]
182 *
183 * *** latest input stored at z[pFIR->state]
184 * *
185 * -------------------------------------------------
186 * z | s14 | s13 | s20 | s19 | s18 | s17 | s16 | s15 |
187 * -------------------------------------------------
188 * h | h0 | h1 | h2 | h3 | h4 | h5 | h6 | h7 |
189 * -------------------------------------------------
190 * 0 1 2 3 4 5 6 7
191 *
192 * (b) Update state to point to newest sample, wrap if < 0
193 *
194 * Since nInput == 1, state for newest sample is still 2
195 * (otherwise, update state -= nInput-1; wrap by adding nTaps if < 0)
196 *
197 * (c) Accumulate "end part" first
198 *
199 * z: start with latest sample at z[pFIR->state], then advance to right
200 * h: start with 1st filter coefficient, then advance to right
201 *
202 * acc = h0*s20 + h1*s19 + h2*s18 + h3*s17 + h4*s16 + h5*s15
203 *
204 * (d) Accumulate "beginning part"
205
206 * z: start with sample at beginning of delay buffer, then advance
207 * to sample before latest one at z[pFIR->state]
208 * h: continue with next filter coefficient from step (a)
209 *
210 * acc += (h6*s14 + h7*s13) FIR filter output
211 * *pOutput = acc
212 *
213 * (e) Update FIR state
214 *
215 * state--, wrapping if < 0 to simulate circular buffer
216 *
217 * INPUT
218 *
219 * nInput - number of new input samples; evaluate FIR at this point
220 * pInput - pointer to input sample buffer
221 * pOutput - pointer to output sample buffer
222 * pFIR - pointer to FIR structure
223 *
224 * OUTPUT
225 *
226 *****************************************************************************/
227
FIR_downsample(unsigned int nInput,typeSample * pInput,typeSample * pOutput,FIR_struct * pFIR)228 void FIR_downsample(unsigned int nInput, typeSample *pInput,
229 typeSample *pOutput, FIR_struct *pFIR)
230 {
231 typeAccum accum;
232 typeCoeff const *ph; // pointer to coefficients
233 typeSample *pz; // pointer to delay line
234 typeSample *pz_beg; // pointer to beginning of delay line
235 typeSample *pz_end; // pointer to last slot in delay line
236 unsigned int nTaps_end;
237 unsigned int nTaps_beg;
238 unsigned int i;
239
240 // initialize
241 accum = 0;
242 ph = pFIR->h; // point to coefficients
243 pz_beg = pFIR->z; // start of delay line
244 pz_end = pFIR->z + pFIR->nTaps - 1; // end of delay line
245
246 // (a) Store new input samples in delay line (circular addressing would help a lot)
247 pz = pFIR->z + pFIR->state; // point to next empty slot in delay line
248 for (i = 0; i < nInput; i++)
249 {
250 *pz-- = *pInput++;
251 if (pz < pz_beg)
252 pz = pz_end; // wrap around (circular buffer)
253 }
254
255 // (b) adjust state to reflect addition of samples
256 pFIR->state -= nInput-1;
257 if (pFIR->state < 0)
258 pFIR->state += pFIR->nTaps; // wrap
259
260 // (c) Accumulate "end part"
261 pz = pFIR->z + pFIR->state;
262 nTaps_end = pFIR->nTaps - pFIR->state;
263 for (i = 0; i < nTaps_end; i++)
264 {
265 accum += *ph++ * *pz++;
266 }
267
268 // (d) Accumulate "beginning part"
269 pz = pFIR->z;
270 nTaps_beg = pFIR->state;
271 for (i = 0; i < nTaps_beg; i++)
272 {
273 accum += *ph++ * *pz++;
274 }
275
276 // (e) Update FIR state for next batch of incoming samples
277 pFIR->state--;
278 if (pFIR->state < 0)
279 pFIR->state += pFIR->nTaps; // wrap
280
281 #ifdef _ROUNDOFF
282 if (accum >= 0)
283 accum += pFIR->round;
284 else
285 accum -= pFIR->round;
286 #endif
287
288 *pOutput = (typeSample) (accum >> pFIR->scale);
289 }
290
291 #if 0
292 /*****************************************************************************
293 * main
294 *****************************************************************************/
295
296 int main(int argc, char* argv[])
297 {
298 FILE *fpInputSamples; // input raw PCM file
299 FILE *fpOutputSamples; // output raw PCM file
300
301 typeSample s_in[FACTOR_DOWN]; // input samples
302 typeSample s_out; // filtered sample
303 FIR_struct *pFIR; // pointer to FIR structure
304 int nSampleGet; // number of samples to read from input speech file
305 int nSampleRead; // number of samples read from input speech file
306 unsigned long nSampleTot; // total number of samples read so far
307 time_t t0; // time upon entry
308
309 t0 = time(NULL); // get time upon entry
310
311 // Check Command-line Parameters
312 if (argc != 3)
313 {
314 fprintf(stderr, "pcm44pcm11 v1.0\n");
315 fprintf(stderr, " - downsamples 44 kHz to 11 kHz (16-bit PCM, Intel byte order)\n");
316 fprintf(stderr, "Usage: pcm44pcm11 <input PCM file> <output PCM file>\n\n");
317 return 0;
318 }
319
320 // Open input sample file
321 if ((fpInputSamples = fopen(argv[1], "rb")) == NULL)
322 {
323 fprintf(stderr, "Error reading input sample file: %s\n\n", argv[1]);
324 exit(1); // abnormal exit
325 }
326
327 // Create output sample file
328 if ((fpOutputSamples = fopen(argv[2], "wb")) == NULL)
329 {
330 fprintf(stderr, "Error creating output file: %s\n\n", argv[2]);
331 fclose(fpInputSamples);
332 exit(1); // abnormal exit
333 }
334
335 // **************************************************************************************
336 // Begin filtering...
337 // **************************************************************************************
338
339 pFIR = FIR_construct(filter_length,
340 ps16FilterCoeff_up1_down4,
341 u16ScaleFilterCoeff_up1_down4,
342 FACTOR_UP,
343 FACTOR_DOWN);
344
345 fprintf(stdout, "Filtering...\n");
346
347 FIR_reset(pFIR);
348
349 nSampleTot = 0;
350 while (!feof(fpInputSamples))
351 {
352 nSampleGet = pFIR->factor_down; // if downsampling, only filter every factor_down samples
353 nSampleRead = fread(s_in, sizeof(typeSample), nSampleGet, fpInputSamples);
354 if (feof(fpInputSamples) || (nSampleRead != nSampleGet))
355 break; // done with input file
356 nSampleTot += nSampleRead;
357
358 FIR_downsample(nSampleRead, s_in, &s_out, pFIR);
359
360 if (nSampleTot < pFIR->nTaps)
361 continue; // wait until delay buffer has been filled to skip transients
362
363 fwrite(&s_out, sizeof(typeSample), 1, fpOutputSamples);
364 }
365
366 fprintf(stdout, "\n\nTime elapsed: %d sec\n", time(NULL)-t0);
367
368 FIR_deconstruct(pFIR);
369
370 fclose(fpInputSamples);
371 fclose(fpOutputSamples);
372
373 return 0;
374 }
375
376 #endif
377