fir.cpp revision 2a967b3fff07b8711aef41f839ad7521323bb64d 
1/*
2 * Copyright (C) 2007 Google Inc.
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 *      http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16
17#include <math.h>
18#include <stdio.h>
19
20static double sinc(double x) {
21    if (fabs(x) == 0.0f) return 1.0f;
22    return sin(x) / x;
23}
24
25static double sqr(double x) {
26    return x*x;
27}
28
29static double I0(double x) {
30    // from the Numerical Recipes in C p. 237
31    double ax,ans,y;
32    ax=fabs(x);
33    if (ax < 3.75) {
34        y=x/3.75;
35        y*=y;
36        ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
37            +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
38    } else {
39        y=3.75/ax;
40        ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1
41            +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
42            +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
43            +y*0.392377e-2))))))));
44    }
45    return ans;
46}
47
48static double kaiser(int k, int N, double alpha) {
49    if (k < 0 || k > N)
50        return 0;
51    return I0(M_PI*alpha * sqrt(1.0 - sqr((2.0*k)/N - 1.0))) / I0(M_PI*alpha);
52}
53
54int main(int argc, char** argv)
55{
56    // nc is the number of bits to store the coefficients
57    int nc = 32;
58
59    // ni is the minimum number of bits needed for interpolation
60    // (not used for generating the coefficients)
61    const int ni = nc / 2;
62
63    // cut off frequency ratio Fc/Fs
64    // The bigger the stop-band, the less coefficients we'll need.
65    double Fcr = 22050.0 / 48000.0;
66
67    // nzc is the number of zero-crossing on one half of the filter
68    int nzc = 16;
69
70    // alpha parameter of the kaiser window
71    // Larger numbers reduce ripples in the rejection band but increase
72    // the width of the transition band.
73    // the table below gives some value of alpha for a given
74    // stop-band attenuation.
75    //    30 dB    2.210
76    //    40 dB    3.384
77    //    50 dB    4.538
78    //    60 dB    5.658
79    //    70 dB    6.764
80    //    80 dB    7.865
81    //    90 dB    8.960
82    //   100 dB   10.056
83    double alpha = 10.056;	// -100dB stop-band attenuation
84
85    // 2^nz is the number coefficients per zero-crossing
86    // (int theory this should be 1<<(nc/2))
87    const int nz = 4;
88
89    // total number of coefficients
90    const int N = (1 << nz) * nzc;
91
92    // generate the right half of the filter
93    printf("const int32_t RESAMPLE_FIR_SIZE           = %d;\n", N);
94    printf("const int32_t RESAMPLE_FIR_NUM_COEF       = %d;\n", nzc);
95    printf("const int32_t RESAMPLE_FIR_COEF_BITS      = %d;\n", nc);
96    printf("const int32_t RESAMPLE_FIR_LERP_FRAC_BITS = %d;\n", ni);
97    printf("const int32_t RESAMPLE_FIR_LERP_INT_BITS  = %d;\n", nz);
98    printf("\n");
99    printf("static int16_t resampleFIR[%d] = {", N);
100    for (int i=0 ; i<N ; i++)
101    {
102        double x = (2.0 * M_PI * i * Fcr) / (1 << nz);
103        double y = kaiser(i+N, 2*N, alpha) * sinc(x);
104
105        long yi = floor(y * ((1ULL<<(nc-1))) + 0.5);
106        if (yi >= (1LL<<(nc-1))) yi = (1LL<<(nc-1))-1;
107
108        if ((i % (1 << 4)) == 0) printf("\n    ");
109        if (nc > 16)
110        	printf("0x%08x, ", int(yi));
111        else
112        	printf("0x%04x, ", int(yi)&0xFFFF);
113    }
114    printf("\n};\n");
115    return 0;
116 }
117
118// http://www.dsptutor.freeuk.com/KaiserFilterDesign/KaiserFilterDesign.html
119// http://www.csee.umbc.edu/help/sound/AFsp-V2R1/html/audio/ResampAudio.html
120
121