Change CWRS indexing to use Pyramid VQ's magnitude ordering.
[opus.git] / libcelt / cwrs.c
1 /* (C) 2007-2008 Timothy B. Terriberry
2    (C) 2008 Jean-Marc Valin */
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    - Redistributions of source code must retain the above copyright
9    notice, this list of conditions and the following disclaimer.
10
11    - Redistributions in binary form must reproduce the above copyright
12    notice, this list of conditions and the following disclaimer in the
13    documentation and/or other materials provided with the distribution.
14
15    - Neither the name of the Xiph.org Foundation nor the names of its
16    contributors may be used to endorse or promote products derived from
17    this software without specific prior written permission.
18
19    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
20    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
21    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
22    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
23    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
24    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
25    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
26    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
27    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
28    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
29    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
30 */
31
32 /* Functions for encoding and decoding pulse vectors.
33    These are based on the function
34      U(n,m) = U(n-1,m) + U(n,m-1) + U(n-1,m-1),
35      U(n,1) = U(1,m) = 2,
36     which counts the number of ways of placing m pulses in n dimensions, where
37      at least one pulse lies in dimension 0.
38    For more details, see: http://people.xiph.org/~tterribe/notes/cwrs.html
39 */
40
41 #ifdef HAVE_CONFIG_H
42 #include "config.h"
43 #endif
44
45 #include "os_support.h"
46 #include <stdlib.h>
47 #include <string.h>
48 #include "cwrs.h"
49 #include "mathops.h"
50
51 #if 0
52 int log2_frac(ec_uint32 val, int frac)
53 {
54    int i;
55    /* EC_ILOG() actually returns log2()+1, go figure */
56    int L = EC_ILOG(val)-1;
57    /*printf ("in: %d %d ", val, L);*/
58    if (L>14)
59       val >>= L-14;
60    else if (L<14)
61       val <<= 14-L;
62    L <<= frac;
63    /*printf ("%d\n", val);*/
64    for (i=0;i<frac;i++)
65 {
66       val = (val*val) >> 15;
67       /*printf ("%d\n", val);*/
68       if (val > 16384)
69          L |= (1<<(frac-i-1));
70       else   
71          val <<= 1;
72 }
73    return L;
74 }
75 #endif
76
77 int log2_frac64(ec_uint64 val, int frac)
78 {
79    int i;
80    /* EC_ILOG64() actually returns log2()+1, go figure */
81    int L = EC_ILOG64(val)-1;
82    /*printf ("in: %d %d ", val, L);*/
83    if (L>14)
84       val >>= L-14;
85    else if (L<14)
86       val <<= 14-L;
87    L <<= frac;
88    /*printf ("%d\n", val);*/
89    for (i=0;i<frac;i++)
90    {
91       val = (val*val) >> 15;
92       /*printf ("%d\n", val);*/
93       if (val > 16384)
94          L |= (1<<(frac-i-1));
95       else   
96          val <<= 1;
97    }
98    return L;
99 }
100
101 int fits_in32(int _n, int _m)
102 {
103    static const celt_int16_t maxN[15] = {
104       255, 255, 255, 255, 255, 109,  60,  40,
105        29,  24,  20,  18,  16,  14,  13};
106    static const celt_int16_t maxM[15] = {
107       255, 255, 255, 255, 255, 238,  95,  53,
108        36,  27,  22,  18,  16,  15,  13};
109    if (_n>=14)
110    {
111       if (_m>=14)
112          return 0;
113       else
114          return _n <= maxN[_m];
115    } else {
116       return _m <= maxM[_n];
117    }   
118 }
119 int fits_in64(int _n, int _m)
120 {
121    static const celt_int16_t maxN[28] = {
122       255, 255, 255, 255, 255, 255, 255, 255,
123       255, 255, 178, 129, 100,  81,  68,  58,
124        51,  46,  42,  38,  36,  33,  31,  30,
125        28, 27, 26, 25};
126    static const celt_int16_t maxM[28] = {
127       255, 255, 255, 255, 255, 255, 255, 255, 
128       255, 255, 245, 166, 122,  94,  77,  64, 
129        56,  49,  44,  40,  37,  34,  32,  30,
130        29,  27,  26,  25};
131    if (_n>=27)
132    {
133       if (_m>=27)
134          return 0;
135       else
136          return _n <= maxN[_m];
137    } else {
138       return _m <= maxM[_n];
139    }
140 }
141
142 #define MASK32 (0xFFFFFFFF)
143
144 /*INV_TABLE[i] holds the multiplicative inverse of (2*i-1) mod 2**32.*/
145 static const unsigned INV_TABLE[128]={
146   0x00000001,0xAAAAAAAB,0xCCCCCCCD,0xB6DB6DB7,
147   0x38E38E39,0xBA2E8BA3,0xC4EC4EC5,0xEEEEEEEF,
148   0xF0F0F0F1,0x286BCA1B,0x3CF3CF3D,0xE9BD37A7,
149   0xC28F5C29,0x684BDA13,0x4F72C235,0xBDEF7BDF,
150   0x3E0F83E1,0x8AF8AF8B,0x914C1BAD,0x96F96F97,
151   0xC18F9C19,0x2FA0BE83,0xA4FA4FA5,0x677D46CF,
152   0x1A1F58D1,0xFAFAFAFB,0x8C13521D,0x586FB587,
153   0xB823EE09,0xA08AD8F3,0xC10C9715,0xBEFBEFBF,
154   0xC0FC0FC1,0x07A44C6B,0xA33F128D,0xE327A977,
155   0xC7E3F1F9,0x962FC963,0x3F2B3885,0x613716AF,
156   0x781948B1,0x2B2E43DB,0xFCFCFCFD,0x6FD0EB67,
157   0xFA3F47E9,0xD2FD2FD3,0x3F4FD3F5,0xD4E25B9F,
158   0x5F02A3A1,0xBF5A814B,0x7C32B16D,0xD3431B57,
159   0xD8FD8FD9,0x8D28AC43,0xDA6C0965,0xDB195E8F,
160   0x0FDBC091,0x61F2A4BB,0xDCFDCFDD,0x46FDD947,
161   0x56BE69C9,0xEB2FDEB3,0x26E978D5,0xEFDFBF7F,
162   0x0FE03F81,0xC9484E2B,0xE133F84D,0xE1A8C537,
163   0x077975B9,0x70586723,0xCD29C245,0xFAA11E6F,
164   0x0FE3C071,0x08B51D9B,0x8CE2CABD,0xBF937F27,
165   0xA8FE53A9,0x592FE593,0x2C0685B5,0x2EB11B5F,
166   0xFCD1E361,0x451AB30B,0x72CFE72D,0xDB35A717,
167   0xFB74A399,0xE80BFA03,0x0D516325,0x1BCB564F,
168   0xE02E4851,0xD962AE7B,0x10F8ED9D,0x95AEDD07,
169   0xE9DC0589,0xA18A4473,0xEA53FA95,0xEE936F3F,
170   0x90948F41,0xEAFEAFEB,0x3D137E0D,0xEF46C0F7,
171   0x028C1979,0x791064E3,0xC04FEC05,0xE115062F,
172   0x32385831,0x6E68575B,0xA10D387D,0x6FECF2E7,
173   0x3FB47F69,0xED4BFB53,0x74FED775,0xDB43BB1F,
174   0x87654321,0x9BA144CB,0x478BBCED,0xBFB912D7,
175   0x1FDCD759,0x14B2A7C3,0xCB125CE5,0x437B2E0F,
176   0x10FEF011,0xD2B3183B,0x386CAB5D,0xEF6AC0C7,
177   0x0E64C149,0x9A020A33,0xE6B41C55,0xFEFEFEFF
178 };
179
180 /*Computes (_a*_b-_c)/(2*_d-1) when the quotient is known to be exact.
181   _a, _b, _c, and _d may be arbitrary so long as the arbitrary precision result
182    fits in 32 bits, but currently the table for multiplicative inverses is only
183    valid for _d<128.*/
184 static inline celt_uint32_t imusdiv32odd(celt_uint32_t _a,celt_uint32_t _b,
185  celt_uint32_t _c,celt_uint32_t _d){
186   return (_a*_b-_c)*INV_TABLE[_d]&MASK32;
187 }
188
189 /*Computes (_a*_b-_c)/_d when the quotient is known to be exact.
190   _d does not actually have to be even, but imusdiv32odd will be faster when
191    it's odd, so you should use that instead.
192   _a and _d are assumed to be small (e.g., _a*_d fits in 32 bits; currently the
193    table for multiplicative inverses is only valid for _d<256).
194   _b and _c may be arbitrary so long as the arbitrary precision reuslt fits in
195    32 bits.*/
196 static inline celt_uint32_t imusdiv32even(celt_uint32_t _a,celt_uint32_t _b,
197  celt_uint32_t _c,celt_uint32_t _d){
198   unsigned inv;
199   int      mask;
200   int      shift;
201   int      one;
202   shift=EC_ILOG(_d^_d-1);
203   inv=INV_TABLE[_d-1>>shift];
204   shift--;
205   one=1<<shift;
206   mask=one-1;
207   return (_a*(_b>>shift)-(_c>>shift)+
208    (_a*(_b&mask)+one-(_c&mask)>>shift)-1)*inv&MASK32;
209 }
210
211 /*Computes the next row/column of any recurrence that obeys the relation
212    u[i][j]=u[i-1][j]+u[i][j-1]+u[i-1][j-1].
213   _ui0 is the base case for the new row/column.*/
214 static inline void unext32(celt_uint32_t *_ui,int _len,celt_uint32_t _ui0){
215   celt_uint32_t ui1;
216   int           j;
217   /* doing a do-while would overrun the array if we had less than 2 samples */
218   j=1; do {
219     ui1=_ui[j]+_ui[j-1]+_ui0;
220     _ui[j-1]=_ui0;
221     _ui0=ui1;
222   } while (++j<_len);
223   _ui[j-1]=_ui0;
224 }
225
226 static inline void unext64(celt_uint64_t *_ui,int _len,celt_uint64_t _ui0){
227   celt_uint64_t ui1;
228   int           j;
229   /* doing a do-while would overrun the array if we had less than 2 samples */
230   j=1; do {
231     ui1=_ui[j]+_ui[j-1]+_ui0;
232     _ui[j-1]=_ui0;
233     _ui0=ui1;
234   } while (++j<_len);
235   _ui[j-1]=_ui0;
236 }
237
238 /*Computes the previous row/column of any recurrence that obeys the relation
239    u[i-1][j]=u[i][j]-u[i][j-1]-u[i-1][j-1].
240   _ui0 is the base case for the new row/column.*/
241 static inline void uprev32(celt_uint32_t *_ui,int _n,celt_uint32_t _ui0){
242   celt_uint32_t ui1;
243   int           j;
244   /* doing a do-while would overrun the array if we had less than 2 samples */
245   j=1; do {
246     ui1=_ui[j]-_ui[j-1]-_ui0;
247     _ui[j-1]=_ui0;
248     _ui0=ui1;
249   } while (++j<_n);
250   _ui[j-1]=_ui0;
251 }
252
253 static inline void uprev64(celt_uint64_t *_ui,int _n,celt_uint64_t _ui0){
254   celt_uint64_t ui1;
255   int           j;
256   /* doing a do-while would overrun the array if we had less than 2 samples */
257   j=1; do {
258     ui1=_ui[j]-_ui[j-1]-_ui0;
259     _ui[j-1]=_ui0;
260     _ui0=ui1;
261   } while (++j<_n);
262   _ui[j-1]=_ui0;
263 }
264
265 /*Returns the number of ways of choosing _m elements from a set of size _n with
266    replacement when a sign bit is needed for each unique element.
267   On input, _u should be initialized to column (_m-1) of U(n,m).
268   On exit, _u will be initialized to column _m of U(n,m).*/
269 celt_uint32_t ncwrs_unext32(int _n,celt_uint32_t *_ui){
270   celt_uint32_t ret;
271   celt_uint32_t ui0;
272   celt_uint32_t ui1;
273   int           j;
274   ret=ui0=2;
275   celt_assert(_n>=2);
276   j=1; do {
277     ui1=_ui[j]+_ui[j-1]+ui0;
278     _ui[j-1]=ui0;
279     ui0=ui1;
280     ret+=ui0;
281   } while (++j<_n);
282   _ui[j-1]=ui0;
283   return ret;
284 }
285
286 celt_uint64_t ncwrs_unext64(int _n,celt_uint64_t *_ui){
287   celt_uint64_t ret;
288   celt_uint64_t ui0;
289   celt_uint64_t ui1;
290   int           j;
291   ret=ui0=1;
292   celt_assert(_n>=2);
293   j=1; do {
294     ui1=_ui[j]+_ui[j-1]+ui0;
295     _ui[j-1]=ui0;
296     ui0=ui1;
297     ret+=ui0;
298   } while (++j<_n);
299   _ui[j-1]=ui0;
300   return ret<<=1;
301 }
302
303 /*Returns the number of ways of choosing _m elements from a set of size _n with
304    replacement when a sign bit is needed for each unique element.
305   _u: On exit, _u[i] contains U(i+1,_m).*/
306 celt_uint32_t ncwrs_u32(int _n,int _m,celt_uint32_t *_u){
307   celt_uint32_t ret;
308   celt_uint32_t um2;
309   int           k;
310   /*If _m==0, _u[] should be set to zero and the return should be 1.*/
311   celt_assert(_m>0);
312   /*We'll overflow our buffer unless _n>=2.*/
313   celt_assert(_n>=2);
314   um2=_u[0]=1;
315   if(_m<=6){
316     if(_m<2){
317       k=1;
318       do _u[k]=1;
319       while(++k<_n);
320     }
321     else{
322       k=1;
323       do _u[k]=(k<<1)+1;
324       while(++k<_n);
325       for(k=2;k<_m;k++)unext32(_u,_n,1);
326     }
327   }
328   else{
329     celt_uint32_t um1;
330     celt_uint32_t n2m1;
331     _u[1]=n2m1=um1=(_m<<1)-1;
332     for(k=2;k<_n;k++){
333       /*U(n,m) = ((2*n-1)*U(n,m-1)-U(n,m-2))/(m-1) + U(n,m-2)*/
334       _u[k]=um2=imusdiv32even(n2m1,um1,um2,k)+um2;
335       if(++k>=_n)break;
336       _u[k]=um1=imusdiv32odd(n2m1,um2,um1,k>>1)+um1;
337     }
338   }
339   ret=1;
340   k=1;
341   do ret+=_u[k];
342   while(++k<_n);
343   return ret<<1;
344 }
345
346 celt_uint64_t ncwrs_u64(int _n,int _m,celt_uint64_t *_u){
347   int k;
348   CELT_MEMSET(_u,0,_n);
349   if(_m<=0)return 1;
350   if(_n<=0)return 0;
351   for(k=1;k<_m;k++)unext64(_u,_n,1);
352   return ncwrs_unext64(_n,_u);
353 }
354
355
356 /*Returns the _i'th combination of _m elements chosen from a set of size _n
357    with associated sign bits.
358   _y: Returns the vector of pulses.
359   _u: Must contain entries [1..._n] of column _m of U() on input.
360       Its contents will be destructively modified.*/
361 void cwrsi32(int _n,int _m,celt_uint32_t _i,celt_uint32_t _nc,int *_y,
362  celt_uint32_t *_u){
363   celt_uint32_t p;
364   celt_uint32_t q;
365   int           j;
366   int           k;
367   celt_assert(_n>0);
368   p=_nc;
369   q=0;
370   j=0;
371   k=_m;
372   do{
373     int s;
374     int yj;
375     p-=q;
376     q=_u[_n-j-1];
377     p-=q;
378     s=_i>=p;
379     if(s)_i-=p;
380     yj=k;
381     while(q>_i){
382       uprev32(_u,_n-j,--k>0);
383       p=q;
384       q=_u[_n-j-1];
385     }
386     _i-=q;
387     yj-=k;
388     _y[j]=yj-(yj<<1&-s);
389   }
390   while(++j<_n);
391 }
392
393 /*Returns the _i'th combination of _m elements chosen from a set of size _n
394    with associated sign bits.
395   _y: Returns the vector of pulses.
396   _u: Must contain entries [1..._n] of column _m of U() on input.
397       Its contents will be destructively modified.*/
398 void cwrsi64(int _n,int _m,celt_uint64_t _i,celt_uint64_t _nc,int *_y,
399  celt_uint64_t *_u){
400   celt_uint64_t p;
401   celt_uint64_t q;
402   int           j;
403   int           k;
404   celt_assert(_n>0);
405   p=_nc;
406   q=0;
407   j=0;
408   k=_m;
409   do{
410     int s;
411     int yj;
412     p-=q;
413     q=_u[_n-j-1];
414     p-=q;
415     s=_i>=p;
416     if(s)_i-=p;
417     yj=k;
418     while(q>_i){
419       uprev64(_u,_n-j,--k>0);
420       p=q;
421       q=_u[_n-j-1];
422     }
423     _i-=q;
424     yj-=k;
425     _y[j]=yj-(yj<<1&-s);
426   }
427   while(++j<_n);
428 }
429
430 /*Returns the index of the given combination of _m elements chosen from a set
431    of size _n with associated sign bits.
432   _y:  The vector of pulses, whose sum of absolute values must be _m.
433   _nc: Returns V(_n,_m).*/
434 celt_uint32_t icwrs32(int _n,int _m,celt_uint32_t *_nc,const int *_y,
435  celt_uint32_t *_u){
436   celt_uint32_t nc;
437   celt_uint32_t i;
438   int           j;
439   int           k;
440   /*We can't unroll the first two iterations of the loop unless _n>=2.*/
441   celt_assert(_n>=2);
442   nc=1;
443   i=_y[_n-1]<0;
444   _u[0]=0;
445   for(k=1;k<=_m+1;k++)_u[k]=(k<<1)-1;
446   k=abs(_y[_n-1]);
447   j=_n-2;
448   nc+=_u[_m];
449   i+=_u[k];
450   k+=abs(_y[j]);
451   if(_y[j]<0)i+=_u[k+1];
452   while(j-->0){
453     unext32(_u,_m+2,0);
454     nc+=_u[_m];
455     i+=_u[k];
456     k+=abs(_y[j]);
457     if(_y[j]<0)i+=_u[k+1];
458   }
459   /*If _m==0, nc should not be doubled.*/
460   celt_assert(_m>0);
461   *_nc=nc<<1;
462   return i;
463 }
464
465 /*Returns the index of the given combination of _m elements chosen from a set
466    of size _n with associated sign bits.
467   _y:  The vector of pulses, whose sum of absolute values must be _m.
468   _nc: Returns V(_n,_m).*/
469 celt_uint64_t icwrs64(int _n,int _m,celt_uint64_t *_nc,const int *_y,
470  celt_uint64_t *_u){
471   celt_uint64_t nc;
472   celt_uint64_t i;
473   int           j;
474   int           k;
475   /*We can't unroll the first two iterations of the loop unless _n>=2.*/
476   celt_assert(_n>=2);
477   nc=1;
478   i=_y[_n-1]<0;
479   _u[0]=0;
480   for(k=1;k<=_m+1;k++)_u[k]=(k<<1)-1;
481   k=abs(_y[_n-1]);
482   j=_n-2;
483   nc+=_u[_m];
484   i+=_u[k];
485   k+=abs(_y[j]);
486   if(_y[j]<0)i+=_u[k+1];
487   while(j-->0){
488     unext64(_u,_m+2,0);
489     nc+=_u[_m];
490     i+=_u[k];
491     k+=abs(_y[j]);
492     if(_y[j]<0)i+=_u[k+1];
493   }
494   /*If _m==0, nc should not be doubled.*/
495   celt_assert(_m>0);
496   *_nc=nc<<1;
497   return i;
498 }
499
500 static inline void encode_pulse32(int _n,int _m,const int *_y,ec_enc *_enc){
501   VARDECL(celt_uint32_t,u);
502   celt_uint32_t nc;
503   celt_uint32_t i;
504   SAVE_STACK;
505   ALLOC(u,_m+2,celt_uint32_t);
506   i=icwrs32(_n,_m,&nc,_y,u);
507   ec_enc_uint(_enc,i,nc);
508   RESTORE_STACK;
509 }
510
511 static inline void encode_pulse64(int _n,int _m,const int *_y,ec_enc *_enc){
512   VARDECL(celt_uint64_t,u);
513   celt_uint64_t nc;
514   celt_uint64_t i;
515   SAVE_STACK;
516   ALLOC(u,_m+2,celt_uint64_t);
517   i=icwrs64(_n,_m,&nc,_y,u);
518   ec_enc_uint64(_enc,i,nc);
519   RESTORE_STACK;
520 }
521
522 int get_required_bits(int N, int K, int frac)
523 {
524    int nbits = 0;
525    if(fits_in64(N,K))
526    {
527       VARDECL(celt_uint64_t,u);
528       SAVE_STACK;
529       ALLOC(u,N,celt_uint64_t);
530       nbits = log2_frac64(ncwrs_u64(N,K,u), frac);
531       RESTORE_STACK;
532    } else {
533       nbits = log2_frac64(N, frac);
534       nbits += get_required_bits(N/2+1, (K+1)/2, frac);
535       nbits += get_required_bits(N/2+1, K/2, frac);
536    }
537    return nbits;
538 }
539
540
541 void encode_pulses(int *_y, int N, int K, ec_enc *enc)
542 {
543    if (K==0) {
544    } else if (N==1)
545    {
546       ec_enc_bits(enc, _y[0]<0, 1);
547    } else if(fits_in32(N,K))
548    {
549       encode_pulse32(N, K, _y, enc);
550    } else if(fits_in64(N,K)) {
551       encode_pulse64(N, K, _y, enc);
552    } else {
553      int i;
554      int count=0;
555      int split;
556      split = (N+1)/2;
557      for (i=0;i<split;i++)
558         count += abs(_y[i]);
559      ec_enc_uint(enc,count,K+1);
560      encode_pulses(_y, split, count, enc);
561      encode_pulses(_y+split, N-split, K-count, enc);
562    }
563 }
564
565 static inline void decode_pulse32(int _n,int _m,int *_y,ec_dec *_dec){
566   VARDECL(celt_uint32_t,u);
567   celt_uint32_t nc;
568   SAVE_STACK;
569   ALLOC(u,_n,celt_uint32_t);
570   nc=ncwrs_u32(_n,_m,u);
571   cwrsi32(_n,_m,ec_dec_uint(_dec,nc),nc,_y,u);
572   RESTORE_STACK;
573 }
574
575 static inline void decode_pulse64(int _n,int _m,int *_y,ec_dec *_dec){
576   VARDECL(celt_uint64_t,u);
577   celt_uint64_t nc;
578   SAVE_STACK;
579   ALLOC(u,_n,celt_uint64_t);
580   nc=ncwrs_u64(_n,_m,u);
581   cwrsi64(_n,_m,ec_dec_uint64(_dec,nc),nc,_y,u);
582   RESTORE_STACK;
583 }
584
585 void decode_pulses(int *_y, int N, int K, ec_dec *dec)
586 {
587    if (K==0) {
588       int i;
589       for (i=0;i<N;i++)
590          _y[i] = 0;
591    } else if (N==1)
592    {
593       int s = ec_dec_bits(dec, 1);
594       if (s==0)
595          _y[0] = K;
596       else
597          _y[0] = -K;
598    } else if(fits_in32(N,K))
599    {
600       decode_pulse32(N, K, _y, dec);
601    } else if(fits_in64(N,K)) {
602       decode_pulse64(N, K, _y, dec);
603    } else {
604      int split;
605      int count = ec_dec_uint(dec,K+1);
606      split = (N+1)/2;
607      decode_pulses(_y, split, count, dec);
608      decode_pulses(_y+split, N-split, K-count, dec);
609    }
610 }