• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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