Initial Skype commit taken from FreeSwitch, which got it from the IETF draft.
[opus.git] / test / signalCompare.c
1 /***********************************************************************\r
2 Copyright (c) 2006-2010, 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 <kHz>     : Sampling rate in kHz, max: %d; default: 24", MAX_FS_KHZ);\r
93         printf("\n");\r
94 }\r
95 \r
96 int main(int argc, char* argv[])\r
97 {\r
98         SKP_int   args, n, i, counterRef, counterTest;\r
99         char      testInFileName[150], refInFileName[150];\r
100         FILE      *refInFile, *testInFile;\r
101         SKP_int   nFrames = 0, isUnequal = 0;\r
102         SKP_int   diff = 0, Fs_kHz = 24;                /* defaults */\r
103         SKP_float c, refWhtnd, testWhtnd, refNrg, diffNrg;\r
104         double    SNR = 0.0;\r
105         SKP_int16 refIn[WIN_LENGTH_MS * MAX_FS_KHZ], testIn[WIN_LENGTH_MS * MAX_FS_KHZ];\r
106         SKP_float refWin[WIN_LENGTH_MS * MAX_FS_KHZ], testWin[WIN_LENGTH_MS * MAX_FS_KHZ];\r
107         SKP_float autoCorr[LPC_ORDER + 1], LPC_Coef[LPC_ORDER];\r
108 \r
109         if (argc < 3) {\r
110                 print_usage(argv);\r
111                 exit(0);\r
112         } \r
113 \r
114         /* get arguments */\r
115         args = 1;\r
116         strcpy(refInFileName, argv[args]);\r
117         args++;\r
118         strcpy(testInFileName, argv[args]);\r
119         args++;\r
120         while(args < argc ) {\r
121                 if (SKP_STR_CASEINSENSITIVE_COMPARE(argv[args], "-diff") == 0) {\r
122                         diff = 1;\r
123                         args++;\r
124                 }else if (SKP_STR_CASEINSENSITIVE_COMPARE(argv[args], "-fs") == 0) {\r
125                         sscanf(argv[args+1], "%d", &Fs_kHz);\r
126                         args += 2;\r
127                 } else {\r
128                         printf("Error: unrecognized setting: %s\n\n", argv[args]);\r
129                         print_usage(argv);\r
130                         exit(0);\r
131                 }\r
132         }\r
133 \r
134         if( Fs_kHz > MAX_FS_KHZ ) {\r
135                 printf("Error: sampling rate too high: %d\n\n", Fs_kHz);\r
136                 print_usage(argv);\r
137                 exit(0);\r
138         }\r
139 \r
140         printf("Reference:  %s\n", refInFileName);\r
141         //printf("Test:       %s\n", testInFileName);\r
142 \r
143         /* open files */\r
144         refInFile = fopen(refInFileName, "rb");\r
145         if (refInFile==NULL) {\r
146                 printf("Error: could not open input file %s\n", refInFileName);\r
147                 exit(0);\r
148         } \r
149         testInFile = fopen(testInFileName, "rb");\r
150         if (testInFile==NULL) {\r
151                 printf("Error: could not open input file %s\n", testInFileName);\r
152                 exit(0);\r
153         }\r
154 \r
155         SKP_memset( refIn,  0, sizeof(refIn) );\r
156         SKP_memset( testIn, 0, sizeof(testIn) );\r
157 \r
158         while(1) {\r
159                 /* Read inputs */\r
160                 counterRef  = (SKP_int)fread(&refIn[(WIN_LENGTH_MS - FRAME_LENGTH_MS) * Fs_kHz], \r
161                         sizeof(SKP_int16), FRAME_LENGTH_MS * Fs_kHz, refInFile);\r
162                 counterTest = (SKP_int)fread(&testIn[(WIN_LENGTH_MS - FRAME_LENGTH_MS) * Fs_kHz], \r
163                         sizeof(SKP_int16), FRAME_LENGTH_MS * Fs_kHz, testInFile);\r
164                 if(counterRef != FRAME_LENGTH_MS * Fs_kHz || counterTest != FRAME_LENGTH_MS * Fs_kHz){\r
165                         break;\r
166                 }\r
167 \r
168                 /* test for bit-exactness */\r
169                 for( n = 0; n < FRAME_LENGTH_MS * Fs_kHz; n++ ) {\r
170                         if( refIn[(WIN_LENGTH_MS - FRAME_LENGTH_MS) * Fs_kHz + n] != \r
171                                 testIn[(WIN_LENGTH_MS - FRAME_LENGTH_MS) * Fs_kHz + n] ) {\r
172                                         isUnequal = 1;\r
173                                         break;\r
174                         }\r
175                 }\r
176 \r
177                 /* apply sine window */\r
178                 for( n = 0; n < WIN_LENGTH_MS * Fs_kHz; n++ ) {\r
179                         c = (SKP_float)sin( 3.14159265 * (n + 1) / (WIN_LENGTH_MS * Fs_kHz + 1) );\r
180                         refWin[n]  = refIn[n]  * c;\r
181                         testWin[n] = testIn[n] * c;\r
182                 }\r
183 \r
184                 /* LPC analysis on reference signal */\r
185 \r
186                 /* Calculate auto correlation */\r
187                 Autocorrelation(autoCorr, refWin, WIN_LENGTH_MS * Fs_kHz, LPC_ORDER + 1);\r
188 \r
189                 /* Add white noise */\r
190                 autoCorr[ 0 ] += autoCorr[ 0 ] * 1e-6f + 1.0f; \r
191 \r
192                 /* Convert correlations to prediction coefficients */\r
193                 Levinsondurbin(LPC_Coef, autoCorr, LPC_ORDER);\r
194 \r
195                 /* Bandwdith expansion */\r
196                 Bwexpander(LPC_Coef, LPC_ORDER, BW_EXPANSION);\r
197 \r
198                 /* Filter both signals */\r
199                 refNrg = 1.0f;\r
200                 diffNrg = 1e-10f;\r
201                 for( n = (WIN_LENGTH_MS - FRAME_LENGTH_MS) / 2 * Fs_kHz; \r
202                          n < (WIN_LENGTH_MS + FRAME_LENGTH_MS) / 2 * Fs_kHz; n++ ) {\r
203                                 refWhtnd = refIn[n];\r
204                                 testWhtnd = testIn[n];\r
205                                 for( i = 0; i < LPC_ORDER; i++ ) {\r
206                                         refWhtnd  -= LPC_Coef[ i ] * refIn[n - i - 1];\r
207                                         testWhtnd -= LPC_Coef[ i ] * testIn[n - i - 1];\r
208                                 }\r
209                                 refNrg += refWhtnd * refWhtnd;\r
210                                 diffNrg += (refWhtnd - testWhtnd) * (refWhtnd - testWhtnd);\r
211                 }\r
212 \r
213                 /* weighted SNR */\r
214                 if( refNrg > FRAME_LENGTH_MS * Fs_kHz ) {\r
215                         SNR += 10.0 * log10( refNrg / diffNrg );\r
216                         nFrames++;\r
217                 }\r
218 \r
219                 /* Update Buffer */\r
220                 SKP_memmove( refIn,  &refIn[FRAME_LENGTH_MS * Fs_kHz],  (WIN_LENGTH_MS - FRAME_LENGTH_MS) * Fs_kHz * sizeof(SKP_int16));\r
221                 SKP_memmove( testIn, &testIn[FRAME_LENGTH_MS * Fs_kHz], (WIN_LENGTH_MS - FRAME_LENGTH_MS) * Fs_kHz * sizeof(SKP_int16));\r
222         }\r
223 \r
224         if( diff ) {\r
225                 if( isUnequal ) {\r
226                         printf("Signals     DIFFER\n");\r
227                 } else {\r
228                         if(counterRef != counterTest){\r
229                                 printf("Warning: signals differ in length\n");\r
230                         }\r
231                         printf("Signals     BIT-EXACT\n");\r
232                 }\r
233         } else {\r
234                 if( nFrames == 0 ) {\r
235                         printf("At least one signal too short or not loud enough\n");\r
236                         exit(0);\r
237                 }\r
238                 if(counterRef != counterTest){\r
239                         printf("Warning: signals differ in length\n");\r
240                 }\r
241                 if( isUnequal == 0 ) {\r
242                         printf("Signals     BIT-EXACT\n");\r
243                 } else {\r
244                         printf("AVERAGE WEIGHTED SNR: %4.1f dB\n", SNR / nFrames);\r
245                 }\r
246         }\r
247         printf("\n");\r
248 \r
249         /* Close Files */\r
250         fclose(refInFile);\r
251         fclose(testInFile);\r
252 \r
253         return 0;\r
254 }\r
255 \r
256 /* compute autocorrelation */\r
257 void Autocorrelation( \r
258         SKP_float *results,                                     /* o    result (length correlationCount)                        */\r
259         const SKP_float *inputData,                     /* i    input data to correlate                                         */\r
260         SKP_int inputDataSize,                          /* i    length of input                                                         */\r
261         SKP_int correlationCount                        /* i    number of correlation taps to compute           */\r
262 )\r
263 {\r
264         SKP_int i;\r
265 \r
266         if (correlationCount > inputDataSize) {\r
267                 correlationCount = inputDataSize;\r
268         }\r
269 \r
270         for( i = 0; i < correlationCount; i++ ) {\r
271                 results[ i ] =  (SKP_float)Inner_product( inputData, inputData + i, inputDataSize - i );\r
272         }\r
273 }\r
274 \r
275 /* inner product of two SKP_float arrays, with result as double */\r
276 double Inner_product( \r
277         const SKP_float         *data1, \r
278         const SKP_float         *data2, \r
279         SKP_int                         dataSize\r
280 )\r
281 {\r
282         SKP_int  i, dataSize4;\r
283         double   result;\r
284 \r
285         /* 4x unrolled loop */\r
286         result = 0.0f;\r
287         dataSize4 = dataSize & 0xFFFC;\r
288         for( i = 0; i < dataSize4; i += 4 ) {\r
289                 result += data1[ i + 0 ] * data2[ i + 0 ] + \r
290                                   data1[ i + 1 ] * data2[ i + 1 ] +\r
291                                   data1[ i + 2 ] * data2[ i + 2 ] +\r
292                                   data1[ i + 3 ] * data2[ i + 3 ];\r
293         }\r
294 \r
295         /* add any remaining products */\r
296         for( ; i < dataSize; i++ ) {\r
297                 result += data1[ i ] * data2[ i ];\r
298         }\r
299 \r
300         return result;\r
301 }\r
302 \r
303 /* Solve the normal equations using the Levinson-Durbin recursion */\r
304 SKP_float Levinsondurbin(                               /* O    prediction error energy                                         */\r
305         SKP_float               A[],                            /* O    prediction coefficients [order]                         */\r
306         const SKP_float corr[],                         /* I    input auto-correlations [order + 1]                     */\r
307         const SKP_int   order                           /* I    prediction order                                                        */\r
308 )\r
309 {\r
310         SKP_int   i, mHalf, m;\r
311         SKP_float min_nrg, nrg, t, km, Atmp1, Atmp2;\r
312         \r
313         min_nrg = 1e-12f * corr[ 0 ] + 1e-9f;\r
314         nrg = corr[ 0 ];\r
315         nrg = SKP_max(min_nrg, nrg);\r
316         A[ 0 ] = corr[ 1 ] / nrg;\r
317         nrg -= A[ 0 ] * corr[ 1 ];\r
318         nrg = SKP_max(min_nrg, nrg);\r
319 \r
320         for( m = 1; m < order; m++ )\r
321         {\r
322                 t = corr[ m + 1 ];\r
323                 for( i = 0; i < m; i++ ) {\r
324                         t -= A[ i ] * corr[ m - i ];\r
325                 }\r
326 \r
327                 /* reflection coefficient */\r
328                 km = t / nrg;\r
329 \r
330                 /* residual energy */\r
331                 nrg -= km * t;\r
332                 nrg = SKP_max(min_nrg, nrg);\r
333 \r
334                 mHalf = m >> 1;\r
335                 for( i = 0; i < mHalf; i++ ) {\r
336                         Atmp1 = A[ i ];\r
337                         Atmp2 = A[ m - i - 1 ];\r
338                         A[ m - i - 1 ] -= km * Atmp1;\r
339                         A[ i ]         -= km * Atmp2;\r
340                 }\r
341                 if( m & 1 ) {\r
342                         A[ mHalf ]     -= km * A[ mHalf ];\r
343                 }\r
344                 A[ m ] = km;\r
345         }\r
346 \r
347         /* return the residual energy */\r
348         return nrg;\r
349 }\r
350 \r
351 /* Chirp (bw expand) LP AR filter */\r
352 void Bwexpander( \r
353         SKP_float *ar,                                          /* io   AR filter to be expanded (without leading 1)   */\r
354         const SKP_int d,                                        /* i    length of ar                                                                    */\r
355         const SKP_float chirp                           /* i    chirp factor (typically in range (0..1) )               */\r
356 )\r
357 {\r
358         SKP_int   i;\r
359         SKP_float cfac = chirp;\r
360 \r
361         for( i = 0; i < d - 1; i++ ) {\r
362                 ar[ i ] *=  cfac;\r
363                 cfac    *=  chirp;\r
364         }\r
365         ar[ d - 1 ] *=  cfac;\r
366 }\r