SILK update with LBRR and some bugfixes
[opus.git] / test / signalCompare.c
1 /***********************************************************************\r
2 Copyright (c) 2006-2011, Skype Limited. All rights reserved. \r
3 Redistribution and use in source and binary forms, with or without \r
4 modification, (subject to the limitations in the disclaimer below) \r
5 are permitted provided that the following conditions are met:\r
6 - Redistributions of source code must retain the above copyright notice,\r
7 this list of conditions and the following disclaimer.\r
8 - Redistributions in binary form must reproduce the above copyright \r
9 notice, this list of conditions and the following disclaimer in the \r
10 documentation and/or other materials provided with the distribution.\r
11 - Neither the name of Skype Limited, nor the names of specific \r
12 contributors, may be used to endorse or promote products derived from \r
13 this software without specific prior written permission.\r
14 NO EXPRESS OR IMPLIED LICENSES TO ANY PARTY'S PATENT RIGHTS ARE GRANTED \r
15 BY THIS LICENSE. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND \r
16 CONTRIBUTORS ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,\r
17 BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND \r
18 FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE \r
19 COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, \r
20 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT\r
21 NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF \r
22 USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON \r
23 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT \r
24 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE \r
25 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.\r
26 ***********************************************************************/\r
27 \r
28 /*\r
29 * Compare two audio signals and compute weighted SNR difference\r
30 */\r
31 \r
32 #ifdef _WIN32\r
33 #define _CRT_SECURE_NO_DEPRECATE    1\r
34 #endif\r
35 \r
36 #include <stdio.h>\r
37 #include <stdlib.h>\r
38 #include <string.h>\r
39 #include <math.h>\r
40 \r
41 #include "SKP_Silk_SigProc_FIX.h"\r
42 \r
43 #define FRAME_LENGTH_MS         10\r
44 #define WIN_LENGTH_MS           20\r
45 #define BW_EXPANSION            0.7f\r
46 \r
47 #define MAX_FS_KHZ              24\r
48 #define LPC_ORDER               10\r
49 \r
50 #ifdef  __cplusplus\r
51 extern "C"\r
52 {\r
53 #endif\r
54 /* Internally used functions */\r
55 void Autocorrelation( \r
56     SKP_float *results,                 /* o    result (length correlationCount)            */\r
57     const SKP_float *inputData,         /* i    input data to correlate                     */\r
58     SKP_int inputDataSize,              /* i    length of input                             */\r
59     SKP_int correlationCount            /* i    number of correlation taps to compute       */\r
60 );\r
61 \r
62 /* inner product of two SKP_float arrays, with result as double */\r
63 double Inner_product( \r
64     const SKP_float     *data1, \r
65     const SKP_float     *data2, \r
66     SKP_int             dataSize\r
67 );\r
68 /* Solve the normal equations using the Levinson-Durbin recursion */\r
69 SKP_float Levinsondurbin(               /* O    prediction error energy                     */\r
70     SKP_float       A[],                /* O    prediction coefficients [order]             */\r
71     const SKP_float corr[],             /* I    input auto-correlations [order + 1]         */\r
72     const SKP_int   order               /* I    prediction order                            */\r
73 );\r
74 \r
75 /* Chirp (bw expand) LP AR filter */\r
76 void Bwexpander( \r
77     SKP_float *ar,                      /* io   AR filter to be expanded (without leading 1)    */\r
78     const SKP_int d,                    /* i    length of ar                                    */\r
79     const SKP_float chirp               /* i    chirp factor (typically in range (0..1) )       */\r
80 );\r
81 \r
82 #ifdef  __cplusplus\r
83 }\r
84 #endif\r
85 \r
86 static void print_usage(char* argv[]) {\r
87     printf("\nusage: %s ref.pcm test.pcm [settings]\n", argv[ 0 ]);\r
88     printf("\nref.pcm       : Reference file");\r
89     printf("\ntest.pcm      : File to be tested, should be of same length as ref.pcm");\r
90     printf("\n   settings:");\r
91     printf("\n-diff         : Only determine bit-exactness");\r
92     printf("\n-fs <Hz>      : Sampling rate in Hz, max: %d; default: 24000", MAX_FS_KHZ * 1000 );\r
93     printf("\n");\r
94 }\r
95 \r
96 #ifdef SKP_MACRO_COUNT\r
97     varDefine /* Define and reset global counters */\r
98 #endif\r
99 \r
100 int main(int argc, char* argv[])\r
101 {\r
102     SKP_int   args, n, i, counterRef, counterTest;\r
103     char      testInFileName[150], refInFileName[150];\r
104     FILE      *refInFile, *testInFile;\r
105     SKP_int   nFrames = 0, isUnequal = 0;\r
106     SKP_int   diff = 0, Fs_kHz;\r
107     SKP_int32 Fs_Hz = 24000;\r
108     SKP_float c, refWhtnd, testWhtnd, refNrg, diffNrg;\r
109     double    SNR = 0.0;\r
110     SKP_int16 refIn[WIN_LENGTH_MS * MAX_FS_KHZ], testIn[WIN_LENGTH_MS * MAX_FS_KHZ];\r
111     SKP_float refWin[WIN_LENGTH_MS * MAX_FS_KHZ], testWin[WIN_LENGTH_MS * MAX_FS_KHZ];\r
112     SKP_float autoCorr[LPC_ORDER + 1], LPC_Coef[LPC_ORDER];\r
113 \r
114     if (argc < 3) {\r
115         print_usage(argv);\r
116         exit(0);\r
117     } \r
118 \r
119     /* get arguments */\r
120     args = 1;\r
121     strcpy(refInFileName, argv[args]);\r
122     args++;\r
123     strcpy(testInFileName, argv[args]);\r
124     args++;\r
125     while(args < argc ) {\r
126         if (SKP_STR_CASEINSENSITIVE_COMPARE(argv[args], "-diff") == 0) {\r
127             diff = 1;\r
128             args++;\r
129         }else if (SKP_STR_CASEINSENSITIVE_COMPARE(argv[args], "-fs") == 0) {\r
130             sscanf(argv[args+1], "%d", &Fs_Hz);\r
131             args += 2;\r
132         } else {\r
133             printf("Error: unrecognized setting: %s\n\n", argv[args]);\r
134             print_usage(argv);\r
135             exit(0);\r
136         }\r
137     }\r
138 \r
139     Fs_kHz = SKP_DIV32_16( Fs_Hz, 1000 );\r
140 \r
141     if( Fs_kHz > MAX_FS_KHZ ) {\r
142         printf("Error: sampling rate too high: %d\n\n", Fs_kHz);\r
143         print_usage(argv);\r
144         exit(0);\r
145     }\r
146 \r
147     printf("Reference:  %s\n", refInFileName);\r
148     //printf("Test:       %s\n", testInFileName);\r
149 \r
150     /* open files */\r
151     refInFile = fopen(refInFileName, "rb");\r
152     if (refInFile==NULL) {\r
153         printf("Error: could not open input file %s\n", refInFileName);\r
154         exit(0);\r
155     } \r
156     testInFile = fopen(testInFileName, "rb");\r
157     if (testInFile==NULL) {\r
158         printf("Error: could not open input file %s\n", testInFileName);\r
159         exit(0);\r
160     }\r
161 \r
162     SKP_memset( refIn,  0, sizeof(refIn) );\r
163     SKP_memset( testIn, 0, sizeof(testIn) );\r
164 \r
165     while(1) {\r
166         /* Read inputs */\r
167         counterRef  = (SKP_int)fread(&refIn[(WIN_LENGTH_MS - FRAME_LENGTH_MS) * Fs_kHz], \r
168             sizeof(SKP_int16), FRAME_LENGTH_MS * Fs_kHz, refInFile);\r
169         counterTest = (SKP_int)fread(&testIn[(WIN_LENGTH_MS - FRAME_LENGTH_MS) * Fs_kHz], \r
170             sizeof(SKP_int16), FRAME_LENGTH_MS * Fs_kHz, testInFile);\r
171         if(counterRef != FRAME_LENGTH_MS * Fs_kHz || counterTest != FRAME_LENGTH_MS * Fs_kHz){\r
172             break;\r
173         }\r
174 \r
175         /* test for bit-exactness */\r
176         for( n = 0; n < FRAME_LENGTH_MS * Fs_kHz; n++ ) {\r
177             if( refIn[(WIN_LENGTH_MS - FRAME_LENGTH_MS) * Fs_kHz + n] != \r
178                 testIn[(WIN_LENGTH_MS - FRAME_LENGTH_MS) * Fs_kHz + n] ) {\r
179                     isUnequal = 1;\r
180                     break;\r
181             }\r
182         }\r
183 \r
184         /* apply sine window */\r
185         for( n = 0; n < WIN_LENGTH_MS * Fs_kHz; n++ ) {\r
186             c = (SKP_float)sin( 3.14159265 * (n + 1) / (WIN_LENGTH_MS * Fs_kHz + 1) );\r
187             refWin[n]  = refIn[n]  * c;\r
188             testWin[n] = testIn[n] * c;\r
189         }\r
190 \r
191         /* LPC analysis on reference signal */\r
192 \r
193         /* Calculate auto correlation */\r
194         Autocorrelation(autoCorr, refWin, WIN_LENGTH_MS * Fs_kHz, LPC_ORDER + 1);\r
195 \r
196         /* Add white noise */\r
197         autoCorr[ 0 ] += autoCorr[ 0 ] * 1e-6f + 1.0f; \r
198 \r
199         /* Convert correlations to prediction coefficients */\r
200         Levinsondurbin(LPC_Coef, autoCorr, LPC_ORDER);\r
201 \r
202         /* Bandwdith expansion */\r
203         Bwexpander(LPC_Coef, LPC_ORDER, BW_EXPANSION);\r
204 \r
205         /* Filter both signals */\r
206         refNrg = 1.0f;\r
207         diffNrg = 1e-10f;\r
208         for( n = (WIN_LENGTH_MS - FRAME_LENGTH_MS) / 2 * Fs_kHz; \r
209              n < (WIN_LENGTH_MS + FRAME_LENGTH_MS) / 2 * Fs_kHz; n++ ) {\r
210                 refWhtnd = refIn[n];\r
211                 testWhtnd = testIn[n];\r
212                 for( i = 0; i < LPC_ORDER; i++ ) {\r
213                     refWhtnd  -= LPC_Coef[ i ] * refIn[n - i - 1];\r
214                     testWhtnd -= LPC_Coef[ i ] * testIn[n - i - 1];\r
215                 }\r
216                 refNrg += refWhtnd * refWhtnd;\r
217                 diffNrg += (refWhtnd - testWhtnd) * (refWhtnd - testWhtnd);\r
218         }\r
219 \r
220         /* weighted SNR */\r
221         if( refNrg > FRAME_LENGTH_MS * Fs_kHz ) {\r
222             SNR += 10.0 * log10( refNrg / diffNrg );\r
223             nFrames++;\r
224         }\r
225 \r
226         /* Update Buffer */\r
227         SKP_memmove( refIn,  &refIn[FRAME_LENGTH_MS * Fs_kHz],  (WIN_LENGTH_MS - FRAME_LENGTH_MS) * Fs_kHz * sizeof(SKP_int16));\r
228         SKP_memmove( testIn, &testIn[FRAME_LENGTH_MS * Fs_kHz], (WIN_LENGTH_MS - FRAME_LENGTH_MS) * Fs_kHz * sizeof(SKP_int16));\r
229     }\r
230 \r
231     if( diff ) {\r
232         if( isUnequal ) {\r
233             printf("Signals     DIFFER\n");\r
234         } else {\r
235             if(counterRef != counterTest){\r
236                 printf("Warning: signals differ in length\n");\r
237             }\r
238             printf("Signals     BIT-EXACT\n");\r
239         }\r
240     } else {\r
241         if( nFrames == 0 ) {\r
242             printf("At least one signal too short or not loud enough\n");\r
243             exit(0);\r
244         }\r
245         if(counterRef != counterTest){\r
246             printf("Warning: signals differ in length\n");\r
247         }\r
248         if( isUnequal == 0 ) {\r
249             printf("Signals     BIT-EXACT\n");\r
250         } else {\r
251             printf("AVERAGE WEIGHTED SNR: %4.1f dB\n", SNR / nFrames);\r
252         }\r
253     }\r
254     printf("\n");\r
255 \r
256     /* Close Files */\r
257     fclose(refInFile);\r
258     fclose(testInFile);\r
259 \r
260     return 0;\r
261 }\r
262 \r
263 /* compute autocorrelation */\r
264 void Autocorrelation( \r
265     SKP_float *results,                 /* o    result (length correlationCount)            */\r
266     const SKP_float *inputData,         /* i    input data to correlate                     */\r
267     SKP_int inputDataSize,              /* i    length of input                             */\r
268     SKP_int correlationCount            /* i    number of correlation taps to compute       */\r
269 )\r
270 {\r
271     SKP_int i;\r
272 \r
273     if (correlationCount > inputDataSize) {\r
274         correlationCount = inputDataSize;\r
275     }\r
276 \r
277     for( i = 0; i < correlationCount; i++ ) {\r
278         results[ i ] =  (SKP_float)Inner_product( inputData, inputData + i, inputDataSize - i );\r
279     }\r
280 }\r
281 \r
282 /* inner product of two SKP_float arrays, with result as double */\r
283 double Inner_product( \r
284     const SKP_float     *data1, \r
285     const SKP_float     *data2, \r
286     SKP_int             dataSize\r
287 )\r
288 {\r
289     SKP_int  i, dataSize4;\r
290     double   result;\r
291 \r
292     /* 4x unrolled loop */\r
293     result = 0.0f;\r
294     dataSize4 = dataSize & 0xFFFC;\r
295     for( i = 0; i < dataSize4; i += 4 ) {\r
296         result += data1[ i + 0 ] * data2[ i + 0 ] + \r
297                   data1[ i + 1 ] * data2[ i + 1 ] +\r
298                   data1[ i + 2 ] * data2[ i + 2 ] +\r
299                   data1[ i + 3 ] * data2[ i + 3 ];\r
300     }\r
301 \r
302     /* add any remaining products */\r
303     for( ; i < dataSize; i++ ) {\r
304         result += data1[ i ] * data2[ i ];\r
305     }\r
306 \r
307     return result;\r
308 }\r
309 \r
310 /* Solve the normal equations using the Levinson-Durbin recursion */\r
311 SKP_float Levinsondurbin(               /* O    prediction error energy                     */\r
312     SKP_float       A[],                /* O    prediction coefficients [order]             */\r
313     const SKP_float corr[],             /* I    input auto-correlations [order + 1]         */\r
314     const SKP_int   order               /* I    prediction order                            */\r
315 )\r
316 {\r
317     SKP_int   i, mHalf, m;\r
318     SKP_float min_nrg, nrg, t, km, Atmp1, Atmp2;\r
319     \r
320     min_nrg = 1e-12f * corr[ 0 ] + 1e-9f;\r
321     nrg = corr[ 0 ];\r
322     nrg = SKP_max( min_nrg, nrg );\r
323     A[ 0 ] = corr[ 1 ] / nrg;\r
324     nrg -= A[ 0 ] * corr[ 1 ];\r
325     nrg = SKP_max(min_nrg, nrg);\r
326 \r
327     for( m = 1; m < order; m++ )\r
328     {\r
329         t = corr[ m + 1 ];\r
330         for( i = 0; i < m; i++ ) {\r
331             t -= A[ i ] * corr[ m - i ];\r
332         }\r
333 \r
334         /* reflection coefficient */\r
335         km = t / nrg;\r
336 \r
337         /* residual energy */\r
338         nrg -= km * t;\r
339         nrg = SKP_max(min_nrg, nrg);\r
340 \r
341         mHalf = m >> 1;\r
342         for( i = 0; i < mHalf; i++ ) {\r
343             Atmp1 = A[ i ];\r
344             Atmp2 = A[ m - i - 1 ];\r
345             A[ m - i - 1 ] -= km * Atmp1;\r
346             A[ i ]         -= km * Atmp2;\r
347         }\r
348         if( m & 1 ) {\r
349             A[ mHalf ]     -= km * A[ mHalf ];\r
350         }\r
351         A[ m ] = km;\r
352     }\r
353 \r
354     /* return the residual energy */\r
355     return nrg;\r
356 }\r
357 \r
358 /* Chirp (bw expand) LP AR filter */\r
359 void Bwexpander( \r
360     SKP_float *ar,                      /* io   AR filter to be expanded (without leading 1)    */\r
361     const SKP_int d,                    /* i    length of ar                                    */\r
362     const SKP_float chirp               /* i    chirp factor (typically in range (0..1) )       */\r
363 )\r
364 {\r
365     SKP_int   i;\r
366     SKP_float cfac = chirp;\r
367 \r
368     for( i = 0; i < d - 1; i++ ) {\r
369         ar[ i ] *=  cfac;\r
370         cfac    *=  chirp;\r
371     }\r
372     ar[ d - 1 ] *=  cfac;\r
373 }\r