1 /* ************************************************************************
2 FILE: rounding.c
3 COMMENT: Functions related to filter coefficients rounding to
4 integer values or CSD form.
5 CONTENT:
6 AUTHOR: Lime Microsystems
7 DATE: Jun 18, 2000
8 REVISION: Nov 30, 2007: Tiger project
9 ************************************************************************ */
10 #include <math.h>
11 #include <stdio.h>
12
13 #include "rounding.h"
14
15 void int2csd(int a, int cprec, int csdprec, int * bincode, int * csdcode, int * csdcoder);
16 int csd2int(int cprec, int *code);
17
18 /* ************************************************************************
19 ************************************************************************ */
round2int(a,b,n,cprec)20 void round2int(a, b, n, cprec)
21 double *a, *b;
22 int n, cprec;
23 {
24 int i,k;
25
26 for(i=0; i<n; i++) {
27 k = a[i] > 0.0 ? 1 : -1;
28 b[i] = (int)(a[i]*(1<<cprec) + k*0.5);
29 b[i] /= (double)(1<<cprec);
30 }
31 }
32
33 /* ************************************************************************
34 ************************************************************************ */
round2csd(a,b,n,cprec,csdprec,bincode,csdcode,csdcoder)35 void round2csd(a, b, n, cprec, csdprec, bincode, csdcode, csdcoder)
36 double *a, *b;
37 int n, cprec, csdprec;
38 int **bincode, **csdcode, **csdcoder;
39 {
40 int i,k, ia;
41
42 for(i=0; i<n; i++) {
43 k = a[i] > 0.0 ? 1 : -1;
44 ia = (int)(a[i]*(1<<cprec) + k*0.5);
45 int2csd(ia, cprec, csdprec, bincode[i], csdcode[i], csdcoder[i]);
46 ia = csd2int(cprec, csdcoder[i]);
47 b[i] = (double)(ia)/(double)(1<<cprec);
48 }
49 }
50
51 /* ************************************************************************
52 ************************************************************************ */
printcode(int ** code,int n,int cprec)53 void printcode(int ** code, int n, int cprec)
54 {
55 int i, j;
56 double sumh, sume, sumo;
57 int ih; double h;
58 int negative, sign, shift, csdprec;
59 int symmetry;
60
61 /* Find maximum nonzero bits per coefficient */
62 csdprec = 0;
63 for(i=0; i<n; i++) {
64 negative = 0;
65 for(j=0; j<= cprec; j++) if(code[i][j] != 0) negative++;
66 if(negative > csdprec) csdprec = negative;
67 }
68
69 /* Check symmetry of the filter */
70 if( csd2int(cprec, code[0]) == csd2int(cprec,code[n-1]) ) symmetry = 1;
71 else symmetry = -1;
72
73 sumh=0.0; sume=0.0; sumo=0.0;
74 for(i=0; i<n; i++) {
75 ih = csd2int(cprec, code[i]);
76 h = (double)ih/(double)(1<<cprec);
77 sumh += fabs(h);
78 if(i%2) sumo += fabs(h);
79 else sume += fabs(h);
80
81 if ((ih != 0) && (i<(n+1)/2) ) {
82 negative = 0;
83 for(j=0; j<=cprec; j++) {
84 if(code[i][j] == -1) negative ++;
85 if(code[i][j] != 0 ) shift = j;
86 }
87
88 if(negative <= (csdprec-1)) sign = 1;
89 else sign = -1;
90
91 shift = cprec - shift;
92 if(fabs(h)*(1<<shift) > 1.0) shift --;
93
94 printf("h(%2d) = %11lg = %2d x (", i, h, sign);
95 for(j=cprec; j>=0; j--) {
96 if(sign*code[i][j] == 1) {
97 // printf(" +1/2^%d", cprec-j-shift);
98 printf(" +1/2^%d", cprec-j);
99 } else if(sign*code[i][j] == -1) {
100 // printf(" -1/2^%d", cprec-j-shift);
101 printf(" -1/2^%d", cprec-j);
102 }
103 }
104 // printf(" ) / ( 2^%d )\n", shift);
105 printf(" )\n");
106 }else if (ih != 0) {
107 printf("h(%2d) = %11lg = %2d x h(%2d)\n",
108 i, h, symmetry, n-1-i);
109 } else {
110 printf("h(%2d) = %11lg\n", i, 0.0);
111 }
112 }
113 printf("Sum of all abs(coefficients): %lg\n", sumh);
114 printf("Sum of even coefficients: %lg\n", sume);
115 printf("Sum of odd coefficients: %lg\n\n", sumo);
116 }
117
118 /* ************************************************************************
119 Print CSD code in the form of two common sub-expressions sharing
120 ************************************************************************ */
print_cses_code(xpx,xmx,x,n,cprec)121 void print_cses_code(xpx, xmx, x, n, cprec)
122 int **xpx, **xmx, **x, n, cprec;
123 {
124 int i, j;
125 int symmetry;
126 int ixpx, ixmx, ix;
127 double h;
128
129 /* Check symmetry of the filter */
130 if( (csd2int(cprec, xpx[0]) == csd2int(cprec,xpx[n-1])) &&
131 (csd2int(cprec, xmx[0]) == csd2int(cprec,xmx[n-1])) &&
132 (csd2int(cprec, x[0]) == csd2int(cprec,x[n-1])) ) symmetry = 1;
133 else symmetry = -1;
134
135 for(i=0; i<n; i++) {
136 ixpx = csd2int(cprec, xpx[i]);
137 ixmx = csd2int(cprec, xmx[i]);
138 ix = csd2int(cprec, x[i]);
139 h = (1.0+1.0/4.0)*(double)ixpx/(double)(1<<cprec) +
140 (1.0-1.0/4.0)*(double)ixmx/(double)(1<<cprec) +
141 (double)ix/(double)(1<<cprec);
142
143 if ((fpclassify(h) != FP_ZERO) && (i<(n+1)/2) ) {
144 printf("h(%2d) = %11lg = ", i, h);
145 if( ixpx ) {
146 printf("(1+1/4)x(");
147 for(j=cprec; j>=0; j--) {
148 if(xpx[i][j] == 1) {
149 printf(" +1/2^%d", cprec-j);
150 } else if(xpx[i][j] == -1) {
151 printf(" -1/2^%d", cprec-j);
152 }
153 }
154 printf(" )");
155 }
156
157 if( ixmx ) {
158 if( ixpx ) printf(" + (1-1/4)x(");
159 else printf("(1-1/4)x(");
160
161 for(j=cprec; j>=0; j--) {
162 if(xmx[i][j] == 1) {
163 printf(" +1/2^%d", cprec-j);
164 } else if(xmx[i][j] == -1) {
165 printf(" -1/2^%d", cprec-j);
166 }
167 }
168 printf(" )");
169 }
170
171 if( ix ) {
172 if( ixpx || ixmx ) printf(" + (");
173 else printf("(");
174
175 for(j=cprec; j>=0; j--) {
176 if(x[i][j] == 1) {
177 printf(" +1/2^%d", cprec-j);
178 } else if(x[i][j] == -1) {
179 printf(" -1/2^%d", cprec-j);
180 }
181 }
182 printf(" )");
183 }
184
185 printf("\n");
186 } else if (fpclassify(h) != FP_ZERO) {
187 printf("h(%2d) = %11lg = %2d x h(%2d)\n",
188 i, h, symmetry, n-1-i);
189 } else {
190 printf("h(%2d) = %11lg\n", i, 0.0);
191 }
192 }
193 }
194
195 /* ************************************************************************
196 ************************************************************************ */
int2csd(a,cprec,csdprec,bincode,csdcode,csdcoder)197 void int2csd(a, cprec, csdprec, bincode, csdcode, csdcoder)
198 int a; /* Input integer to be converted into CSD code */
199 int cprec; /* Integer precision */
200 int csdprec; /* CSD precistion */
201 int *bincode; /* Binary code */
202 int *csdcode; /* CSD code */
203 int *csdcoder; /* CSD code rounded to 'csdprec' nonzero bits */
204 {
205 int i, sign, ci, ci1, nzeroes;
206
207 if( a < 0 ) {
208 a *= -1;
209 sign = -1;
210 } else {
211 sign = 1;
212 }
213
214 /* Generate binary code of input */
215 for(i=0; i<cprec; i++) {
216 if( a & (1<<i) ) bincode[i] = 1;
217 else bincode[i] = 0;
218 }
219 bincode[cprec] = 0;
220
221 /* Construct CSD code */
222 ci = 0;
223 for(i=0; i<cprec; i++) {
224 if( (ci + bincode[i] + bincode[i+1]) > 1 ) ci1 = 1;
225 else ci1 = 0;
226
227 csdcode[i] = sign*(bincode[i] + ci -2*ci1);
228 bincode[i] *= sign;
229 ci = ci1;
230 }
231 csdcode[cprec] = sign*ci;
232
233 /* Round CSD code */
234 nzeroes = 0;
235 for(i=cprec; i >= 0 ; i--) {
236 if(csdcode[i] != 0) nzeroes ++;
237
238 if(nzeroes <= csdprec) csdcoder[i] = csdcode[i];
239 else csdcoder[i] = 0;
240 }
241
242 }
243
244 /* ************************************************************************
245 ************************************************************************ */
csd2int(cprec,code)246 int csd2int(cprec, code)
247 int cprec, *code;
248 {
249 int i, a;
250
251 a = 0;
252 for(i=cprec; i>=0; i--) a = a*2 + code[i];
253
254 return(a);
255 }
256
257 /* ************************************************************************
258 Extract x+x>>2 and x-x>>2 subexpressions from the CSD code
259 ************************************************************************ */
csesh(code,n,cprec,xpx,xmx,x)260 void csesh(code, n, cprec, xpx, xmx, x )
261 int n, cprec;
262 int **code, **xpx, **xmx, **x;
263 {
264 int i, k;
265
266 /* Set code matrices to zero */
267 for(i = 0; i < n; i++) {
268 for(k=0; k <= cprec; k++) {
269 xpx[i][k] = 0;
270 xmx[i][k] = 0;
271 x[i][k] = 0;
272 }
273 }
274
275 /* Extract two common subexpressions from all filter coefficients */
276 for(i = 0; i <n; i++) {
277 k = cprec;
278 while(1) {
279 /* Find next nonzero element in CSD code */
280 for(; k >= 0; k--) if(code[i][k] != 0) break;
281
282 if(k == -1) { /* There are no more nonzero digits */
283 break;
284 } else if(k == 0) { /* It is the last digit */
285 x[i][0] = code[i][0];
286 break;
287 } else if(k == 1) { /* Two more digits left */
288 x[i][0] = code[i][0];
289 x[i][1] = code[i][1];
290 break;
291 } else {
292 if((code[i][k] == 1) && (code[i][k-2] == 1) ) {
293 xpx[i][k] = 1;
294 code[i][k] = 0;
295 code[i][k-2] = 0;
296 } else if((code[i][k] == -1) && (code[i][k-2] == -1) ) {
297 xpx[i][k] = -1;
298 code[i][k] = 0;
299 code[i][k-2] = 0;
300 } else if((code[i][k] == 1) && (code[i][k-2] == -1) ) {
301 xmx[i][k] = 1;
302 code[i][k] = 0;
303 code[i][k-2] = 0;
304 } else if((code[i][k] == -1) && (code[i][k-2] == 1) ) {
305 xmx[i][k] = -1;
306 code[i][k] = 0;
307 code[i][k-2] = 0;
308 } else {
309 x[i][k] = code[i][k];
310 code[i][k] = 0;
311 }
312 }
313 }
314 }
315 }
316