1 /*
2  * Copyright (C) 2013 The Android Open Source Project
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 #ifndef ANDROID_AUDIO_RESAMPLER_FIR_PROCESS_H
18 #define ANDROID_AUDIO_RESAMPLER_FIR_PROCESS_H
19 
20 namespace android {
21 
22 // depends on AudioResamplerFirOps.h
23 
24 /* variant for input type TI = int16_t input samples */
25 template<typename TC>
26 static inline
mac(int32_t & l,int32_t & r,TC coef,const int16_t * samples)27 void mac(int32_t& l, int32_t& r, TC coef, const int16_t* samples)
28 {
29     uint32_t rl = *reinterpret_cast<const uint32_t*>(samples);
30     l = mulAddRL(1, rl, coef, l);
31     r = mulAddRL(0, rl, coef, r);
32 }
33 
34 template<typename TC>
35 static inline
mac(int32_t & l,TC coef,const int16_t * samples)36 void mac(int32_t& l, TC coef, const int16_t* samples)
37 {
38     l = mulAdd(samples[0], coef, l);
39 }
40 
41 /* variant for input type TI = float input samples */
42 template<typename TC>
43 static inline
mac(float & l,float & r,TC coef,const float * samples)44 void mac(float& l, float& r, TC coef,  const float* samples)
45 {
46     l += *samples++ * coef;
47     r += *samples * coef;
48 }
49 
50 template<typename TC>
51 static inline
mac(float & l,TC coef,const float * samples)52 void mac(float& l, TC coef,  const float* samples)
53 {
54     l += *samples * coef;
55 }
56 
57 /* variant for output type TO = int32_t output samples */
58 static inline
volumeAdjust(int32_t value,int32_t volume)59 int32_t volumeAdjust(int32_t value, int32_t volume)
60 {
61     return 2 * mulRL(0, value, volume);  // Note: only use top 16b
62 }
63 
64 /* variant for output type TO = float output samples */
65 static inline
volumeAdjust(float value,float volume)66 float volumeAdjust(float value, float volume)
67 {
68     return value * volume;
69 }
70 
71 /*
72  * Helper template functions for loop unrolling accumulator operations.
73  *
74  * Unrolling the loops achieves about 2x gain.
75  * Using a recursive template rather than an array of TO[] for the accumulator
76  * values is an additional 10-20% gain.
77  */
78 
79 template<int CHANNELS, typename TO>
80 class Accumulator : public Accumulator<CHANNELS-1, TO> // recursive
81 {
82 public:
clear()83     inline void clear() {
84         value = 0;
85         Accumulator<CHANNELS-1, TO>::clear();
86     }
87     template<typename TC, typename TI>
acc(TC coef,const TI * & data)88     inline void acc(TC coef, const TI*& data) {
89         mac(value, coef, data++);
90         Accumulator<CHANNELS-1, TO>::acc(coef, data);
91     }
volume(TO * & out,TO gain)92     inline void volume(TO*& out, TO gain) {
93         *out++ = volumeAdjust(value, gain);
94         Accumulator<CHANNELS-1, TO>::volume(out, gain);
95     }
96 
97     TO value; // one per recursive inherited base class
98 };
99 
100 template<typename TO>
101 class Accumulator<0, TO> {
102 public:
clear()103     inline void clear() {
104     }
105     template<typename TC, typename TI>
acc(TC coef __unused,const TI * & data __unused)106     inline void acc(TC coef __unused, const TI*& data __unused) {
107     }
volume(TO * & out __unused,TO gain __unused)108     inline void volume(TO*& out __unused, TO gain __unused) {
109     }
110 };
111 
112 template<typename TC, typename TINTERP>
113 inline
interpolate(TC coef_0,TC coef_1,TINTERP lerp)114 TC interpolate(TC coef_0, TC coef_1, TINTERP lerp)
115 {
116     return lerp * (coef_1 - coef_0) + coef_0;
117 }
118 
119 template<>
120 inline
121 int16_t interpolate<int16_t, uint32_t>(int16_t coef_0, int16_t coef_1, uint32_t lerp)
122 {   // in some CPU architectures 16b x 16b multiplies are faster.
123     return (static_cast<int16_t>(lerp) * static_cast<int16_t>(coef_1 - coef_0) >> 15) + coef_0;
124 }
125 
126 template<>
127 inline
128 int32_t interpolate<int32_t, uint32_t>(int32_t coef_0, int32_t coef_1, uint32_t lerp)
129 {
130     return (lerp * static_cast<int64_t>(coef_1 - coef_0) >> 31) + coef_0;
131 }
132 
133 /* class scope for passing in functions into templates */
134 struct InterpCompute {
135     template<typename TC, typename TINTERP>
136     static inline
interpolatepInterpCompute137     TC interpolatep(TC coef_0, TC coef_1, TINTERP lerp) {
138         return interpolate(coef_0, coef_1, lerp);
139     }
140 
141     template<typename TC, typename TINTERP>
142     static inline
interpolatenInterpCompute143     TC interpolaten(TC coef_0, TC coef_1, TINTERP lerp) {
144         return interpolate(coef_0, coef_1, lerp);
145     }
146 };
147 
148 struct InterpNull {
149     template<typename TC, typename TINTERP>
150     static inline
interpolatepInterpNull151     TC interpolatep(TC coef_0, TC coef_1 __unused, TINTERP lerp __unused) {
152         return coef_0;
153     }
154 
155     template<typename TC, typename TINTERP>
156     static inline
interpolatenInterpNull157     TC interpolaten(TC coef_0 __unused, TC coef_1, TINTERP lerp __unused) {
158         return coef_1;
159     }
160 };
161 
162 /*
163  * Calculates a single output frame (two samples).
164  *
165  * The Process*() functions compute both the positive half FIR dot product and
166  * the negative half FIR dot product, accumulates, and then applies the volume.
167  *
168  * Use fir() to compute the proper coefficient pointers for a polyphase
169  * filter bank.
170  *
171  * ProcessBase() is the fundamental processing template function.
172  *
173  * ProcessL() calls ProcessBase() with TFUNC = InterpNull, for fixed/locked phase.
174  * Process() calls ProcessBase() with TFUNC = InterpCompute, for interpolated phase.
175  */
176 
177 template <int CHANNELS, int STRIDE, typename TFUNC, typename TC, typename TI, typename TO, typename TINTERP>
178 static inline
ProcessBase(TO * const out,size_t count,const TC * coefsP,const TC * coefsN,const TI * sP,const TI * sN,TINTERP lerpP,const TO * const volumeLR)179 void ProcessBase(TO* const out,
180         size_t count,
181         const TC* coefsP,
182         const TC* coefsN,
183         const TI* sP,
184         const TI* sN,
185         TINTERP lerpP,
186         const TO* const volumeLR)
187 {
188     COMPILE_TIME_ASSERT_FUNCTION_SCOPE(CHANNELS > 0)
189 
190     if (CHANNELS > 2) {
191         // TO accum[CHANNELS];
192         Accumulator<CHANNELS, TO> accum;
193 
194         // for (int j = 0; j < CHANNELS; ++j) accum[j] = 0;
195         accum.clear();
196         for (size_t i = 0; i < count; ++i) {
197             TC c = TFUNC::interpolatep(coefsP[0], coefsP[count], lerpP);
198 
199             // for (int j = 0; j < CHANNELS; ++j) mac(accum[j], c, sP + j);
200             const TI *tmp_data = sP; // tmp_ptr seems to work better
201             accum.acc(c, tmp_data);
202 
203             coefsP++;
204             sP -= CHANNELS;
205             c = TFUNC::interpolaten(coefsN[count], coefsN[0], lerpP);
206 
207             // for (int j = 0; j < CHANNELS; ++j) mac(accum[j], c, sN + j);
208             tmp_data = sN; // tmp_ptr seems faster than directly using sN
209             accum.acc(c, tmp_data);
210 
211             coefsN++;
212             sN += CHANNELS;
213         }
214         // for (int j = 0; j < CHANNELS; ++j) out[j] += volumeAdjust(accum[j], volumeLR[0]);
215         TO *tmp_out = out; // may remove if const out definition changes.
216         accum.volume(tmp_out, volumeLR[0]);
217     } else if (CHANNELS == 2) {
218         TO l = 0;
219         TO r = 0;
220         for (size_t i = 0; i < count; ++i) {
221             mac(l, r, TFUNC::interpolatep(coefsP[0], coefsP[count], lerpP), sP);
222             coefsP++;
223             sP -= CHANNELS;
224             mac(l, r, TFUNC::interpolaten(coefsN[count], coefsN[0], lerpP), sN);
225             coefsN++;
226             sN += CHANNELS;
227         }
228         out[0] += volumeAdjust(l, volumeLR[0]);
229         out[1] += volumeAdjust(r, volumeLR[1]);
230     } else { /* CHANNELS == 1 */
231         TO l = 0;
232         for (size_t i = 0; i < count; ++i) {
233             mac(l, TFUNC::interpolatep(coefsP[0], coefsP[count], lerpP), sP);
234             coefsP++;
235             sP -= CHANNELS;
236             mac(l, TFUNC::interpolaten(coefsN[count], coefsN[0], lerpP), sN);
237             coefsN++;
238             sN += CHANNELS;
239         }
240         out[0] += volumeAdjust(l, volumeLR[0]);
241         out[1] += volumeAdjust(l, volumeLR[1]);
242     }
243 }
244 
245 template <int CHANNELS, int STRIDE, typename TC, typename TI, typename TO>
246 static inline
ProcessL(TO * const out,int count,const TC * coefsP,const TC * coefsN,const TI * sP,const TI * sN,const TO * const volumeLR)247 void ProcessL(TO* const out,
248         int count,
249         const TC* coefsP,
250         const TC* coefsN,
251         const TI* sP,
252         const TI* sN,
253         const TO* const volumeLR)
254 {
255     ProcessBase<CHANNELS, STRIDE, InterpNull>(out, count, coefsP, coefsN, sP, sN, 0, volumeLR);
256 }
257 
258 template <int CHANNELS, int STRIDE, typename TC, typename TI, typename TO, typename TINTERP>
259 static inline
Process(TO * const out,int count,const TC * coefsP,const TC * coefsN,const TC * coefsP1 __unused,const TC * coefsN1 __unused,const TI * sP,const TI * sN,TINTERP lerpP,const TO * const volumeLR)260 void Process(TO* const out,
261         int count,
262         const TC* coefsP,
263         const TC* coefsN,
264         const TC* coefsP1 __unused,
265         const TC* coefsN1 __unused,
266         const TI* sP,
267         const TI* sN,
268         TINTERP lerpP,
269         const TO* const volumeLR)
270 {
271     ProcessBase<CHANNELS, STRIDE, InterpCompute>(out, count, coefsP, coefsN, sP, sN, lerpP, volumeLR);
272 }
273 
274 /*
275  * Calculates a single output frame (two samples) from input sample pointer.
276  *
277  * This sets up the params for the accelerated Process() and ProcessL()
278  * functions to do the appropriate dot products.
279  *
280  * @param out should point to the output buffer with space for at least one output frame.
281  *
282  * @param phase is the fractional distance between input frames for interpolation:
283  * phase >= 0  && phase < phaseWrapLimit.  It can be thought of as a rational fraction
284  * of phase/phaseWrapLimit.
285  *
286  * @param phaseWrapLimit is #polyphases<<coefShift, where #polyphases is the number of polyphases
287  * in the polyphase filter. Likewise, #polyphases can be obtained as (phaseWrapLimit>>coefShift).
288  *
289  * @param coefShift gives the bit alignment of the polyphase index in the phase parameter.
290  *
291  * @param halfNumCoefs is the half the number of coefficients per polyphase filter. Since the
292  * overall filterbank is odd-length symmetric, only halfNumCoefs need be stored.
293  *
294  * @param coefs is the polyphase filter bank, starting at from polyphase index 0, and ranging to
295  * and including the #polyphases.  Each polyphase of the filter has half-length halfNumCoefs
296  * (due to symmetry).  The total size of the filter bank in coefficients is
297  * (#polyphases+1)*halfNumCoefs.
298  *
299  * The filter bank coefs should be aligned to a minimum of 16 bytes (preferrably to cache line).
300  *
301  * The coefs should be attenuated (to compensate for passband ripple)
302  * if storing back into the native format.
303  *
304  * @param samples are unaligned input samples.  The position is in the "middle" of the
305  * sample array with respect to the FIR filter:
306  * the negative half of the filter is dot product from samples+1 to samples+halfNumCoefs;
307  * the positive half of the filter is dot product from samples to samples-halfNumCoefs+1.
308  *
309  * @param volumeLR is a pointer to an array of two 32 bit volume values, one per stereo channel,
310  * expressed as a S32 integer.  A negative value inverts the channel 180 degrees.
311  * The pointer volumeLR should be aligned to a minimum of 8 bytes.
312  * A typical value for volume is 0x1000 to align to a unity gain output of 20.12.
313  *
314  * In between calls to filterCoefficient, the phase is incremented by phaseIncrement, where
315  * phaseIncrement is calculated as inputSampling * phaseWrapLimit / outputSampling.
316  *
317  * The filter polyphase index is given by indexP = phase >> coefShift. Due to
318  * odd length symmetric filter, the polyphase index of the negative half depends on
319  * whether interpolation is used.
320  *
321  * The fractional siting between the polyphase indices is given by the bits below coefShift:
322  *
323  * lerpP = phase << 32 - coefShift >> 1;  // for 32 bit unsigned phase multiply
324  * lerpP = phase << 32 - coefShift >> 17; // for 16 bit unsigned phase multiply
325  *
326  * For integer types, this is expressed as:
327  *
328  * lerpP = phase << sizeof(phase)*8 - coefShift
329  *              >> (sizeof(phase)-sizeof(*coefs))*8 + 1;
330  *
331  * For floating point, lerpP is the fractional phase scaled to [0.0, 1.0):
332  *
333  * lerpP = (phase << 32 - coefShift) / (1 << 32); // floating point equivalent
334  */
335 
336 template<int CHANNELS, bool LOCKED, int STRIDE, typename TC, typename TI, typename TO>
337 static inline
fir(TO * const out,const uint32_t phase,const uint32_t phaseWrapLimit,const int coefShift,const int halfNumCoefs,const TC * const coefs,const TI * const samples,const TO * const volumeLR)338 void fir(TO* const out,
339         const uint32_t phase, const uint32_t phaseWrapLimit,
340         const int coefShift, const int halfNumCoefs, const TC* const coefs,
341         const TI* const samples, const TO* const volumeLR)
342 {
343     // NOTE: be very careful when modifying the code here. register
344     // pressure is very high and a small change might cause the compiler
345     // to generate far less efficient code.
346     // Always sanity check the result with objdump or test-resample.
347 
348     if (LOCKED) {
349         // locked polyphase (no interpolation)
350         // Compute the polyphase filter index on the positive and negative side.
351         uint32_t indexP = phase >> coefShift;
352         uint32_t indexN = (phaseWrapLimit - phase) >> coefShift;
353         const TC* coefsP = coefs + indexP*halfNumCoefs;
354         const TC* coefsN = coefs + indexN*halfNumCoefs;
355         const TI* sP = samples;
356         const TI* sN = samples + CHANNELS;
357 
358         // dot product filter.
359         ProcessL<CHANNELS, STRIDE>(out,
360                 halfNumCoefs, coefsP, coefsN, sP, sN, volumeLR);
361     } else {
362         // interpolated polyphase
363         // Compute the polyphase filter index on the positive and negative side.
364         uint32_t indexP = phase >> coefShift;
365         uint32_t indexN = (phaseWrapLimit - phase - 1) >> coefShift; // one's complement.
366         const TC* coefsP = coefs + indexP*halfNumCoefs;
367         const TC* coefsN = coefs + indexN*halfNumCoefs;
368         const TC* coefsP1 = coefsP + halfNumCoefs;
369         const TC* coefsN1 = coefsN + halfNumCoefs;
370         const TI* sP = samples;
371         const TI* sN = samples + CHANNELS;
372 
373         // Interpolation fraction lerpP derived by shifting all the way up and down
374         // to clear the appropriate bits and align to the appropriate level
375         // for the integer multiply.  The constants should resolve in compile time.
376         //
377         // The interpolated filter coefficient is derived as follows for the pos/neg half:
378         //
379         // interpolated[P] = index[P]*lerpP + index[P+1]*(1-lerpP)
380         // interpolated[N] = index[N+1]*lerpP + index[N]*(1-lerpP)
381 
382         // on-the-fly interpolated dot product filter
383         if (is_same<TC, float>::value || is_same<TC, double>::value) {
384             static const TC scale = 1. / (65536. * 65536.); // scale phase bits to [0.0, 1.0)
385             TC lerpP = TC(phase << (sizeof(phase)*8 - coefShift)) * scale;
386 
387             Process<CHANNELS, STRIDE>(out,
388                     halfNumCoefs, coefsP, coefsN, coefsP1, coefsN1, sP, sN, lerpP, volumeLR);
389         } else {
390             uint32_t lerpP = phase << (sizeof(phase)*8 - coefShift)
391                     >> ((sizeof(phase)-sizeof(*coefs))*8 + 1);
392 
393             Process<CHANNELS, STRIDE>(out,
394                     halfNumCoefs, coefsP, coefsN, coefsP1, coefsN1, sP, sN, lerpP, volumeLR);
395         }
396     }
397 }
398 
399 }; // namespace android
400 
401 #endif /*ANDROID_AUDIO_RESAMPLER_FIR_PROCESS_H*/
402