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 ****************************************************************************/ 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 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 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 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
