xref: /aosp_15_r20/external/webrtc/common_audio/signal_processing/splitting_filter.c (revision d9f758449e529ab9291ac668be2861e7a55c2422)
1 /*
2  *  Copyright (c) 2011 The WebRTC project authors. All Rights Reserved.
3  *
4  *  Use of this source code is governed by a BSD-style license
5  *  that can be found in the LICENSE file in the root of the source
6  *  tree. An additional intellectual property rights grant can be found
7  *  in the file PATENTS.  All contributing project authors may
8  *  be found in the AUTHORS file in the root of the source tree.
9  */
10 
11 /*
12  * This file contains the splitting filter functions.
13  *
14  */
15 
16 #include "rtc_base/checks.h"
17 #include "common_audio/signal_processing/include/signal_processing_library.h"
18 
19 // Maximum number of samples in a low/high-band frame.
20 enum
21 {
22     kMaxBandFrameLength = 320  // 10 ms at 64 kHz.
23 };
24 
25 // QMF filter coefficients in Q16.
26 static const uint16_t WebRtcSpl_kAllPassFilter1[3] = {6418, 36982, 57261};
27 static const uint16_t WebRtcSpl_kAllPassFilter2[3] = {21333, 49062, 63010};
28 
29 ///////////////////////////////////////////////////////////////////////////////////////////////
30 // WebRtcSpl_AllPassQMF(...)
31 //
32 // Allpass filter used by the analysis and synthesis parts of the QMF filter.
33 //
34 // Input:
35 //    - in_data             : Input data sequence (Q10)
36 //    - data_length         : Length of data sequence (>2)
37 //    - filter_coefficients : Filter coefficients (length 3, Q16)
38 //
39 // Input & Output:
40 //    - filter_state        : Filter state (length 6, Q10).
41 //
42 // Output:
43 //    - out_data            : Output data sequence (Q10), length equal to
44 //                            `data_length`
45 //
46 
WebRtcSpl_AllPassQMF(int32_t * in_data,size_t data_length,int32_t * out_data,const uint16_t * filter_coefficients,int32_t * filter_state)47 static void WebRtcSpl_AllPassQMF(int32_t* in_data,
48                                  size_t data_length,
49                                  int32_t* out_data,
50                                  const uint16_t* filter_coefficients,
51                                  int32_t* filter_state)
52 {
53     // The procedure is to filter the input with three first order all pass
54     // filters (cascade operations).
55     //
56     //         a_3 + q^-1    a_2 + q^-1    a_1 + q^-1
57     // y[n] =  -----------   -----------   -----------   x[n]
58     //         1 + a_3q^-1   1 + a_2q^-1   1 + a_1q^-1
59     //
60     // The input vector `filter_coefficients` includes these three filter
61     // coefficients. The filter state contains the in_data state, in_data[-1],
62     // followed by the out_data state, out_data[-1]. This is repeated for each
63     // cascade. The first cascade filter will filter the `in_data` and store
64     // the output in `out_data`. The second will the take the `out_data` as
65     // input and make an intermediate storage in `in_data`, to save memory. The
66     // third, and final, cascade filter operation takes the `in_data` (which is
67     // the output from the previous cascade filter) and store the output in
68     // `out_data`. Note that the input vector values are changed during the
69     // process.
70     size_t k;
71     int32_t diff;
72     // First all-pass cascade; filter from in_data to out_data.
73 
74     // Let y_i[n] indicate the output of cascade filter i (with filter
75     // coefficient a_i) at vector position n. Then the final output will be
76     // y[n] = y_3[n]
77 
78     // First loop, use the states stored in memory.
79     // "diff" should be safe from wrap around since max values are 2^25
80     // diff = (x[0] - y_1[-1])
81     diff = WebRtcSpl_SubSatW32(in_data[0], filter_state[1]);
82     // y_1[0] =  x[-1] + a_1 * (x[0] - y_1[-1])
83     out_data[0] = WEBRTC_SPL_SCALEDIFF32(filter_coefficients[0], diff, filter_state[0]);
84 
85     // For the remaining loops, use previous values.
86     for (k = 1; k < data_length; k++)
87     {
88         // diff = (x[n] - y_1[n-1])
89         diff = WebRtcSpl_SubSatW32(in_data[k], out_data[k - 1]);
90         // y_1[n] =  x[n-1] + a_1 * (x[n] - y_1[n-1])
91         out_data[k] = WEBRTC_SPL_SCALEDIFF32(filter_coefficients[0], diff, in_data[k - 1]);
92     }
93 
94     // Update states.
95     filter_state[0] = in_data[data_length - 1]; // x[N-1], becomes x[-1] next time
96     filter_state[1] = out_data[data_length - 1]; // y_1[N-1], becomes y_1[-1] next time
97 
98     // Second all-pass cascade; filter from out_data to in_data.
99     // diff = (y_1[0] - y_2[-1])
100     diff = WebRtcSpl_SubSatW32(out_data[0], filter_state[3]);
101     // y_2[0] =  y_1[-1] + a_2 * (y_1[0] - y_2[-1])
102     in_data[0] = WEBRTC_SPL_SCALEDIFF32(filter_coefficients[1], diff, filter_state[2]);
103     for (k = 1; k < data_length; k++)
104     {
105         // diff = (y_1[n] - y_2[n-1])
106         diff = WebRtcSpl_SubSatW32(out_data[k], in_data[k - 1]);
107         // y_2[0] =  y_1[-1] + a_2 * (y_1[0] - y_2[-1])
108         in_data[k] = WEBRTC_SPL_SCALEDIFF32(filter_coefficients[1], diff, out_data[k-1]);
109     }
110 
111     filter_state[2] = out_data[data_length - 1]; // y_1[N-1], becomes y_1[-1] next time
112     filter_state[3] = in_data[data_length - 1]; // y_2[N-1], becomes y_2[-1] next time
113 
114     // Third all-pass cascade; filter from in_data to out_data.
115     // diff = (y_2[0] - y[-1])
116     diff = WebRtcSpl_SubSatW32(in_data[0], filter_state[5]);
117     // y[0] =  y_2[-1] + a_3 * (y_2[0] - y[-1])
118     out_data[0] = WEBRTC_SPL_SCALEDIFF32(filter_coefficients[2], diff, filter_state[4]);
119     for (k = 1; k < data_length; k++)
120     {
121         // diff = (y_2[n] - y[n-1])
122         diff = WebRtcSpl_SubSatW32(in_data[k], out_data[k - 1]);
123         // y[n] =  y_2[n-1] + a_3 * (y_2[n] - y[n-1])
124         out_data[k] = WEBRTC_SPL_SCALEDIFF32(filter_coefficients[2], diff, in_data[k-1]);
125     }
126     filter_state[4] = in_data[data_length - 1]; // y_2[N-1], becomes y_2[-1] next time
127     filter_state[5] = out_data[data_length - 1]; // y[N-1], becomes y[-1] next time
128 }
129 
WebRtcSpl_AnalysisQMF(const int16_t * in_data,size_t in_data_length,int16_t * low_band,int16_t * high_band,int32_t * filter_state1,int32_t * filter_state2)130 void WebRtcSpl_AnalysisQMF(const int16_t* in_data, size_t in_data_length,
131                            int16_t* low_band, int16_t* high_band,
132                            int32_t* filter_state1, int32_t* filter_state2)
133 {
134     size_t i;
135     int16_t k;
136     int32_t tmp;
137     int32_t half_in1[kMaxBandFrameLength];
138     int32_t half_in2[kMaxBandFrameLength];
139     int32_t filter1[kMaxBandFrameLength];
140     int32_t filter2[kMaxBandFrameLength];
141     const size_t band_length = in_data_length / 2;
142     RTC_DCHECK_EQ(0, in_data_length % 2);
143     RTC_DCHECK_LE(band_length, kMaxBandFrameLength);
144 
145     // Split even and odd samples. Also shift them to Q10.
146     for (i = 0, k = 0; i < band_length; i++, k += 2)
147     {
148         half_in2[i] = ((int32_t)in_data[k]) * (1 << 10);
149         half_in1[i] = ((int32_t)in_data[k + 1]) * (1 << 10);
150     }
151 
152     // All pass filter even and odd samples, independently.
153     WebRtcSpl_AllPassQMF(half_in1, band_length, filter1,
154                          WebRtcSpl_kAllPassFilter1, filter_state1);
155     WebRtcSpl_AllPassQMF(half_in2, band_length, filter2,
156                          WebRtcSpl_kAllPassFilter2, filter_state2);
157 
158     // Take the sum and difference of filtered version of odd and even
159     // branches to get upper & lower band.
160     for (i = 0; i < band_length; i++)
161     {
162         tmp = (filter1[i] + filter2[i] + 1024) >> 11;
163         low_band[i] = WebRtcSpl_SatW32ToW16(tmp);
164 
165         tmp = (filter1[i] - filter2[i] + 1024) >> 11;
166         high_band[i] = WebRtcSpl_SatW32ToW16(tmp);
167     }
168 }
169 
WebRtcSpl_SynthesisQMF(const int16_t * low_band,const int16_t * high_band,size_t band_length,int16_t * out_data,int32_t * filter_state1,int32_t * filter_state2)170 void WebRtcSpl_SynthesisQMF(const int16_t* low_band, const int16_t* high_band,
171                             size_t band_length, int16_t* out_data,
172                             int32_t* filter_state1, int32_t* filter_state2)
173 {
174     int32_t tmp;
175     int32_t half_in1[kMaxBandFrameLength];
176     int32_t half_in2[kMaxBandFrameLength];
177     int32_t filter1[kMaxBandFrameLength];
178     int32_t filter2[kMaxBandFrameLength];
179     size_t i;
180     int16_t k;
181     RTC_DCHECK_LE(band_length, kMaxBandFrameLength);
182 
183     // Obtain the sum and difference channels out of upper and lower-band channels.
184     // Also shift to Q10 domain.
185     for (i = 0; i < band_length; i++)
186     {
187         tmp = (int32_t)low_band[i] + (int32_t)high_band[i];
188         half_in1[i] = tmp * (1 << 10);
189         tmp = (int32_t)low_band[i] - (int32_t)high_band[i];
190         half_in2[i] = tmp * (1 << 10);
191     }
192 
193     // all-pass filter the sum and difference channels
194     WebRtcSpl_AllPassQMF(half_in1, band_length, filter1,
195                          WebRtcSpl_kAllPassFilter2, filter_state1);
196     WebRtcSpl_AllPassQMF(half_in2, band_length, filter2,
197                          WebRtcSpl_kAllPassFilter1, filter_state2);
198 
199     // The filtered signals are even and odd samples of the output. Combine
200     // them. The signals are Q10 should shift them back to Q0 and take care of
201     // saturation.
202     for (i = 0, k = 0; i < band_length; i++)
203     {
204         tmp = (filter2[i] + 512) >> 10;
205         out_data[k++] = WebRtcSpl_SatW32ToW16(tmp);
206 
207         tmp = (filter1[i] + 512) >> 10;
208         out_data[k++] = WebRtcSpl_SatW32ToW16(tmp);
209     }
210 
211 }
212