1 /* K=9 r=1/3 Viterbi decoder for x86 MMX
2  * Aug 2006, Phil Karn, KA9Q
3  * May be used under the terms of the GNU Lesser General Public License (LGPL)
4  */
5 #include <mmintrin.h>
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <memory.h>
9 #include "fec.h"
10 
11 typedef union { unsigned char c[256]; __m64 v[32];} decision_t;
12 typedef union { unsigned short s[256]; __m64 v[64];} metric_t;
13 
14 static union branchtab39 { unsigned short s[128]; __m64 v[32];} Branchtab39[3];
15 static int Init = 0;
16 
17 /* State info for instance of Viterbi decoder */
18 struct v39 {
19   metric_t metrics1; /* path metric buffer 1 */
20   metric_t metrics2; /* path metric buffer 2 */
21   void *dp;          /* Pointer to current decision */
22   metric_t *old_metrics,*new_metrics; /* Pointers to path metrics, swapped on every bit */
23   void *decisions;   /* Beginning of decisions for block */
24 };
25 
26 /* Initialize Viterbi decoder for start of new frame */
init_viterbi39_mmx(void * p,int starting_state)27 int init_viterbi39_mmx(void *p,int starting_state){
28   struct v39 *vp = p;
29   int i;
30 
31   if(p == NULL)
32     return -1;
33   for(i=0;i<256;i++)
34     vp->metrics1.s[i] = 1000;
35 
36   vp->old_metrics = &vp->metrics1;
37   vp->new_metrics = &vp->metrics2;
38   vp->dp = vp->decisions;
39   vp->old_metrics->s[starting_state & 255] = 0; /* Bias known start state */
40   return 0;
41 }
42 
set_viterbi39_polynomial_mmx(int polys[3])43 void set_viterbi39_polynomial_mmx(int polys[3]){
44   int state;
45 
46   for(state=0;state < 128;state++){
47     Branchtab39[0].s[state] = (polys[0] < 0) ^ parity((2*state) & polys[0]) ? 255:0;
48     Branchtab39[1].s[state] = (polys[1] < 0) ^ parity((2*state) & polys[1]) ? 255:0;
49     Branchtab39[2].s[state] = (polys[2] < 0) ^ parity((2*state) & polys[2]) ? 255:0;
50   }
51   Init++;
52 }
53 
54 /* Create a new instance of a Viterbi decoder */
create_viterbi39_mmx(int len)55 void *create_viterbi39_mmx(int len){
56   struct v39 *vp;
57 
58   if(!Init){
59     int polys[3] = { V39POLYA,V39POLYB,V39POLYC };
60     set_viterbi39_polynomial_mmx(polys);
61   }
62   if((vp = (struct v39 *)malloc(sizeof(struct v39))) == NULL)
63     return NULL;
64   if((vp->decisions = malloc((len+8)*sizeof(decision_t))) == NULL){
65     free(vp);
66     return NULL;
67   }
68   init_viterbi39_mmx(vp,0);
69   return vp;
70 }
71 
72 
73 
74 /* Viterbi chainback */
chainback_viterbi39_mmx(void * p,unsigned char * data,unsigned int nbits,unsigned int endstate)75 int chainback_viterbi39_mmx(
76       void *p,
77       unsigned char *data, /* Decoded output data */
78       unsigned int nbits, /* Number of data bits */
79       unsigned int endstate){ /* Terminal encoder state */
80   struct v39 *vp = p;
81   decision_t *d;
82   int path_metric;
83 
84   if(p == NULL)
85     return -1;
86 
87   d = (decision_t *)vp->decisions;
88 
89   endstate %= 256;
90 
91   path_metric = vp->old_metrics->s[endstate];
92 
93   /* The store into data[] only needs to be done every 8 bits.
94    * But this avoids a conditional branch, and the writes will
95    * combine in the cache anyway
96    */
97   d += 8; /* Look past tail */
98   while(nbits-- != 0){
99     int k;
100 
101     k = d[nbits].c[endstate] & 1;
102     endstate = (k << 7) | (endstate >> 1);
103     data[nbits>>3] = endstate;
104   }
105   return path_metric;
106 }
107 
108 /* Delete instance of a Viterbi decoder */
delete_viterbi39_mmx(void * p)109 void delete_viterbi39_mmx(void *p){
110   struct v39 *vp = p;
111 
112   if(vp != NULL){
113     free(vp->decisions);
114     free(vp);
115   }
116 }
117 
118 
update_viterbi39_blk_mmx(void * p,unsigned char * syms,int nbits)119 int update_viterbi39_blk_mmx(void *p,unsigned char *syms,int nbits){
120   struct v39 *vp = p;
121   decision_t *d;
122   int path_metric = 0;
123 
124   if(p == NULL)
125     return -1;
126 
127   d = (decision_t *)vp->dp;
128 
129   while(nbits--){
130     __m64 sym0v,sym1v,sym2v;
131     void *tmp;
132     int i;
133 
134     /* Splat the 0th symbol across sym0v, the 1st symbol across sym1v, etc */
135     sym0v = _mm_set1_pi16(syms[0]);
136     sym1v = _mm_set1_pi16(syms[1]);
137     sym2v = _mm_set1_pi16(syms[2]);
138     syms += 3;
139 
140     for(i=0;i<32;i++){
141       __m64 decision0,decision1,metric,m_metric,m0,m1,m2,m3,survivor0,survivor1;
142 
143       /* Form branch metrics
144        * Because Branchtab takes on values 0 and 255, and the values of sym?v are offset binary in the range 0-255,
145        * the XOR operations constitute conditional negation.
146        * metric and m_metric (-metric) are in the range 0-1530
147        */
148       m0 = _mm_add_pi16(_mm_xor_si64(Branchtab39[0].v[i],sym0v),_mm_xor_si64(Branchtab39[1].v[i],sym1v));
149       metric = _mm_add_pi16(_mm_xor_si64(Branchtab39[2].v[i],sym2v),m0);
150       m_metric = _mm_sub_pi16(_mm_set1_pi16(765),metric);
151 
152       /* Add branch metrics to path metrics */
153       m0 = _mm_add_pi16(vp->old_metrics->v[i],metric);
154       m3 = _mm_add_pi16(vp->old_metrics->v[32+i],metric);
155       m1 = _mm_add_pi16(vp->old_metrics->v[32+i],m_metric);
156       m2 = _mm_add_pi16(vp->old_metrics->v[i],m_metric);
157 
158       /* Compare and select
159        * There's no packed min instruction in MMX, so we use modulo arithmetic
160        * to form the decisions and then do the select the hard way
161        */
162       decision0 = _mm_cmpgt_pi16(_mm_sub_pi16(m0,m1),_mm_setzero_si64());
163       decision1 = _mm_cmpgt_pi16(_mm_sub_pi16(m2,m3),_mm_setzero_si64());
164       survivor0 = _mm_or_si64(_mm_and_si64(decision0,m1),_mm_andnot_si64(decision0,m0));
165       survivor1 = _mm_or_si64(_mm_and_si64(decision1,m3),_mm_andnot_si64(decision1,m2));
166 
167       /* Merge decisions and store as bytes */
168       d->v[i] = _mm_unpacklo_pi8(_mm_packs_pi16(decision0,_mm_setzero_si64()),_mm_packs_pi16(decision1,_mm_setzero_si64()));
169 
170       /* Store surviving metrics */
171       vp->new_metrics->v[2*i] = _mm_unpacklo_pi16(survivor0,survivor1);
172       vp->new_metrics->v[2*i+1] = _mm_unpackhi_pi16(survivor0,survivor1);
173     }
174     if(vp->new_metrics->s[0] < vp->old_metrics->s[0])
175       path_metric += 65536; /* Hack: wraparound probably occured */
176     d++;
177     /* Swap pointers to old and new metrics */
178     tmp = vp->old_metrics;
179     vp->old_metrics = vp->new_metrics;
180     vp->new_metrics = tmp;
181   }
182   vp->dp = d;
183   _mm_empty();
184   return path_metric;
185 }
186