xref: /btstack/src/classic/btstack_cvsd_plc.c (revision f9f2075ceac5e9dc08e9abea437e43d733a3a0ea)
1 /*
2  * Copyright (C) 2016 BlueKitchen GmbH
3  *
4  * Redistribution and use in source and binary forms, with or without
5  * modification, are permitted provided that the following conditions
6  * are met:
7  *
8  * 1. Redistributions of source code must retain the above copyright
9  *    notice, this list of conditions and the following disclaimer.
10  * 2. Redistributions in binary form must reproduce the above copyright
11  *    notice, this list of conditions and the following disclaimer in the
12  *    documentation and/or other materials provided with the distribution.
13  * 3. Neither the name of the copyright holders nor the names of
14  *    contributors may be used to endorse or promote products derived
15  *    from this software without specific prior written permission.
16  * 4. Any redistribution, use, or modification is done solely for
17  *    personal benefit and not for any commercial purpose or for
18  *    monetary gain.
19  *
20  * THIS SOFTWARE IS PROVIDED BY BLUEKITCHEN GMBH AND CONTRIBUTORS
21  * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23  * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL MATTHIAS
24  * RINGWALD OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26  * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
27  * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
28  * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
29  * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
30  * THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
31  * SUCH DAMAGE.
32  *
33  * Please inquire about commercial licensing options at
34  * [email protected]
35  *
36  */
37 
38 /*
39  * btstack_sbc_plc.c
40  *
41  */
42 
43 
44 #include <stdint.h>
45 #include <stdio.h>
46 #include <stdlib.h>
47 #include <string.h>
48 
49 #include "btstack_cvsd_plc.h"
50 #include "btstack_debug.h"
51 
52 #define SAMPLE_FORMAT int16_t
53 
54 static float rcos[CVSD_OLAL] = {
55     0.99148655f,0.96623611f,0.92510857f,0.86950446f,
56     0.80131732f,0.72286918f,0.63683150f,0.54613418f,
57     0.45386582f,0.36316850f,0.27713082f,0.19868268f,
58     0.13049554f,0.07489143f,0.03376389f,0.00851345f};
59 
60 // taken from http://www.codeproject.com/Articles/69941/Best-Square-Root-Method-Algorithm-Function-Precisi
61 // Algorithm: Babylonian Method + some manipulations on IEEE 32 bit floating point representation
62 static float sqrt3(const float x){
63     union {
64         int i;
65         float x;
66     } u;
67     u.x = x;
68     u.i = (1<<29) + (u.i >> 1) - (1<<22);
69 
70     // Two Babylonian Steps (simplified from:)
71     // u.x = 0.5f * (u.x + x/u.x);
72     // u.x = 0.5f * (u.x + x/u.x);
73     u.x =       u.x + x/u.x;
74     u.x = 0.25f*u.x + x/u.x;
75 
76     return u.x;
77 }
78 
79 static float absolute(float x){
80      if (x < 0) x = -x;
81      return x;
82 }
83 
84 static float CrossCorrelation(SAMPLE_FORMAT *x, SAMPLE_FORMAT *y){
85     float num = 0;
86     float den = 0;
87     float x2 = 0;
88     float y2 = 0;
89     int   m;
90     for (m=0;m<CVSD_M;m++){
91         num+=((float)x[m])*y[m];
92         x2+=((float)x[m])*x[m];
93         y2+=((float)y[m])*y[m];
94     }
95     den = (float)sqrt3(x2*y2);
96     return num/den;
97 }
98 
99 static int PatternMatch(SAMPLE_FORMAT *y){
100     float maxCn = -999999.0;  // large negative number
101     int   bestmatch = 0;
102     float Cn;
103     int   n;
104     for (n=0;n<CVSD_N;n++){
105         Cn = CrossCorrelation(&y[CVSD_LHIST-CVSD_M], &y[n]);
106         if (Cn>maxCn){
107             bestmatch=n;
108             maxCn = Cn;
109         }
110     }
111     return bestmatch;
112 }
113 
114 static float AmplitudeMatch(SAMPLE_FORMAT *y, SAMPLE_FORMAT bestmatch) {
115     int   i;
116     float sumx = 0;
117     float sumy = 0.000001f;
118     float sf;
119 
120     for (i=0;i<CVSD_FS;i++){
121         sumx += absolute(y[CVSD_LHIST-CVSD_FS+i]);
122         sumy += absolute(y[bestmatch+i]);
123     }
124     sf = sumx/sumy;
125     // This is not in the paper, but limit the scaling factor to something reasonable to avoid creating artifacts
126     if (sf<0.75f) sf=0.75f;
127     if (sf>1.2f) sf=1.2f;
128     return sf;
129 }
130 
131 static SAMPLE_FORMAT crop_sample(float val){
132     float croped_val = val;
133     if (croped_val > 32767.0)  croped_val= 32767.0;
134     if (croped_val < -32768.0) croped_val=-32768.0;
135     return (SAMPLE_FORMAT) croped_val;
136 }
137 
138 void btstack_cvsd_plc_init(btstack_cvsd_plc_state_t *plc_state){
139     memset(plc_state, 0, sizeof(btstack_cvsd_plc_state_t));
140 }
141 
142 void btstack_cvsd_plc_bad_frame(btstack_cvsd_plc_state_t *plc_state, SAMPLE_FORMAT *out){
143     float val;
144     int   i = 0;
145     float sf = 1;
146     plc_state->nbf++;
147 
148     if (plc_state->nbf==1){
149         // Perform pattern matching to find where to replicate
150         plc_state->bestlag = PatternMatch(plc_state->hist);
151         // the replication begins after the template match
152         plc_state->bestlag += CVSD_M;
153 
154         // Compute Scale Factor to Match Amplitude of Substitution Packet to that of Preceding Packet
155         sf = AmplitudeMatch(plc_state->hist, plc_state->bestlag);
156         for (i=0;i<CVSD_OLAL;i++){
157             val = sf*plc_state->hist[plc_state->bestlag+i];
158             plc_state->hist[CVSD_LHIST+i] = crop_sample(val);
159         }
160 
161         for (;i<CVSD_FS;i++){
162             val = sf*plc_state->hist[plc_state->bestlag+i];
163             plc_state->hist[CVSD_LHIST+i] = crop_sample(val);
164         }
165 
166         for (;i<CVSD_FS+CVSD_OLAL;i++){
167             float left  = sf*plc_state->hist[plc_state->bestlag+i];
168             float right = plc_state->hist[plc_state->bestlag+i];
169             val = left*rcos[i-CVSD_FS] + right*rcos[CVSD_OLAL-1-i+CVSD_FS];
170             plc_state->hist[CVSD_LHIST+i] = crop_sample(val);
171         }
172 
173         for (;i<CVSD_FS+CVSD_RT+CVSD_OLAL;i++){
174             plc_state->hist[CVSD_LHIST+i] = plc_state->hist[plc_state->bestlag+i];
175         }
176     } else {
177         for (;i<CVSD_FS+CVSD_RT+CVSD_OLAL;i++){
178             plc_state->hist[CVSD_LHIST+i] = plc_state->hist[plc_state->bestlag+i];
179         }
180     }
181     for (i=0;i<CVSD_FS;i++){
182         out[i] = plc_state->hist[CVSD_LHIST+i];
183     }
184 
185    // shift the history buffer
186    for (i=0;i<CVSD_LHIST+CVSD_RT+CVSD_OLAL;i++){
187         plc_state->hist[i] = plc_state->hist[i+CVSD_FS];
188     }
189 }
190 
191 void btstack_cvsd_plc_good_frame(btstack_cvsd_plc_state_t *plc_state, SAMPLE_FORMAT *in, SAMPLE_FORMAT *out){
192     float val;
193     int i = 0;
194     if (plc_state->nbf>0){
195         for (i=0;i<CVSD_RT;i++){
196             out[i] = plc_state->hist[CVSD_LHIST+i];
197         }
198 
199         for (i=CVSD_RT;i<CVSD_RT+CVSD_OLAL;i++){
200             float left  = plc_state->hist[CVSD_LHIST+i];
201             float right = in[i];
202             val = left * rcos[i-CVSD_RT] + right *rcos[CVSD_OLAL+CVSD_RT-1-i];
203             out[i] = (SAMPLE_FORMAT)val;
204         }
205     }
206 
207     for (;i<CVSD_FS;i++){
208         out[i] = in[i];
209     }
210     // Copy the output to the history buffer
211     for (i=0;i<CVSD_FS;i++){
212         plc_state->hist[CVSD_LHIST+i] = out[i];
213     }
214     // shift the history buffer
215     for (i=0;i<CVSD_LHIST;i++){
216         plc_state->hist[i] = plc_state->hist[i+CVSD_FS];
217     }
218     plc_state->nbf=0;
219 }
220 
221 static int count_equal_samples(SAMPLE_FORMAT * packet, uint16_t size){
222     int count = 0;
223     int temp_count = 1;
224     int i;
225     for (i = 0; i < size-1; i++){
226         if (packet[i] == packet[i+1]){
227             temp_count++;
228             continue;
229         }
230         if (count < temp_count){
231             count = temp_count;
232         }
233         temp_count = 1;
234     }
235     if (temp_count > count + 1){
236         count = temp_count;
237     }
238     return count;
239 }
240 
241 // @assumption frame len 24 samples
242 static int bad_frame(SAMPLE_FORMAT * frame, uint16_t size){
243     return count_equal_samples(frame, size) > 20;
244 }
245 
246 void btstack_cvsd_plc_process_data(btstack_cvsd_plc_state_t * state, SAMPLE_FORMAT * in, uint16_t size, SAMPLE_FORMAT * out){
247     if (size != 24){
248         log_error("btstack_cvsd_plc_process_data: audio frame size is incorrect. Expected %d, got %d", CVSD_FS, size);
249         return;
250     }
251     state->frame_count++;
252     if (bad_frame(in,size)){
253         memcpy(out, in, size);
254         if (state->good_frames_nr > CVSD_LHIST/CVSD_FS){
255             btstack_cvsd_plc_bad_frame(state, out);
256             state->bad_frames_nr++;
257         } else {
258             memset(out, 0, CVSD_FS);
259         }
260     } else {
261         btstack_cvsd_plc_good_frame(state, in, out);
262         state->good_frames_nr++;
263         if (state->good_frames_nr == 1){
264             log_info("First good frame at index %d\n", state->frame_count-1);
265         }
266     }
267 }
268 
269 void btstack_cvsd_dump_statistics(btstack_cvsd_plc_state_t * state){
270     log_info("Good frames: %d\n", state->good_frames_nr);
271     log_info("Bad frames: %d\n", state->bad_frames_nr);
272 }
273