/* <![CDATA[ */
2 * Copyright (C) 2013 The Android Open Source Project
3 *
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
9 *
10 * Unless required by applicable law or agreed to in writing, software
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#ifndef ANDROID_AUDIO_RESAMPLER_FIR_GEN_H
18#define ANDROID_AUDIO_RESAMPLER_FIR_GEN_H
19
20namespace android {
21
22/*
23 * generates a sine wave at equal steps.
24 *
25 * As most of our functions use sine or cosine at equal steps,
26 * it is very efficient to compute them that way (single multiply and subtract),
27 * rather than invoking the math library sin() or cos() each time.
28 *
29 * SineGen uses Goertzel's Algorithm (as a generator not a filter)
30 * to calculate sine(wstart + n * wstep) or cosine(wstart + n * wstep)
31 * by stepping through 0, 1, ... n.
32 *
33 * e^i(wstart+wstep) = 2cos(wstep) * e^i(wstart) - e^i(wstart-wstep)
34 *
35 * or looking at just the imaginary sine term, as the cosine follows identically:
36 *
37 * sin(wstart+wstep) = 2cos(wstep) * sin(wstart) - sin(wstart-wstep)
38 *
39 * Goertzel's algorithm is more efficient than the angle addition formula,
40 * e^i(wstart+wstep) = e^i(wstart) * e^i(wstep), which takes up to
41 * 4 multiplies and 2 adds (or 3* and 3+) and requires both sine and
42 * cosine generation due to the complex * complex multiply (full rotation).
43 *
44 * See: http://en.wikipedia.org/wiki/Goertzel_algorithm
45 *
46 */
47
48class SineGen {
49public:
50    SineGen(double wstart, double wstep, bool cosine = false) {
51        if (cosine) {
52            mCurrent = cos(wstart);
53            mPrevious = cos(wstart - wstep);
54        } else {
55            mCurrent = sin(wstart);
56            mPrevious = sin(wstart - wstep);
57        }
58        mTwoCos = 2.*cos(wstep);
59    }
60    SineGen(double expNow, double expPrev, double twoCosStep) {
61        mCurrent = expNow;
62        mPrevious = expPrev;
63        mTwoCos = twoCosStep;
64    }
65    inline double value() const {
66        return mCurrent;
67    }
69        double tmp = mCurrent;
70        mCurrent = mCurrent*mTwoCos - mPrevious;
71        mPrevious = tmp;
72    }
74        double tmp = mCurrent;
75        mCurrent = mCurrent*mTwoCos - mPrevious;
76        mPrevious = tmp;
77        return tmp;
78    }
79
80private:
81    double mCurrent; // current value of sine/cosine
82    double mPrevious; // previous value of sine/cosine
83    double mTwoCos; // stepping factor
84};
85
86/*
87 * generates a series of sine generators, phase offset by fixed steps.
88 *
89 * This is used to generate polyphase sine generators, one per polyphase
90 * in the filter code below.
91 *
92 * The SineGen returned by value() starts at innerStart = outerStart + n*outerStep;
93 * increments by innerStep.
94 *
95 */
96
97class SineGenGen {
98public:
99    SineGenGen(double outerStart, double outerStep, double innerStep, bool cosine = false)
100            : mSineInnerCur(outerStart, outerStep, cosine),
101              mSineInnerPrev(outerStart-innerStep, outerStep, cosine)
102    {
103        mTwoCos = 2.*cos(innerStep);
104    }
105    inline SineGen value() {
106        return SineGen(mSineInnerCur.value(), mSineInnerPrev.value(), mTwoCos);
107    }
111    }
114    }
115
116private:
117    SineGen mSineInnerCur; // generate the inner sine values (stepped by outerStep).
118    SineGen mSineInnerPrev; // generate the inner sine previous values
119                            // (behind by innerStep, stepped by outerStep).
120    double mTwoCos; // the inner stepping factor for the returned SineGen.
121};
122
123static inline double sqr(double x) {
124    return x * x;
125}
126
127/*
128 * rounds a double to the nearest integer for FIR coefficients.
129 *
130 * One variant uses noise shaping, which must keep error history
131 * to work (the err parameter, initialized to 0).
132 * The other variant is a non-noise shaped version for
133 * S32 coefficients (noise shaping doesn't gain much).
134 *
135 * Caution: No bounds saturation is applied, but isn't needed in this case.
136 *
137 * @param x is the value to round.
138 *
139 * @param maxval is the maximum integer scale factor expressed as an int64 (for headroom).
140 * Typically this may be the maximum positive integer+1 (using the fact that double precision
141 * FIR coefficients generated here are never that close to 1.0 to pose an overflow condition).
142 *
143 * @param err is the previous error (actual - rounded) for the previous rounding op.
144 * For 16b coefficients this can improve stopband dB performance by up to 2dB.
145 *
146 * Many variants exist for the noise shaping: http://en.wikipedia.org/wiki/Noise_shaping
147 *
148 */
149
150static inline int64_t toint(double x, int64_t maxval, double& err) {
151    double val = x * maxval;
152    double ival = floor(val + 0.5 + err*0.2);
153    err = val - ival;
154    return static_cast<int64_t>(ival);
155}
156
157static inline int64_t toint(double x, int64_t maxval) {
158    return static_cast<int64_t>(floor(x * maxval + 0.5));
159}
160
161/*
162 * Modified Bessel function of the first kind
163 * http://en.wikipedia.org/wiki/Bessel_function
164 *
165 * The formulas are taken from Abramowitz and Stegun,
166 * _Handbook of Mathematical Functions_ (links below):
167 *
168 * http://people.math.sfu.ca/~cbm/aands/page_375.htm
169 * http://people.math.sfu.ca/~cbm/aands/page_378.htm
170 *
171 * http://dlmf.nist.gov/10.25
172 * http://dlmf.nist.gov/10.40
173 *
174 * Note we assume x is nonnegative (the function is symmetric,
175 * pass in the absolute value as needed).
176 *
177 * Constants are compile time derived with templates I0Term<> and
178 * I0ATerm<> to the precision of the compiler.  The series can be expanded
179 * to any precision needed, but currently set around 24b precision.
180 *
181 * We use a bit of template math here, constexpr would probably be
182 * more appropriate for a C++11 compiler.
183 *
184 * For the intermediate range 3.75 < x < 15, we use minimax polynomial fit.
185 *
186 */
187
188template <int N>
189struct I0Term {
190    static const double value = I0Term<N-1>::value / (4. * N * N);
191};
192
193template <>
194struct I0Term<0> {
195    static const double value = 1.;
196};
197
198template <int N>
199struct I0ATerm {
200    static const double value = I0ATerm<N-1>::value * (2.*N-1.) * (2.*N-1.) / (8. * N);
201};
202
203template <>
204struct I0ATerm<0> { // 1/sqrt(2*PI);
205    static const double value = 0.398942280401432677939946059934381868475858631164934657665925;
206};
207
208#if USE_HORNERS_METHOD
209/* Polynomial evaluation of A + Bx + Cx^2 + Dx^3 + ...
210 * using Horner's Method: http://en.wikipedia.org/wiki/Horner's_method
211 *
212 * This has fewer multiplications than Estrin's method below, but has back to back
213 * floating point dependencies.
214 *
215 * On ARM this appears to work slower, so USE_HORNERS_METHOD is not default enabled.
216 */
217
218inline double Poly2(double A, double B, double x) {
219    return A + x * B;
220}
221
222inline double Poly4(double A, double B, double C, double D, double x) {
223    return A + x * (B + x * (C + x * (D)));
224}
225
226inline double Poly7(double A, double B, double C, double D, double E, double F, double G,
227        double x) {
228    return A + x * (B + x * (C + x * (D + x * (E + x * (F + x * (G))))));
229}
230
231inline double Poly9(double A, double B, double C, double D, double E, double F, double G,
232        double H, double I, double x) {
233    return A + x * (B + x * (C + x * (D + x * (E + x * (F + x * (G + x * (H + x * (I))))))));
234}
235
236#else
237/* Polynomial evaluation of A + Bx + Cx^2 + Dx^3 + ...
238 * using Estrin's Method: http://en.wikipedia.org/wiki/Estrin's_scheme
239 *
240 * This is typically faster, perhaps gains about 5-10% overall on ARM processors
241 * over Horner's method above.
242 */
243
244inline double Poly2(double A, double B, double x) {
245    return A + B * x;
246}
247
248inline double Poly3(double A, double B, double C, double x, double x2) {
249    return Poly2(A, B, x) + C * x2;
250}
251
252inline double Poly3(double A, double B, double C, double x) {
253    return Poly2(A, B, x) + C * x * x;
254}
255
256inline double Poly4(double A, double B, double C, double D, double x, double x2) {
257    return Poly2(A, B, x) + Poly2(C, D, x) * x2; // same as poly2(poly2, poly2, x2);
258}
259
260inline double Poly4(double A, double B, double C, double D, double x) {
261    return Poly4(A, B, C, D, x, x * x);
262}
263
264inline double Poly7(double A, double B, double C, double D, double E, double F, double G,
265        double x) {
266    double x2 = x * x;
267    return Poly4(A, B, C, D, x, x2) + Poly3(E, F, G, x, x2) * (x2 * x2);
268}
269
270inline double Poly8(double A, double B, double C, double D, double E, double F, double G,
271        double H, double x, double x2, double x4) {
272    return Poly4(A, B, C, D, x, x2) + Poly4(E, F, G, H, x, x2) * x4;
273}
274
275inline double Poly9(double A, double B, double C, double D, double E, double F, double G,
276        double H, double I, double x) {
277    double x2 = x * x;
278#if 1
279    // It does not seem faster to explicitly decompose Poly8 into Poly4, but
280    // could depend on compiler floating point scheduling.
281    double x4 = x2 * x2;
282    return Poly8(A, B, C, D, E, F, G, H, x, x2, x4) + I * (x4 * x4);
283#else
284    double val = Poly4(A, B, C, D, x, x2);
285    double x4 = x2 * x2;
286    return val + Poly4(E, F, G, H, x, x2) * x4 + I * (x4 * x4);
287#endif
288}
289#endif
290
291static inline double I0(double x) {
292    if (x < 3.75) {
293        x *= x;
294        return Poly7(I0Term<0>::value, I0Term<1>::value,
295                I0Term<2>::value, I0Term<3>::value,
296                I0Term<4>::value, I0Term<5>::value,
297                I0Term<6>::value, x); // e < 1.6e-7
298    }
299    if (1) {
300        /*
301         * Series expansion coefs are easy to calculate, but are expanded around 0,
302         * so error is unequal over the interval 0 < x < 3.75, the error being
303         * significantly better near 0.
304         *
305         * A better solution is to use precise minimax polynomial fits.
306         *
307         * We use a slightly more complicated solution for 3.75 < x < 15, based on
308         * the tables in Blair and Edwards, "Stable Rational Minimax Approximations
309         * to the Modified Bessel Functions I0(x) and I1(x)", Chalk Hill Nuclear Laboratory,
310         * AECL-4928.
311         *
312         * http://www.iaea.org/inis/collection/NCLCollectionStore/_Public/06/178/6178667.pdf
313         *
314         * See Table 11 for 0 < x < 15; e < 10^(-7.13).
315         *
316         * Note: Beta cannot exceed 15 (hence Stopband cannot exceed 144dB = 24b).
317         *
318         * This speeds up overall computation by about 40% over using the else clause below,
319         * which requires sqrt and exp.
320         *
321         */
322
323        x *= x;
324        double num = Poly9(-0.13544938430e9, -0.33153754512e8,
325                -0.19406631946e7, -0.48058318783e5,
326                -0.63269783360e3, -0.49520779070e1,
327                -0.24970910370e-1, -0.74741159550e-4,
328                -0.18257612460e-6, x);
329        double y = x - 225.; // reflection around 15 (squared)
330        double den = Poly4(-0.34598737196e8, 0.23852643181e6,
331                -0.70699387620e3, 0.10000000000e1, y);
332        return num / den;
333
334#if IO_EXTENDED_BETA
335        /* Table 42 for x > 15; e < 10^(-8.11).
336         * This is used for Beta>15, but is disabled here as
337         * we never use Beta that high.
338         *
339         * NOTE: This should be enabled only for x > 15.
340         */
341
342        double y = 1./x;
343        double z = y - (1./15);
344        double num = Poly2(0.415079861746e1, -0.5149092496e1, z);
345        double den = Poly3(0.103150763823e2, -0.14181687413e2,
346                0.1000000000e1, z);
347        return exp(x) * sqrt(y) * num / den;
348#endif
349    } else {
350        /*
351         * NOT USED, but reference for large Beta.
352         *
353         * Abramowitz and Stegun asymptotic formula.
354         * works for x > 3.75.
355         */
356        double y = 1./x;
357        return exp(x) * sqrt(y) *
358                // note: reciprocal squareroot may be easier!
359                // http://en.wikipedia.org/wiki/Fast_inverse_square_root
360                Poly9(I0ATerm<0>::value, I0ATerm<1>::value,
361                        I0ATerm<2>::value, I0ATerm<3>::value,
362                        I0ATerm<4>::value, I0ATerm<5>::value,
363                        I0ATerm<6>::value, I0ATerm<7>::value,
364                        I0ATerm<8>::value, y); // (... e) < 1.9e-7
365    }
366}
367
368/* A speed optimized version of the Modified Bessel I0() which incorporates
369 * the sqrt and numerator multiply and denominator divide into the computation.
370 * This speeds up filter computation by about 10-15%.
371 */
372static inline double I0SqrRat(double x2, double num, double den) {
373    if (x2 < (3.75 * 3.75)) {
374        return Poly7(I0Term<0>::value, I0Term<1>::value,
375                I0Term<2>::value, I0Term<3>::value,
376                I0Term<4>::value, I0Term<5>::value,
377                I0Term<6>::value, x2) * num / den; // e < 1.6e-7
378    }
379    num *= Poly9(-0.13544938430e9, -0.33153754512e8,
380            -0.19406631946e7, -0.48058318783e5,
381            -0.63269783360e3, -0.49520779070e1,
382            -0.24970910370e-1, -0.74741159550e-4,
383            -0.18257612460e-6, x2); // e < 10^(-7.13).
384    double y = x2 - 225.; // reflection around 15 (squared)
385    den *= Poly4(-0.34598737196e8, 0.23852643181e6,
386            -0.70699387620e3, 0.10000000000e1, y);
387    return num / den;
388}
389
390/*
391 * calculates the transition bandwidth for a Kaiser filter
392 *
393 * Formula 3.2.8, Vaidyanathan, _Multirate Systems and Filter Banks_, p. 48
394 * Formula 7.76, Oppenheim and Schafer, _Discrete-time Signal Processing, 3e_, p. 542
395 *
396 * @param halfNumCoef is half the number of coefficients per filter phase.
397 *
398 * @param stopBandAtten is the stop band attenuation desired.
399 *
400 * @return the transition bandwidth in normalized frequency (0 <= f <= 0.5)
401 */
402static inline double firKaiserTbw(int halfNumCoef, double stopBandAtten) {
403    return (stopBandAtten - 7.95)/((2.*14.36)*halfNumCoef);
404}
405
406/*
407 * calculates the fir transfer response of the overall polyphase filter at w.
408 *
409 * Calculates the DTFT transfer coefficient H(w) for 0 <= w <= PI, utilizing the
410 * fact that h[n] is symmetric (cosines only, no complex arithmetic).
411 *
412 * We use Goertzel's algorithm to accelerate the computation to essentially
413 * a single multiply and 2 adds per filter coefficient h[].
414 *
415 * Be careful be careful to consider that h[n] is the overall polyphase filter,
416 * with L phases, so rescaling H(w)/L is probably what you expect for "unity gain",
417 * as you only use one of the polyphases at a time.
418 */
419template <typename T>
420static inline double firTransfer(const T* coef, int L, int halfNumCoef, double w) {
421    double accum = static_cast<double>(coef[0])*0.5;  // "center coefficient" from first bank
422    coef += halfNumCoef;    // skip first filterbank (picked up by the last filterbank).
423#if SLOW_FIRTRANSFER
424    /* Original code for reference.  This is equivalent to the code below, but slower. */
425    for (int i=1 ; i<=L ; ++i) {
426        for (int j=0, ix=i ; j<halfNumCoef ; ++j, ix+=L) {
427            accum += cos(ix*w)*static_cast<double>(*coef++);
428        }
429    }
430#else
431    /*
432     * Our overall filter is stored striped by polyphases, not a contiguous h[n].
433     * We could fetch coefficients in a non-contiguous fashion
434     * but that will not scale to vector processing.
435     *
436     * We apply Goertzel's algorithm directly to each polyphase filter bank instead of
437     * using cosine generation/multiplication, thereby saving one multiply per inner loop.
438     *
439     * See: http://en.wikipedia.org/wiki/Goertzel_algorithm
440     * Also: Oppenheim and Schafer, _Discrete Time Signal Processing, 3e_, p. 720.
441     *
442     * We use the basic recursion to incorporate the cosine steps into real sequence x[n]:
443     * s[n] = x[n] + (2cosw)*s[n-1] + s[n-2]
444     *
445     * y[n] = s[n] - e^(iw)s[n-1]
446     *      = sum_{k=-\infty}^{n} x[k]e^(-iw(n-k))
447     *      = e^(-iwn) sum_{k=0}^{n} x[k]e^(iwk)
448     *
449     * The summation contains the frequency steps we want multiplied by the source
450     * (similar to a DTFT).
451     *
452     * Using symmetry, and just the real part (be careful, this must happen
453     * after any internal complex multiplications), the polyphase filterbank
454     * transfer function is:
455     *
456     * Hpp[n, w, w_0] = sum_{k=0}^{n} x[k] * cos(wk + w_0)
457     *                = Re{ e^(iwn + iw_0) y[n]}
458     *                = cos(wn+w_0) * s[n] - cos(w(n+1)+w_0) * s[n-1]
459     *
460     * using the fact that s[n] of real x[n] is real.
461     *
462     */
463    double dcos = 2. * cos(L*w);
464    int start = ((halfNumCoef)*L + 1);
465    SineGen cc((start - L) * w, w, true); // cosine
466    SineGen cp(start * w, w, true); // cosine
467    for (int i=1 ; i<=L ; ++i) {
468        double sc = 0;
469        double sp = 0;
470        for (int j=0 ; j<halfNumCoef ; ++j) {
471            double tmp = sc;
472            sc  = static_cast<double>(*coef++) + dcos*sc - sp;
473            sp = tmp;
474        }
475        // If we are awfully clever, we can apply Goertzel's algorithm
476        // again on the sc and sp sequences returned here.
478    }
479#endif
480    return accum*2.;
481}
482
483/*
484 * evaluates the minimum and maximum |H(f)| bound in a band region.
485 *
486 * This is usually done with equally spaced increments in the target band in question.
487 * The passband is often very small, and sampled that way. The stopband is often much
488 * larger.
489 *
490 * We use the fact that the overall polyphase filter has an additional bank at the end
491 * for interpolation; hence it is overspecified for the H(f) computation.  Thus the
492 * first polyphase is never actually checked, excepting its first term.
493 *
494 * In this code we use the firTransfer() evaluator above, which uses Goertzel's
495 * algorithm to calculate the transfer function at each point.
496 *
497 * TODO: An alternative with equal spacing is the FFT/DFT.  An alternative with unequal
498 * spacing is a chirp transform.
499 *
500 * @param coef is the designed polyphase filter banks
501 *
502 * @param L is the number of phases (for interpolation)
503 *
504 * @param halfNumCoef should be half the number of coefficients for a single
505 * polyphase.
506 *
507 * @param fstart is the normalized frequency start.
508 *
509 * @param fend is the normalized frequency end.
510 *
511 * @param steps is the number of steps to take (sampling) between frequency start and end
512 *
513 * @param firMin returns the minimum transfer |H(f)| found
514 *
515 * @param firMax returns the maximum transfer |H(f)| found
516 *
517 * 0 <= f <= 0.5.
518 * This is used to test passband and stopband performance.
519 */
520template <typename T>
521static void testFir(const T* coef, int L, int halfNumCoef,
522        double fstart, double fend, int steps, double &firMin, double &firMax) {
523    double wstart = fstart*(2.*M_PI);
524    double wend = fend*(2.*M_PI);
525    double wstep = (wend - wstart)/steps;
526    double fmax, fmin;
527    double trf = firTransfer(coef, L, halfNumCoef, wstart);
528    if (trf<0) {
529        trf = -trf;
530    }
531    fmin = fmax = trf;
532    wstart += wstep;
533    for (int i=1; i<steps; ++i) {
534        trf = firTransfer(coef, L, halfNumCoef, wstart);
535        if (trf<0) {
536            trf = -trf;
537        }
538        if (trf>fmax) {
539            fmax = trf;
540        }
541        else if (trf<fmin) {
542            fmin = trf;
543        }
544        wstart += wstep;
545    }
546    // renormalize - this is only needed for integer filter types
547    double norm = 1./((1ULL<<(sizeof(T)*8-1))*L);
548
549    firMin = fmin * norm;
550    firMax = fmax * norm;
551}
552
553/*
554 * evaluates the |H(f)| lowpass band characteristics.
555 *
556 * This function tests the lowpass characteristics for the overall polyphase filter,
557 * and is used to verify the design.  For this case, fp should be set to the
558 * passband normalized frequency from 0 to 0.5 for the overall filter (thus it
559 * is the designed polyphase bank value / L).  Likewise for fs.
560 *
561 * @param coef is the designed polyphase filter banks
562 *
563 * @param L is the number of phases (for interpolation)
564 *
565 * @param halfNumCoef should be half the number of coefficients for a single
566 * polyphase.
567 *
568 * @param fp is the passband normalized frequency, 0 < fp < fs < 0.5.
569 *
570 * @param fs is the stopband normalized frequency, 0 < fp < fs < 0.5.
571 *
572 * @param passSteps is the number of passband sampling steps.
573 *
574 * @param stopSteps is the number of stopband sampling steps.
575 *
576 * @param passMin is the minimum value in the passband
577 *
578 * @param passMax is the maximum value in the passband (useful for scaling).  This should
579 * be less than 1., to avoid sine wave test overflow.
580 *
581 * @param passRipple is the passband ripple.  Typically this should be less than 0.1 for
582 * an audio filter.  Generally speaker/headphone device characteristics will dominate
583 * the passband term.
584 *
585 * @param stopMax is the maximum value in the stopband.
586 *
587 * @param stopRipple is the stopband ripple, also known as stopband attenuation.
588 * Typically this should be greater than ~80dB for low quality, and greater than
589 * ~100dB for full 16b quality, otherwise aliasing may become noticeable.
590 *
591 */
592template <typename T>
593static void testFir(const T* coef, int L, int halfNumCoef,
594        double fp, double fs, int passSteps, int stopSteps,
595        double &passMin, double &passMax, double &passRipple,
596        double &stopMax, double &stopRipple) {
597    double fmin, fmax;
598    testFir(coef, L, halfNumCoef, 0., fp, passSteps, fmin, fmax);
599    double d1 = (fmax - fmin)/2.;
600    passMin = fmin;
601    passMax = fmax;
602    passRipple = -20.*log10(1. - d1); // passband ripple
603    testFir(coef, L, halfNumCoef, fs, 0.5, stopSteps, fmin, fmax);
604    // fmin is really not important for the stopband.
605    stopMax = fmax;
606    stopRipple = -20.*log10(fmax); // stopband ripple/attenuation
607}
608
609/*
610 * Calculates the overall polyphase filter based on a windowed sinc function.
611 *
612 * The windowed sinc is an odd length symmetric filter of exactly L*halfNumCoef*2+1
613 * taps for the entire kernel.  This is then decomposed into L+1 polyphase filterbanks.
614 * The last filterbank is used for interpolation purposes (and is mostly composed
615 * of the first bank shifted by one sample), and is unnecessary if one does
616 * not do interpolation.
617 *
618 * We use the last filterbank for some transfer function calculation purposes,
619 * so it needs to be generated anyways.
620 *
621 * @param coef is the caller allocated space for coefficients.  This should be
622 * exactly (L+1)*halfNumCoef in size.
623 *
624 * @param L is the number of phases (for interpolation)
625 *
626 * @param halfNumCoef should be half the number of coefficients for a single
627 * polyphase.
628 *
629 * @param stopBandAtten is the stopband value, should be >50dB.
630 *
631 * @param fcr is cutoff frequency/sampling rate (<0.5).  At this point, the energy
632 * should be 6dB less. (fcr is where the amplitude drops by half).  Use the
633 * firKaiserTbw() to calculate the transition bandwidth.  fcr is the midpoint
634 * between the stop band and the pass band (fstop+fpass)/2.
635 *
636 * @param atten is the attenuation (generally slightly less than 1).
637 */
638
639template <typename T>
640static inline void firKaiserGen(T* coef, int L, int halfNumCoef,
641        double stopBandAtten, double fcr, double atten) {
642    //
643    // Formula 3.2.5, 3.2.7, Vaidyanathan, _Multirate Systems and Filter Banks_, p. 48
644    // Formula 7.75, Oppenheim and Schafer, _Discrete-time Signal Processing, 3e_, p. 542
645    //
647    //
648    // Kaiser window and beta parameter
649    //
650    //         | 0.1102*(A - 8.7)                         A > 50
651    //  beta = | 0.5842*(A - 21)^0.4 + 0.07886*(A - 21)   21 <= A <= 50
652    //         | 0.                                       A < 21
653    //
654    // with A is the desired stop-band attenuation in dBFS
655    //
656    //    30 dB    2.210
657    //    40 dB    3.384
658    //    50 dB    4.538
659    //    60 dB    5.658
660    //    70 dB    6.764
661    //    80 dB    7.865
662    //    90 dB    8.960
663    //   100 dB   10.056
664
665    const int N = L * halfNumCoef; // non-negative half
666    const double beta = 0.1102 * (stopBandAtten - 8.7); // >= 50dB always
667    const double xstep = (2. * M_PI) * fcr / L;
668    const double xfrac = 1. / N;
669    const double yscale = atten * L / (I0(beta) * M_PI);
670    const double sqrbeta = sqr(beta);
671
672    // We use sine generators, which computes sines on regular step intervals.
673    // This speeds up overall computation about 40% from computing the sine directly.
674
675    SineGenGen sgg(0., xstep, L*xstep); // generates sine generators (one per polyphase)
676
677    for (int i=0 ; i<=L ; ++i) { // generate an extra set of coefs for interpolation
678
679        // computation for a single polyphase of the overall filter.
680        SineGen sg = sgg.valueAdvance(); // current sine generator for "j" inner loop.
681        double err = 0; // for noise shaping on int16_t coefficients (over each polyphase)
682
683        for (int j=0, ix=i ; j<halfNumCoef ; ++j, ix+=L) {
684            double y;
685            if (CC_LIKELY(ix)) {
686                double x = static_cast<double>(ix);
687
688                // sine generator: sg.valueAdvance() returns sin(ix*xstep);
689                // y = I0(beta * sqrt(1.0 - sqr(x * xfrac))) * yscale * sg.valueAdvance() / x;
690                y = I0SqrRat(sqrbeta * (1.0 - sqr(x * xfrac)), yscale * sg.valueAdvance(), x);
691            } else {
692                y = 2. * atten * fcr; // center of filter, sinc(0) = 1.
694            }
695
696            if (is_same<T, int16_t>::value) { // int16_t needs noise shaping
697                *coef++ = static_cast<T>(toint(y, 1ULL<<(sizeof(T)*8-1), err));
698            } else if (is_same<T, int32_t>::value) {
699                *coef++ = static_cast<T>(toint(y, 1ULL<<(sizeof(T)*8-1)));
700            } else { // assumed float or double
701                *coef++ = static_cast<T>(y);
702            }
703        }
704    }
705}
706
707}; // namespace android
708
709#endif /*ANDROID_AUDIO_RESAMPLER_FIR_GEN_H*/
710