1 /*
2 * Musepack audio compression
3 * Copyright (c) 2005-2009, The Musepack Development Team
4 * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf
5 *
6 * This library is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU Lesser General Public
8 * License as published by the Free Software Foundation; either
9 * version 2.1 of the License, or (at your option) any later version.
10 *
11 * This library is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with this library; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 */
20
21 #include "mpcenc.h"
22
23
24 #define MAX_LPC_ORDER 35
25 #define log2(x) ( log (x) * (1./M_LN2) )
26 #define ORDER_PENALTY 0
27
28
29 static int // best prediction order model
CalculateLPCCoeffs(Int32_t * buf,size_t nbuf,Int32_t offset,double * lpcout,int nlpc,float * psigbit,float * presbit)30 CalculateLPCCoeffs ( Int32_t* buf, // Samples
31 size_t nbuf, // Number of samples
32 Int32_t offset, //
33 double* lpcout, // quantized prediction coefficients
34 int nlpc, // max. prediction order
35 float* psigbit, // expected number of bits per original signal sample
36 float* presbit ) // expected number of bits per residual signal sample
37 {
38 static double* fbuf = NULL;
39 static int nflpc = 0;
40 static int nfbuf = 0;
41 int nbit;
42 int i;
43 int j;
44 int bestnbit;
45 int bestnlpc;
46 double e;
47 double bestesize;
48 double ci;
49 double esize;
50 double acf [MAX_LPC_ORDER + 1];
51 double ref [MAX_LPC_ORDER + 1];
52 double lpc [MAX_LPC_ORDER + 1];
53 double tmp [MAX_LPC_ORDER + 1];
54 double escale = 0.5 * M_LN2 * M_LN2 / nbuf;
55 double sum;
56
57 if ( nlpc >= nbuf ) // if necessary, limit the LPC order to the number of samples available
58 nlpc = nbuf - 1;
59
60 if ( nlpc > nflpc || nbuf > nfbuf ) { // grab some space for a 'zero mean' buffer of floats if needed
61 if ( fbuf != NULL )
62 free ( fbuf - nflpc );
63 fbuf = nlpc + ((double*) calloc ( nlpc+nbuf, sizeof (*fbuf) ));
64 nfbuf = nbuf;
65 nflpc = nlpc;
66 }
67
68 e = 0.;
69 for ( j = 0; j < nbuf; j++ ) { // zero mean signal and compute energy
70 sum = fbuf [j] = buf[j] - (double)offset;
71 e += sum * sum;
72 }
73
74 esize = e > 0. ? 0.5 * log2 (escale * e) : 0.;
75 *psigbit = esize; // return the expected number of bits per original signal sample
76
77 acf [0] = e; // store the best values so far (the zeroth order predictor)
78 bestnlpc = 0;
79 bestnbit = nbuf * esize;
80 bestesize = esize;
81
82 for ( i = 1; i <= nlpc && e > 0. && i < bestnlpc + 4; i++ ) { // just check two more than bestnlpc
83
84 sum = 0.;
85 for ( j = i; j < nbuf; j++ ) // compute the jth autocorrelation coefficient
86 sum += fbuf [j] * fbuf [j-i];
87 acf [i] = sum;
88
89 ci = 0.; // compute the reflection and LP coeffients for order j predictor
90 for ( j = 1; j < i; j++ )
91 ci += lpc [j] * acf [i-j];
92 lpc [i] = ref [i] = ci = (acf [i] - ci) / e;
93 for ( j = 1; j < i; j++ )
94 tmp [j] = lpc [j] - ci * lpc [i-j];
95 for ( j = 1; j < i; j++ )
96 lpc [j] = tmp [j];
97
98 e *= 1 - ci*ci; // compute the new energy in the prediction residual
99 esize = e > 0. ? 0.5 * log2 (escale * e) : 0.;
100
101 nbit = nbuf * esize + i * ORDER_PENALTY;
102 if ( nbit < bestnbit ) { // store this model if it is the best so far
103 bestnlpc = i; // store best model order
104 bestnbit = nbit;
105 bestesize = esize;
106
107 for ( j = 0; j < bestnlpc; j++ ) // store the quantized LP coefficients
108 lpcout [j] = lpc [j+1];
109 }
110 }
111
112 *presbit = bestesize; // return the expected number of bits per residual signal sample
113 return bestnlpc; // return the best model order
114 }
115
116
117 static void
Pred(const unsigned int * new,unsigned int * old)118 Pred ( const unsigned int* new,
119 unsigned int* old )
120 {
121 static Double DOUBLE [36];
122 Float org;
123 Float pred;
124 int i;
125 int j;
126 int sum = 18;
127 int order;
128 double oldeff = 0.;
129 double neweff = 0.;
130
131 for ( i = 0; i < 36; i++ )
132 sum += old [i];
133 sum = (int) floor (sum / 36.);
134
135 order = CalculateLPCCoeffs ( old, 36, sum*0, DOUBLE, 35, &org, &pred );
136
137 printf ("avg: %4u [%2u] %.2f %.2f\n\n", sum, order, org, pred );
138 if ( order < 1 )
139 return;
140
141 for ( i = 0; i < order; i++ )
142 printf ("%f ", DOUBLE[i] );
143 printf ("\n");
144
145 for ( i = 0; i < 36; i++ ) {
146 double sum = 0.;
147 for ( j = 1; j <= order; j++ ) {
148 sum += (i-j < 0 ? old[i-j+36] : new [i-j]) * DOUBLE [j-1];
149 }
150 printf ("%2u: %6.2f %3d\n", i, sum, new [i] );
151 oldeff += new[i] * new[i];
152 neweff += (sum-new[i]) * (sum-new[i]);
153 }
154 printf ("%6.2f %6.2f\n", sqrt(oldeff), sqrt(neweff) );
155 }
156
157
158 void
Predicate(int Channel,int Band,unsigned int * x,int * scf)159 Predicate ( int Channel, int Band, unsigned int* x, int* scf )
160 {
161 static Int32_t OLD [2] [32] [36];
162 int i;
163
164 printf ("=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=\n");
165 for ( i = 0; i < 36; i++ )
166 printf ("%2d ", OLD [Channel][Band][i] );
167 printf ("\n");
168 for ( i = 0; i < 36; i++ )
169 printf ("%2d ", x[i] );
170 printf ("\n");
171 printf ("%2u-%2u-%2u ", scf[0], scf[1], scf[2] );
172 Pred ( x, OLD [Channel][Band] );
173 for ( i = 0; i < 36; i++ )
174 OLD [Channel][Band][i] = x[i];
175 }
176
177 /* end of predict.c */
178