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
47const 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};
63unsigned 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 ****************************************************************************/
79FIR_struct* FIR_construct(unsigned int nTaps, const typeCoeff *pCoeffs, int scale, int factor_up, int factor_down)
80{
81FIR_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
124int 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
152void 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
228void FIR_downsample(unsigned int nInput, typeSample *pInput,
229                    typeSample *pOutput, FIR_struct *pFIR)
230{
231typeAccum        accum;
232typeCoeff const *ph;      // pointer to coefficients
233typeSample      *pz;      // pointer to delay line
234typeSample      *pz_beg;  // pointer to beginning of delay line
235typeSample      *pz_end;  // pointer to last slot in delay line
236unsigned int     nTaps_end;
237unsigned int     nTaps_beg;
238unsigned 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
296int main(int argc, char* argv[])
297{
298FILE         *fpInputSamples;    // input raw PCM file
299FILE         *fpOutputSamples;   // output raw PCM file
300
301typeSample    s_in[FACTOR_DOWN]; // input samples
302typeSample    s_out;             // filtered sample
303FIR_struct   *pFIR;              // pointer to FIR structure
304int           nSampleGet;        // number of samples to read from input speech file
305int           nSampleRead;       // number of samples read from input speech file
306unsigned long nSampleTot;        // total number of samples read so far
307time_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