##// END OF EJS Templates
HG_REPOSITORIES » LPP » INSTRUMENTATION » SOLO_LFR » LFR_basic-parameters
RSS

Clone from http://pc-instru.lpp.polytechnique.fr:5000/LFR_basic-parameters

Sync
paul -
r4:294fc10efc0e default
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,840 +1,858
1 1 // In the frame of RPW LFR Sofware ICD Issue1 Rev8 (05/07/2013)
2 2 // version 1: 31/07/2013
3 3
4 4 #include "basic_parameters.h"
5 5 #include <math.h>
6 6 #include <stdio.h>
7 7
8 8 #define K44_PE 0
9 9 #define K55_PE 1
10 10 #define K45_PE_RE 2
11 11 #define K45_PE_IM 3
12 12
13 13 #define K14_SX_RE 4
14 14 #define K14_SX_IM 5
15 15 #define K15_SX_RE 6
16 16 #define K15_SX_IM 7
17 17 #define K24_SX_RE 8
18 18 #define K24_SX_IM 9
19 19 #define K25_SX_RE 10
20 20 #define K25_SX_IM 11
21 21 #define K34_SX_RE 12
22 22 #define K34_SX_IM 13
23 23 #define K35_SX_RE 14
24 24 #define K35_SX_IM 15
25 25
26 26 #define K24_NY_RE 16
27 27 #define K24_NY_IM 17
28 28 #define K25_NY_RE 18
29 29 #define K25_NY_IM 19
30 30 #define K34_NY_RE 20
31 31 #define K34_NY_IM 21
32 32 #define K35_NY_RE 22
33 33 #define K35_NY_IM 23
34 34
35 35 #define K24_NZ_RE 24
36 36 #define K24_NZ_IM 25
37 37 #define K25_NZ_RE 26
38 38 #define K25_NZ_IM 27
39 39 #define K34_NZ_RE 28
40 40 #define K34_NZ_IM 29
41 41 #define K35_NZ_RE 30
42 42 #define K35_NZ_IM 31
43 43
44 44 float k_f0[NB_BINS_COMPRESSED_MATRIX_f0][32];
45 45
46 void init_k_f0(){
46 void init_k_f0( void )
47 {
47 48 unsigned char i;
48 49
49 50 for(i=0; i<NB_BINS_COMPRESSED_MATRIX_f0; i++){
50 51 k_f0[i][K44_PE] = 1;
51 52 k_f0[i][K55_PE] = 1;
52 53 k_f0[i][K45_PE_RE] = 1;
53 54 k_f0[i][K45_PE_IM] = 1;
54 55
55 56 k_f0[i][K14_SX_RE] = 1;
56 57 k_f0[i][K14_SX_IM] = 1;
57 58 k_f0[i][K15_SX_RE] = 1;
58 59 k_f0[i][K15_SX_IM] = 1;
59 60 k_f0[i][K24_SX_RE] = 1;
60 61 k_f0[i][K24_SX_IM] = 1;
61 62 k_f0[i][K25_SX_RE] = 1;
62 63 k_f0[i][K25_SX_IM] = 1;
63 64 k_f0[i][K34_SX_RE] = 1;
64 65 k_f0[i][K34_SX_IM] = 1;
65 66 k_f0[i][K35_SX_RE] = 1;
66 67 k_f0[i][K35_SX_IM] = 1;
67 68
68 69 k_f0[i][K24_NY_RE] = 1;
69 70 k_f0[i][K24_NY_IM] = 1;
70 71 k_f0[i][K25_NY_RE] = 1;
71 72 k_f0[i][K25_NY_IM] = 1;
72 73 k_f0[i][K34_NY_RE] = 1;
73 74 k_f0[i][K34_NY_IM] = 1;
74 75 k_f0[i][K35_NY_RE] = 1;
75 76 k_f0[i][K35_NY_IM] = 1;
76 77
77 78 k_f0[i][K24_NZ_RE] = 1;
78 79 k_f0[i][K24_NZ_IM] = 1;
79 80 k_f0[i][K25_NZ_RE] = 1;
80 81 k_f0[i][K25_NZ_IM] = 1;
81 82 k_f0[i][K34_NZ_RE] = 1;
82 83 k_f0[i][K34_NZ_IM] = 1;
83 84 k_f0[i][K35_NZ_RE] = 1;
84 85 k_f0[i][K35_NZ_IM] = 1;
85 86 }
86 87 }
87 88
88 89 float alpha_M = 45 * (3.1415927/180);
89 90
90 void BP1_set(){
91 void BP1_set( float * compressed_spec_mat, unsigned char nb_bins_compressed_spec_mat, unsigned char * LFR_BP1 ){
91 92 int i, exponent;
92 float PSDB, PSDE, tmp, NVEC_V0, NVEC_V1, NVEC_V2, aux, tr_SB_SB,
93 e_cross_b_re, e_cross_b_im,
94 n_cross_e_scal_b_re, n_cross_e_scal_b_im,
95 ny, nz, bx_bx_star, vphi,
96 significand;
97 signed char nbitexp, nbitsig, expmin, expmax; // 8 bits
93 float PSDB;
94 float PSDE;
95 float tmp;
96 float NVEC_V0;
97 float NVEC_V1;
98 float NVEC_V2;
99 float aux;
100 float tr_SB_SB;
101 float e_cross_b_re;
102 float e_cross_b_im;
103 float n_cross_e_scal_b_re;
104 float n_cross_e_scal_b_im;
105 float ny;
106 float nz;
107 float bx_bx_star;
108 float vphi;
109 float significand;
110 signed char nbitexp;
111 signed char nbitsig;
112 signed char expmin;
113 signed char expmax; // 8 bits
98 114 short int rangesig; // 16 bits
99 unsigned short int psd, tmp_u_short_int; // 16 bits
115 unsigned short int psd;
116 unsigned short int tmp_u_short_int; // 16 bits
100 117 unsigned short int *pt_u_short_int; // pointer on unsigned 16-bit words
101 118 unsigned char tmp_u_char; // 8 bits
102 119 unsigned char *pt_u_char; // pointer on unsigned 8-bit bytes
103 120
104 121 init_k_f0();
105 122
106 123 #ifdef DEBUG_TCH
107 printf("Number of bins: %d\n", NB_BINS_COMPRESSED_MATRIX_f0);
124 printf("Number of bins: %d\n", nb_bins_compressed_spec_mat);
108 125 printf("BP1 : \n");
109 126 #endif
110 127
111 128 // initialization for managing the exponents of the floating point data:
112 129 nbitexp = 5; // number of bits for the exponent
113 130 expmax = 30; // maximum value of the exponent
114 131 expmin = expmax - (1 << nbitexp) + 1; // accordingly the minimum exponent value
115 132 // for floating point data to be recorded on 12-bit words:
116 133 nbitsig = 12 - nbitexp; // number of bits for the significand
117 134 rangesig = (1 << nbitsig)-1; // == 2^nbitsig - 1
118 135
119 136 #ifdef DEBUG_TCH
120 137 printf("nbitexp : %d, expmax : %d, expmin : %d\n", nbitexp, expmax, expmin);
121 138 printf("nbitsig : %d, rangesig : %d\n", nbitsig, rangesig);
122 139 #endif
123 140
124 for(i=0; i<NB_BINS_COMPRESSED_MATRIX_f0; i++){
141 for(i=0; i<nb_bins_compressed_spec_mat; i++){
125 142 //==============================================
126 143 // BP1 PSDB == PA_LFR_SC_BP1_PB_F0 == 12 bits = 5 bits (exponent) + 7 bits (significand)
127 PSDB = compressed_spectral_matrix_f0[i*30] // S11
128 + compressed_spectral_matrix_f0[i*30+9] // S22
129 + compressed_spectral_matrix_f0[i*30+16]; // S33
144 PSDB = compressed_spec_mat[i*30] // S11
145 + compressed_spec_mat[i*30+9] // S22
146 + compressed_spec_mat[i*30+16]; // S33
130 147
131 148 significand = frexpf(PSDB/3, &exponent); // 0.5 <= significand < 1
132 149 // PSDB/3 = significand * 2^exponent
133 150 // the division by 3 is to ensure that max value <= 2^30
134 151
135 152 if (exponent < expmin) { // value should be >= 0.5 * 2^expmin
136 153 exponent = expmin;
137 154 significand = 0.5; // min value that can be recorded
138 155 }
139 156 if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)
140 157 exponent = expmax;
141 158 significand = 1.0; // max value that can be recorded
142 159 }
143 160 if (significand == 0) {// in that case exponent == 0 too
144 161 exponent = expmin;
145 162 significand = 0.5; // min value that can be recorded
146 163 }
147 164
148 165 psd = (unsigned short int) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding
149 166 // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)
150 167 tmp_u_short_int = (unsigned short int) (exponent-expmin); // Shift and cast into a 16-bit unsigned int
151 168 // where just the first nbitexp bits are used (0, ..., 2^nbitexp-1)
152 pt_u_short_int = (unsigned short int*) &LFR_BP1_F0[i*9+2]; // Affect an unsigned short int pointer with the
169 pt_u_short_int = (unsigned short int*) &LFR_BP1[i*9+2]; // Affect an unsigned short int pointer with the
153 170 // adress where the 16-bit word result will be stored
154 171 *pt_u_short_int = psd | (tmp_u_short_int << nbitsig); // Put the exponent bits (nbitexp) next to the
155 172 // left place of the significand bits (nbitsig), making
156 173 // the 16-bit word to be recorded, and record it using the pointer
157 174 #ifdef DEBUG_TCH
158 175 printf("PSDB / 3 : %16.8e\n",PSDB/3);
159 176 printf("significand : %16.8e\n",significand);
160 177 printf("exponent : %d\n" ,exponent);
161 178 printf("psd for PSDB significand : %d\n",psd);
162 179 printf("tmp_u_short_int for PSDB exponent : %d\n",tmp_u_short_int);
163 180 printf("*pt_u_short_int for PSDB exponent + significand: %.3d or %.4x\n",*pt_u_short_int, *pt_u_short_int);
164 printf("LFR_BP1_F0[i*9+3] : %.3d or %.2x\n",LFR_BP1_F0[i*9+3], LFR_BP1_F0[i*9+3]);
165 printf("LFR_BP1_F0[i*9+2] : %.3d or %.2x\n",LFR_BP1_F0[i*9+2], LFR_BP1_F0[i*9+2]);
181 printf("LFR_BP1[i*9+3] : %.3d or %.2x\n",LFR_BP1[i*9+3], LFR_BP1[i*9+3]);
182 printf("LFR_BP1[i*9+2] : %.3d or %.2x\n",LFR_BP1[i*9+2], LFR_BP1[i*9+2]);
166 183 #endif
167 184 //==============================================
168 185 // BP1 PSDE == PA_LFR_SC_BP1_PE_F0 == 12 bits = 5 bits (exponent) + 7 bits (significand)
169 PSDE = compressed_spectral_matrix_f0[i*30+21] * k_f0[i][K44_PE] // S44
170 + compressed_spectral_matrix_f0[i*30+24] * k_f0[i][K55_PE] // S55
171 + compressed_spectral_matrix_f0[i*30+22] * k_f0[i][K45_PE_RE] // S45 Re
172 - compressed_spectral_matrix_f0[i*30+23] * k_f0[i][K45_PE_IM]; // S45 Im
186 PSDE = compressed_spec_mat[i*30+21] * k_f0[i][K44_PE] // S44
187 + compressed_spec_mat[i*30+24] * k_f0[i][K55_PE] // S55
188 + compressed_spec_mat[i*30+22] * k_f0[i][K45_PE_RE] // S45 Re
189 - compressed_spec_mat[i*30+23] * k_f0[i][K45_PE_IM]; // S45 Im
173 190
174 191 significand = frexpf(PSDE/2, &exponent); // 0.5 <= significand < 1
175 192 // PSDE/2 = significand * 2^exponent
176 193 // the division by 2 is to ensure that max value <= 2^30
177 194 // should be reconsidered by taking into account the k-coefficients ...
178 195
179 196 if (exponent < expmin) { // value should be >= 0.5 * 2^expmin
180 197 exponent = expmin;
181 198 significand = 0.5; // min value that can be recorded
182 199 }
183 200 if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)
184 201 exponent = expmax;
185 202 significand = 1.0; // max value that can be recorded
186 203 }
187 204 if (significand == 0) {// in that case exponent == 0 too
188 205 exponent = expmin;
189 206 significand = 0.5; // min value that can be recorded
190 207 }
191 208
192 209 psd = (unsigned short int) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding
193 210 // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)
194 211 tmp_u_short_int = (unsigned short int) (exponent-expmin); // Shift and cast into a 16-bit unsigned int
195 212 // where just the first nbitexp bits are used (0, ..., 2^nbitexp-1)
196 pt_u_short_int = (unsigned short int*) &LFR_BP1_F0[i*9+0]; // Affect an unsigned short int pointer with the
213 pt_u_short_int = (unsigned short int*) &LFR_BP1[i*9+0]; // Affect an unsigned short int pointer with the
197 214 // adress where the 16-bit word result will be stored
198 215 *pt_u_short_int = psd | (tmp_u_short_int << nbitsig); // Put the exponent bits (nbitexp) next to the
199 216 // left place of the significand bits (nbitsig), making
200 217 // the 16-bit word to be recorded, and record it using the pointer
201 218 #ifdef DEBUG_TCH
202 219 printf("PSDE / 2 : %16.8e\n",PSDE/2);
203 220 printf("significand : %16.8e\n",significand);
204 221 printf("exponent : %d\n" ,exponent);
205 222 printf("psd for PSDE significand : %d\n",psd);
206 223 printf("tmp_u_short_int for PSDE exponent : %d\n",tmp_u_short_int);
207 224 printf("*pt_u_short_int for PSDE exponent + significand: %.3d or %.4x\n",*pt_u_short_int, *pt_u_short_int);
208 printf("LFR_BP1_F0[i*9+1] : %.3d or %.2x\n",LFR_BP1_F0[i*9+1], LFR_BP1_F0[i*9+1]);
209 printf("LFR_BP1_F0[i*9+0] : %.3d or %.2x\n",LFR_BP1_F0[i*9+0], LFR_BP1_F0[i*9+0]);
225 printf("LFR_BP1[i*9+1] : %.3d or %.2x\n",LFR_BP1[i*9+1], LFR_BP1[i*9+1]);
226 printf("LFR_BP1[i*9+0] : %.3d or %.2x\n",LFR_BP1[i*9+0], LFR_BP1[i*9+0]);
210 227 #endif
211 228 //==============================================================================
212 229 // BP1 normal wave vector == PA_LFR_SC_BP1_NVEC_V0_F0 == 8 bits
213 230 // == PA_LFR_SC_BP1_NVEC_V1_F0 == 8 bits
214 231 // == PA_LFR_SC_BP1_NVEC_V2_F0 == 1 sign bit
215 tmp = sqrt( compressed_spectral_matrix_f0[i*30+2] *compressed_spectral_matrix_f0[i*30+2] //Im S12
216 +compressed_spectral_matrix_f0[i*30+4] *compressed_spectral_matrix_f0[i*30+4] //Im S13
217 +compressed_spectral_matrix_f0[i*30+11]*compressed_spectral_matrix_f0[i*30+11] //Im S23
232 tmp = sqrt( compressed_spec_mat[i*30+2] *compressed_spec_mat[i*30+2] //Im S12
233 +compressed_spec_mat[i*30+4] *compressed_spec_mat[i*30+4] //Im S13
234 +compressed_spec_mat[i*30+11]*compressed_spec_mat[i*30+11] //Im S23
218 235 );
219 NVEC_V0 = compressed_spectral_matrix_f0[i*30+11]/ tmp; // S23 Im => n1
220 NVEC_V1 = -compressed_spectral_matrix_f0[i*30+4] / tmp; // S13 Im => n2
221 NVEC_V2 = compressed_spectral_matrix_f0[i*30+2] / tmp; // S12 Im => n3
236 NVEC_V0 = compressed_spec_mat[i*30+11]/ tmp; // S23 Im => n1
237 NVEC_V1 = -compressed_spec_mat[i*30+4] / tmp; // S13 Im => n2
238 NVEC_V2 = compressed_spec_mat[i*30+2] / tmp; // S12 Im => n3
222 239
223 LFR_BP1_F0[i*9+4] = (unsigned char) (NVEC_V0*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
224 LFR_BP1_F0[i*9+5] = (unsigned char) (NVEC_V1*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
240 LFR_BP1[i*9+4] = (unsigned char) (NVEC_V0*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
241 LFR_BP1[i*9+5] = (unsigned char) (NVEC_V1*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
225 242 pt_u_char = (unsigned char*) &NVEC_V2; // affect an unsigned char pointer with the adress of NVEC_V2
226 243 #ifdef LSB_FIRST_TCH
227 LFR_BP1_F0[i*9+6] = pt_u_char[3] & 0x80; // extract the sign bit of NVEC_V2 (32-bit float, sign bit in the 4th octet:PC convention)
228 // record it at the 8th bit position (from the right to the left) of LFR_BP1_F0[i*9+6]
244 LFR_BP1[i*9+6] = pt_u_char[3] & 0x80; // extract the sign bit of NVEC_V2 (32-bit float, sign bit in the 4th octet:PC convention)
245 // record it at the 8th bit position (from the right to the left) of LFR_BP1[i*9+6]
229 246 #endif
230 247 #ifdef MSB_FIRST_TCH
231 LFR_BP1_F0[i*9+6] = pt_u_char[0] & 0x80; // extract the sign bit of NVEC_V2 (32-bit float, sign bit in the 0th octet:SPARC convention)
232 // record it at the 8th bit position (from the right to the left) of LFR_BP1_F0[i*9+6]
248 LFR_BP1[i*9+6] = pt_u_char[0] & 0x80; // extract the sign bit of NVEC_V2 (32-bit float, sign bit in the 0th octet:SPARC convention)
249 // record it at the 8th bit position (from the right to the left) of LFR_BP1[i*9+6]
233 250 #endif
234 251 #ifdef DEBUG_TCH
235 252 printf("NVEC_V0 : %16.8e\n",NVEC_V0);
236 253 printf("NVEC_V1 : %16.8e\n",NVEC_V1);
237 254 printf("NVEC_V2 : %16.8e\n",NVEC_V2);
238 printf("LFR_BP1_F0[i*9+4] for NVEC_V0 : %u\n",LFR_BP1_F0[i*9+4]);
239 printf("LFR_BP1_F0[i*9+5] for NVEC_V1 : %u\n",LFR_BP1_F0[i*9+5]);
240 printf("LFR_BP1_F0[i*9+6] for NVEC_V2 : %u\n",LFR_BP1_F0[i*9+6]);
255 printf("LFR_BP1[i*9+4] for NVEC_V0 : %u\n",LFR_BP1[i*9+4]);
256 printf("LFR_BP1[i*9+5] for NVEC_V1 : %u\n",LFR_BP1[i*9+5]);
257 printf("LFR_BP1[i*9+6] for NVEC_V2 : %u\n",LFR_BP1[i*9+6]);
241 258 #endif
242 259 //=======================================================
243 260 // BP1 ellipticity == PA_LFR_SC_BP1_ELLIP_F0 == 4 bits
244 261 aux = 2*tmp / PSDB; // compute the ellipticity
245 262
246 263 tmp_u_char = (unsigned char) (aux*15 + 0.5); // shift and cast into a 8-bit unsigned char with rounding
247 264 // where just the first 4 bits are used (0, ..., 15)
248 LFR_BP1_F0[i*9+6] = LFR_BP1_F0[i*9+6] | (tmp_u_char << 3); // put these 4 bits next to the right place
265 LFR_BP1[i*9+6] = LFR_BP1[i*9+6] | (tmp_u_char << 3); // put these 4 bits next to the right place
249 266 // of the sign bit of NVEC_V2 (recorded
250 // previously in LFR_BP1_F0[i*9+6])
267 // previously in LFR_BP1[i*9+6])
251 268 #ifdef DEBUG_TCH
252 269 printf("ellipticity : %16.8e\n",aux);
253 270 printf("tmp_u_char for ellipticity : %u\n",tmp_u_char);
254 printf("LFR_BP1_F0[i*9+6] for NVEC_V2 + ellipticity : %u\n",LFR_BP1_F0[i*9+6]);
271 printf("LFR_BP1[i*9+6] for NVEC_V2 + ellipticity : %u\n",LFR_BP1[i*9+6]);
255 272 #endif
256 273 //==============================================================
257 274 // BP1 degree of polarization == PA_LFR_SC_BP1_DOP_F0 == 3 bits
258 tr_SB_SB = compressed_spectral_matrix_f0[i*30] *compressed_spectral_matrix_f0[i*30]
259 + compressed_spectral_matrix_f0[i*30+9] *compressed_spectral_matrix_f0[i*30+9]
260 + compressed_spectral_matrix_f0[i*30+16] *compressed_spectral_matrix_f0[i*30+16]
261 + 2 * compressed_spectral_matrix_f0[i*30+1] *compressed_spectral_matrix_f0[i*30+1]
262 + 2 * compressed_spectral_matrix_f0[i*30+2] *compressed_spectral_matrix_f0[i*30+2]
263 + 2 * compressed_spectral_matrix_f0[i*30+3] *compressed_spectral_matrix_f0[i*30+3]
264 + 2 * compressed_spectral_matrix_f0[i*30+4] *compressed_spectral_matrix_f0[i*30+4]
265 + 2 * compressed_spectral_matrix_f0[i*30+10]*compressed_spectral_matrix_f0[i*30+10]
266 + 2 * compressed_spectral_matrix_f0[i*30+11]*compressed_spectral_matrix_f0[i*30+11];
275 tr_SB_SB = compressed_spec_mat[i*30] *compressed_spec_mat[i*30]
276 + compressed_spec_mat[i*30+9] *compressed_spec_mat[i*30+9]
277 + compressed_spec_mat[i*30+16] *compressed_spec_mat[i*30+16]
278 + 2 * compressed_spec_mat[i*30+1] *compressed_spec_mat[i*30+1]
279 + 2 * compressed_spec_mat[i*30+2] *compressed_spec_mat[i*30+2]
280 + 2 * compressed_spec_mat[i*30+3] *compressed_spec_mat[i*30+3]
281 + 2 * compressed_spec_mat[i*30+4] *compressed_spec_mat[i*30+4]
282 + 2 * compressed_spec_mat[i*30+10]*compressed_spec_mat[i*30+10]
283 + 2 * compressed_spec_mat[i*30+11]*compressed_spec_mat[i*30+11];
267 284 aux = PSDB*PSDB;
268 285 tmp = ( 3*tr_SB_SB - aux ) / ( 2 * aux ); // compute the degree of polarisation
269 286
270 287 tmp_u_char = (unsigned char) (tmp*7 + 0.5);// shift and cast into a 8-bit unsigned char with rounding
271 288 // where just the first 3 bits are used (0, ..., 7)
272 LFR_BP1_F0[i*9+6] = LFR_BP1_F0[i*9+6] | tmp_u_char; // record these 3 bits at the 3 first bit positions
273 // (from the right to the left) of LFR_BP1_F0[i*9+6]
289 LFR_BP1[i*9+6] = LFR_BP1[i*9+6] | tmp_u_char; // record these 3 bits at the 3 first bit positions
290 // (from the right to the left) of LFR_BP1[i*9+6]
274 291 #ifdef DEBUG_TCH
275 292 printf("DOP : %16.8e\n",tmp);
276 293 printf("tmp_u_char for DOP : %u\n",tmp_u_char);
277 printf("LFR_BP1_F0[i*9+6] for NVEC_V2 + ellipticity + DOP : %u\n",LFR_BP1_F0[i*9+6]);
294 printf("LFR_BP1[i*9+6] for NVEC_V2 + ellipticity + DOP : %u\n",LFR_BP1[i*9+6]);
278 295 #endif
279 296 //=======================================================================================
280 297 // BP1 X_SO-component of the Poynting flux == PA_LFR_SC_BP1_SX_F0 == 8 (+ 2) bits
281 298 // = 5 bits (exponent) + 3 bits (significand)
282 299 // + 1 sign bit + 1 argument bit (two sectors)
283 e_cross_b_re = compressed_spectral_matrix_f0[i*30+17]*k_f0[i][K34_SX_RE] //S34 Re
284 + compressed_spectral_matrix_f0[i*30+19]*k_f0[i][K35_SX_RE] //S35 Re
285 + compressed_spectral_matrix_f0[i*30+5] *k_f0[i][K14_SX_RE] //S14 Re
286 + compressed_spectral_matrix_f0[i*30+7] *k_f0[i][K15_SX_RE] //S15 Re
287 + compressed_spectral_matrix_f0[i*30+12]*k_f0[i][K24_SX_RE] //S24 Re
288 + compressed_spectral_matrix_f0[i*30+14]*k_f0[i][K25_SX_RE] //S25 Re
289 + compressed_spectral_matrix_f0[i*30+18]*k_f0[i][K34_SX_IM] //S34 Im
290 + compressed_spectral_matrix_f0[i*30+20]*k_f0[i][K35_SX_IM] //S35 Im
291 + compressed_spectral_matrix_f0[i*30+6] *k_f0[i][K14_SX_IM] //S14 Im
292 + compressed_spectral_matrix_f0[i*30+8] *k_f0[i][K15_SX_IM] //S15 Im
293 + compressed_spectral_matrix_f0[i*30+13]*k_f0[i][K24_SX_IM] //S24 Im
294 + compressed_spectral_matrix_f0[i*30+15]*k_f0[i][K25_SX_IM]; //S25 Im
300 e_cross_b_re = compressed_spec_mat[i*30+17]*k_f0[i][K34_SX_RE] //S34 Re
301 + compressed_spec_mat[i*30+19]*k_f0[i][K35_SX_RE] //S35 Re
302 + compressed_spec_mat[i*30+5] *k_f0[i][K14_SX_RE] //S14 Re
303 + compressed_spec_mat[i*30+7] *k_f0[i][K15_SX_RE] //S15 Re
304 + compressed_spec_mat[i*30+12]*k_f0[i][K24_SX_RE] //S24 Re
305 + compressed_spec_mat[i*30+14]*k_f0[i][K25_SX_RE] //S25 Re
306 + compressed_spec_mat[i*30+18]*k_f0[i][K34_SX_IM] //S34 Im
307 + compressed_spec_mat[i*30+20]*k_f0[i][K35_SX_IM] //S35 Im
308 + compressed_spec_mat[i*30+6] *k_f0[i][K14_SX_IM] //S14 Im
309 + compressed_spec_mat[i*30+8] *k_f0[i][K15_SX_IM] //S15 Im
310 + compressed_spec_mat[i*30+13]*k_f0[i][K24_SX_IM] //S24 Im
311 + compressed_spec_mat[i*30+15]*k_f0[i][K25_SX_IM]; //S25 Im
295 312 // Im(S_ji) = -Im(S_ij)
296 313 // k_ji = k_ij
297 e_cross_b_im = compressed_spectral_matrix_f0[i*30+17]*k_f0[i][K34_SX_IM] //S34 Re
298 + compressed_spectral_matrix_f0[i*30+19]*k_f0[i][K35_SX_IM] //S35 Re
299 + compressed_spectral_matrix_f0[i*30+5] *k_f0[i][K14_SX_IM] //S14 Re
300 + compressed_spectral_matrix_f0[i*30+7] *k_f0[i][K15_SX_IM] //S15 Re
301 + compressed_spectral_matrix_f0[i*30+12]*k_f0[i][K24_SX_IM] //S24 Re
302 + compressed_spectral_matrix_f0[i*30+14]*k_f0[i][K25_SX_IM] //S25 Re
303 - compressed_spectral_matrix_f0[i*30+18]*k_f0[i][K34_SX_RE] //S34 Im
304 - compressed_spectral_matrix_f0[i*30+20]*k_f0[i][K35_SX_RE] //S35 Im
305 - compressed_spectral_matrix_f0[i*30+6] *k_f0[i][K14_SX_RE] //S14 Im
306 - compressed_spectral_matrix_f0[i*30+8] *k_f0[i][K15_SX_RE] //S15 Im
307 - compressed_spectral_matrix_f0[i*30+13]*k_f0[i][K24_SX_RE] //S24 Im
308 - compressed_spectral_matrix_f0[i*30+15]*k_f0[i][K25_SX_RE]; //S25 Im
314 e_cross_b_im = compressed_spec_mat[i*30+17]*k_f0[i][K34_SX_IM] //S34 Re
315 + compressed_spec_mat[i*30+19]*k_f0[i][K35_SX_IM] //S35 Re
316 + compressed_spec_mat[i*30+5] *k_f0[i][K14_SX_IM] //S14 Re
317 + compressed_spec_mat[i*30+7] *k_f0[i][K15_SX_IM] //S15 Re
318 + compressed_spec_mat[i*30+12]*k_f0[i][K24_SX_IM] //S24 Re
319 + compressed_spec_mat[i*30+14]*k_f0[i][K25_SX_IM] //S25 Re
320 - compressed_spec_mat[i*30+18]*k_f0[i][K34_SX_RE] //S34 Im
321 - compressed_spec_mat[i*30+20]*k_f0[i][K35_SX_RE] //S35 Im
322 - compressed_spec_mat[i*30+6] *k_f0[i][K14_SX_RE] //S14 Im
323 - compressed_spec_mat[i*30+8] *k_f0[i][K15_SX_RE] //S15 Im
324 - compressed_spec_mat[i*30+13]*k_f0[i][K24_SX_RE] //S24 Im
325 - compressed_spec_mat[i*30+15]*k_f0[i][K25_SX_RE]; //S25 Im
309 326 #ifdef DEBUG_TCH
310 327 printf("ReaSX / 2 : %16.8e\n",e_cross_b_re/2);
311 328 #endif
312 329 pt_u_char = (unsigned char*) &e_cross_b_re; // Affect an unsigned char pointer with the adress of e_cross_b_re
313 330 #ifdef LSB_FIRST_TCH
314 LFR_BP1_F0[i*9+1] = LFR_BP1_F0[i*9+1] | (pt_u_char[3] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 4th octet:PC convention)
331 LFR_BP1[i*9+1] = LFR_BP1[i*9+1] | (pt_u_char[3] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 4th octet:PC convention)
315 332 // Record it at the 8th bit position (from the right to the left)
316 // of LFR_BP1_F0[i*9+1]
333 // of LFR_BP1[i*9+1]
317 334 pt_u_char[3] = (pt_u_char[3] & 0x7f); // Make e_cross_b_re be positive in any case: |ReaSX|
318 335 #endif
319 336 #ifdef MSB_FIRST_TCH
320 LFR_BP1_F0[i*9+1] = LFR_BP1_F0[i*9+1] | (pt_u_char[0] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 0th octet:SPARC convention)
337 LFR_BP1[i*9+1] = LFR_BP1[i*9+1] | (pt_u_char[0] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 0th octet:SPARC convention)
321 338 // Record it at the 8th bit position (from the right to the left)
322 // of LFR_BP1_F0[i*9+1]
339 // of LFR_BP1[i*9+1]
323 340 pt_u_char[0] = (pt_u_char[0] & 0x7f); // Make e_cross_b_re be positive in any case: |ReaSX|
324 341 #endif
325 342 significand = frexpf(e_cross_b_re/2, &exponent);// 0.5 <= significand < 1
326 343 // ReaSX/2 = significand * 2^exponent
327 344 // The division by 2 is to ensure that max value <= 2^30 (rough estimate)
328 345 // Should be reconsidered by taking into account the k-coefficients ...
329 346
330 347 if (exponent < expmin) { // value should be >= 0.5 * 2^expmin
331 348 exponent = expmin;
332 349 significand = 0.5; // min value that can be recorded
333 350 }
334 351 if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)
335 352 exponent = expmax;
336 353 significand = 1.0; // max value that can be recorded
337 354 }
338 355 if (significand == 0) {// in that case exponent == 0 too
339 356 exponent = expmin;
340 357 significand = 0.5; // min value that can be recorded
341 358 }
342 359
343 LFR_BP1_F0[i*9+7] = (unsigned char) ((significand*2-1)*7 + 0.5); // Shift and cast into a 8-bit unsigned char with rounding
360 LFR_BP1[i*9+7] = (unsigned char) ((significand*2-1)*7 + 0.5); // Shift and cast into a 8-bit unsigned char with rounding
344 361 // where just the first 3 bits are used (0, ..., 7)
345 362 tmp_u_char = (unsigned char) (exponent-expmin); // Shift and cast into a 8-bit unsigned char where
346 363 // just the first 5 bits are used (0, ..., 2^5-1)
347 364 #ifdef DEBUG_TCH
348 365 printf("|ReaSX| / 2 : %16.8e\n",e_cross_b_re/2);
349 366 printf("significand : %16.8e\n",significand);
350 367 printf("exponent : %d\n" ,exponent);
351 printf("LFR_BP1_F0[i*9+7] for ReaSX significand : %u\n",LFR_BP1_F0[i*9+7]);
368 printf("LFR_BP1[i*9+7] for ReaSX significand : %u\n",LFR_BP1[i*9+7]);
352 369 printf("tmp_u_char for ReaSX exponent : %d\n",tmp_u_char);
353 370 #endif
354 LFR_BP1_F0[i*9+7] = LFR_BP1_F0[i*9+7] | (tmp_u_char << 3); // shift these 5 bits to the left before logical addition
355 // with LFR_BP1_F0[i*9+7]
371 LFR_BP1[i*9+7] = LFR_BP1[i*9+7] | (tmp_u_char << 3); // shift these 5 bits to the left before logical addition
372 // with LFR_BP1[i*9+7]
356 373 #ifdef DEBUG_TCH
357 printf("LFR_BP1_F0[i*9+7] for ReaSX exponent + significand : %u\n",LFR_BP1_F0[i*9+7]);
358 printf("LFR_BP1_F0[i*9+1] for ReaSX sign + PSDE 'exponent' : %u\n",LFR_BP1_F0[i*9+1]);
374 printf("LFR_BP1[i*9+7] for ReaSX exponent + significand : %u\n",LFR_BP1[i*9+7]);
375 printf("LFR_BP1[i*9+1] for ReaSX sign + PSDE 'exponent' : %u\n",LFR_BP1[i*9+1]);
359 376 printf("ImaSX / 2 : %16.8e\n",e_cross_b_im/2);
360 377 #endif
361 378 pt_u_char = (unsigned char*) &e_cross_b_im; // Affect an unsigned char pointer with the adress of e_cross_b_im
362 379 #ifdef LSB_FIRST_TCH
363 380 pt_u_char[3] = pt_u_char[3] & 0x7f; // Make e_cross_b_im be positive in any case: |ImaSX|
364 381 #endif
365 382 #ifdef MSB_FIRST_TCH
366 383 pt_u_char[0] = pt_u_char[0] & 0x7f; // Make e_cross_b_im be positive in any case: |ImaSX|
367 384 #endif
368 385 tmp_u_char = (e_cross_b_im > e_cross_b_re) ? 0x40 : 0x00; // Determine the sector argument of SX. If |Im| > |Re| affect
369 386 // an unsigned 8-bit char with 01000000; otherwise with null.
370 LFR_BP1_F0[i*9+1] = LFR_BP1_F0[i*9+1] | tmp_u_char; // Record it as a sign bit at the 7th bit position (from the right
371 // to the left) of LFR_BP1_F0[i*9+1], by simple logical addition.
387 LFR_BP1[i*9+1] = LFR_BP1[i*9+1] | tmp_u_char; // Record it as a sign bit at the 7th bit position (from the right
388 // to the left) of LFR_BP1[i*9+1], by simple logical addition.
372 389 #ifdef DEBUG_TCH
373 390 printf("|ImaSX| / 2 : %16.8e\n",e_cross_b_im/2);
374 391 printf("ArgSX sign : %u\n",tmp_u_char);
375 printf("LFR_BP1_F0[i*9+1] for ReaSX & ArgSX signs + PSDE 'exponent' : %u\n",LFR_BP1_F0[i*9+1]);
392 printf("LFR_BP1[i*9+1] for ReaSX & ArgSX signs + PSDE 'exponent' : %u\n",LFR_BP1[i*9+1]);
376 393 #endif
377 394 //======================================================================
378 395 // BP1 phase velocity estimator == PA_LFR_SC_BP1_VPHI_F0 == 8 (+ 2) bits
379 396 // = 5 bits (exponent) + 3 bits (significand)
380 397 // + 1 sign bit + 1 argument bit (two sectors)
381 398 ny = sin(alpha_M)*NVEC_V1 + cos(alpha_M)*NVEC_V2;
382 399 nz = NVEC_V0;
383 bx_bx_star = cos(alpha_M)*cos(alpha_M)*compressed_spectral_matrix_f0[i*30+9] // S22 Re
384 + sin(alpha_M)*sin(alpha_M)*compressed_spectral_matrix_f0[i*30+16] // S33 Re
385 - 2*sin(alpha_M)*cos(alpha_M)*compressed_spectral_matrix_f0[i*30+10]; // S23 Re
400 bx_bx_star = cos(alpha_M)*cos(alpha_M)*compressed_spec_mat[i*30+9] // S22 Re
401 + sin(alpha_M)*sin(alpha_M)*compressed_spec_mat[i*30+16] // S33 Re
402 - 2*sin(alpha_M)*cos(alpha_M)*compressed_spec_mat[i*30+10]; // S23 Re
386 403
387 n_cross_e_scal_b_re = ny * (compressed_spectral_matrix_f0[i*30+12]*k_f0[i][K24_NY_RE] //S24 Re
388 +compressed_spectral_matrix_f0[i*30+14]*k_f0[i][K25_NY_RE] //S25 Re
389 +compressed_spectral_matrix_f0[i*30+17]*k_f0[i][K34_NY_RE] //S34 Re
390 +compressed_spectral_matrix_f0[i*30+19]*k_f0[i][K35_NY_RE] //S35 Re
391 +compressed_spectral_matrix_f0[i*30+13]*k_f0[i][K24_NY_IM] //S24 Im
392 +compressed_spectral_matrix_f0[i*30+15]*k_f0[i][K25_NY_IM] //S25 Im
393 +compressed_spectral_matrix_f0[i*30+18]*k_f0[i][K34_NY_IM] //S34 Im
394 +compressed_spectral_matrix_f0[i*30+20]*k_f0[i][K35_NY_IM]) //S35 Im
395 + nz * (compressed_spectral_matrix_f0[i*30+12]*k_f0[i][K24_NZ_RE] //S24 Re
396 +compressed_spectral_matrix_f0[i*30+14]*k_f0[i][K25_NZ_RE] //S25 Re
397 +compressed_spectral_matrix_f0[i*30+17]*k_f0[i][K34_NZ_RE] //S34 Re
398 +compressed_spectral_matrix_f0[i*30+19]*k_f0[i][K35_NZ_RE] //S35 Re
399 +compressed_spectral_matrix_f0[i*30+13]*k_f0[i][K24_NZ_IM] //S24 Im
400 +compressed_spectral_matrix_f0[i*30+15]*k_f0[i][K25_NZ_IM] //S25 Im
401 +compressed_spectral_matrix_f0[i*30+18]*k_f0[i][K34_NZ_IM] //S34 Im
402 +compressed_spectral_matrix_f0[i*30+20]*k_f0[i][K35_NZ_IM]);//S35 Im
404 n_cross_e_scal_b_re = ny * (compressed_spec_mat[i*30+12]*k_f0[i][K24_NY_RE] //S24 Re
405 +compressed_spec_mat[i*30+14]*k_f0[i][K25_NY_RE] //S25 Re
406 +compressed_spec_mat[i*30+17]*k_f0[i][K34_NY_RE] //S34 Re
407 +compressed_spec_mat[i*30+19]*k_f0[i][K35_NY_RE] //S35 Re
408 +compressed_spec_mat[i*30+13]*k_f0[i][K24_NY_IM] //S24 Im
409 +compressed_spec_mat[i*30+15]*k_f0[i][K25_NY_IM] //S25 Im
410 +compressed_spec_mat[i*30+18]*k_f0[i][K34_NY_IM] //S34 Im
411 +compressed_spec_mat[i*30+20]*k_f0[i][K35_NY_IM]) //S35 Im
412 + nz * (compressed_spec_mat[i*30+12]*k_f0[i][K24_NZ_RE] //S24 Re
413 +compressed_spec_mat[i*30+14]*k_f0[i][K25_NZ_RE] //S25 Re
414 +compressed_spec_mat[i*30+17]*k_f0[i][K34_NZ_RE] //S34 Re
415 +compressed_spec_mat[i*30+19]*k_f0[i][K35_NZ_RE] //S35 Re
416 +compressed_spec_mat[i*30+13]*k_f0[i][K24_NZ_IM] //S24 Im
417 +compressed_spec_mat[i*30+15]*k_f0[i][K25_NZ_IM] //S25 Im
418 +compressed_spec_mat[i*30+18]*k_f0[i][K34_NZ_IM] //S34 Im
419 +compressed_spec_mat[i*30+20]*k_f0[i][K35_NZ_IM]);//S35 Im
403 420 // Im(S_ji) = -Im(S_ij)
404 421 // k_ji = k_ij
405 n_cross_e_scal_b_im = ny * (compressed_spectral_matrix_f0[i*30+12]*k_f0[i][K24_NY_IM] //S24 Re
406 +compressed_spectral_matrix_f0[i*30+14]*k_f0[i][K25_NY_IM] //S25 Re
407 +compressed_spectral_matrix_f0[i*30+17]*k_f0[i][K34_NY_IM] //S34 Re
408 +compressed_spectral_matrix_f0[i*30+19]*k_f0[i][K35_NY_IM] //S35 Re
409 -compressed_spectral_matrix_f0[i*30+13]*k_f0[i][K24_NY_RE] //S24 Im
410 -compressed_spectral_matrix_f0[i*30+15]*k_f0[i][K25_NY_RE] //S25 Im
411 -compressed_spectral_matrix_f0[i*30+18]*k_f0[i][K34_NY_RE] //S34 Im
412 -compressed_spectral_matrix_f0[i*30+20]*k_f0[i][K35_NY_RE]) //S35 Im
413 + nz * (compressed_spectral_matrix_f0[i*30+12]*k_f0[i][K24_NZ_IM] //S24 Re
414 +compressed_spectral_matrix_f0[i*30+14]*k_f0[i][K25_NZ_IM] //S25 Re
415 +compressed_spectral_matrix_f0[i*30+17]*k_f0[i][K34_NZ_IM] //S34 Re
416 +compressed_spectral_matrix_f0[i*30+19]*k_f0[i][K35_NZ_IM] //S35 Re
417 -compressed_spectral_matrix_f0[i*30+13]*k_f0[i][K24_NZ_RE] //S24 Im
418 -compressed_spectral_matrix_f0[i*30+15]*k_f0[i][K25_NZ_RE] //S25 Im
419 -compressed_spectral_matrix_f0[i*30+18]*k_f0[i][K34_NZ_RE] //S34 Im
420 -compressed_spectral_matrix_f0[i*30+20]*k_f0[i][K35_NZ_RE]);//S35 Im
422 n_cross_e_scal_b_im = ny * (compressed_spec_mat[i*30+12]*k_f0[i][K24_NY_IM] //S24 Re
423 +compressed_spec_mat[i*30+14]*k_f0[i][K25_NY_IM] //S25 Re
424 +compressed_spec_mat[i*30+17]*k_f0[i][K34_NY_IM] //S34 Re
425 +compressed_spec_mat[i*30+19]*k_f0[i][K35_NY_IM] //S35 Re
426 -compressed_spec_mat[i*30+13]*k_f0[i][K24_NY_RE] //S24 Im
427 -compressed_spec_mat[i*30+15]*k_f0[i][K25_NY_RE] //S25 Im
428 -compressed_spec_mat[i*30+18]*k_f0[i][K34_NY_RE] //S34 Im
429 -compressed_spec_mat[i*30+20]*k_f0[i][K35_NY_RE]) //S35 Im
430 + nz * (compressed_spec_mat[i*30+12]*k_f0[i][K24_NZ_IM] //S24 Re
431 +compressed_spec_mat[i*30+14]*k_f0[i][K25_NZ_IM] //S25 Re
432 +compressed_spec_mat[i*30+17]*k_f0[i][K34_NZ_IM] //S34 Re
433 +compressed_spec_mat[i*30+19]*k_f0[i][K35_NZ_IM] //S35 Re
434 -compressed_spec_mat[i*30+13]*k_f0[i][K24_NZ_RE] //S24 Im
435 -compressed_spec_mat[i*30+15]*k_f0[i][K25_NZ_RE] //S25 Im
436 -compressed_spec_mat[i*30+18]*k_f0[i][K34_NZ_RE] //S34 Im
437 -compressed_spec_mat[i*30+20]*k_f0[i][K35_NZ_RE]);//S35 Im
421 438 #ifdef DEBUG_TCH
422 439 printf("n_cross_e_scal_b_re : %16.8e\n",n_cross_e_scal_b_re);
423 440 printf("n_cross_e_scal_b_im : %16.8e\n",n_cross_e_scal_b_im);
424 441 #endif
425 442 // vphi = n_cross_e_scal_b_re / bx_bx_star => sign(VPHI) = sign(n_cross_e_scal_b_re)
426 443 pt_u_char = (unsigned char*) &n_cross_e_scal_b_re; // Affect an unsigned char pointer with the adress of n_cross_e_scal_b_re
427 444 #ifdef LSB_FIRST_TCH
428 LFR_BP1_F0[i*9+3] = LFR_BP1_F0[i*9+3] | (pt_u_char[3] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 4th octet:PC convention)
445 LFR_BP1[i*9+3] = LFR_BP1[i*9+3] | (pt_u_char[3] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 4th octet:PC convention)
429 446 // Record it at the 8th bit position (from the right to the left)
430 // of LFR_BP1_F0[i*9+3]
447 // of LFR_BP1[i*9+3]
431 448 pt_u_char[3] = (pt_u_char[3] & 0x7f); // Make n_cross_e_scal_b_re be positive in any case: |n_cross_e_scal_b_re|
432 449 #endif
433 450 #ifdef MSB_FIRST_TCH
434 LFR_BP1_F0[i*9+3] = LFR_BP1_F0[i*9+3] | (pt_u_char[0] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 0th octet:SPARC convention)
451 LFR_BP1[i*9+3] = LFR_BP1[i*9+3] | (pt_u_char[0] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 0th octet:SPARC convention)
435 452 // Record it at the 8th bit position (from the right to the left)
436 // of LFR_BP1_F0[i*9+3]
453 // of LFR_BP1[i*9+3]
437 454 pt_u_char[0] = (pt_u_char[0] & 0x7f); // Make n_cross_e_scal_b_re be positive in any case: |n_cross_e_scal_b_re|
438 455 #endif
439 456 vphi = n_cross_e_scal_b_re / bx_bx_star; // Compute |VPHI|
440 457
441 458 significand = frexpf(vphi/2, &exponent); // 0.5 <= significand < 1
442 459 // vphi/2 = significand * 2^exponent
443 460 // The division by 2 is to ensure that max value <= 2^30 (rough estimate)
444 461 // Should be reconsidered by taking into account the k-coefficients ...
445 462
446 463 if (exponent < expmin) { // value should be >= 0.5 * 2^expmin
447 464 exponent = expmin;
448 465 significand = 0.5; // min value that can be recorded
449 466 }
450 467 if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)
451 468 exponent = expmax;
452 469 significand = 1.0; // max value that can be recorded
453 470 }
454 471 if (significand == 0) {// in that case exponent == 0 too
455 472 exponent = expmin;
456 473 significand = 0.5; // min value that can be recorded
457 474 }
458 475 #ifdef DEBUG_TCH
459 476 printf("|VPHI| / 2 : %16.8e\n",vphi/2);
460 477 printf("significand : %16.8e\n",significand);
461 478 printf("exponent : %d\n" ,exponent);
462 479 #endif
463 LFR_BP1_F0[i*9+8] = (unsigned char) ((significand*2-1)*7 + 0.5); // Shift and cast into a 8-bit unsigned char with rounding
480 LFR_BP1[i*9+8] = (unsigned char) ((significand*2-1)*7 + 0.5); // Shift and cast into a 8-bit unsigned char with rounding
464 481 // where just the first 3 bits are used (0, ..., 7)
465 482 tmp_u_char = (unsigned char) (exponent-expmin); // Shift and cast into a 8-bit unsigned char where
466 483 // just the first 5 bits are used (0, ..., 2^5-1)
467 484 #ifdef DEBUG_TCH
468 printf("LFR_BP1_F0[i*9+8] for VPHI significand : %u\n",LFR_BP1_F0[i*9+8]);
485 printf("LFR_BP1[i*9+8] for VPHI significand : %u\n",LFR_BP1[i*9+8]);
469 486 printf("tmp_u_char for VPHI exponent : %d\n",tmp_u_char);
470 487 #endif
471 LFR_BP1_F0[i*9+8] = LFR_BP1_F0[i*9+8] | (tmp_u_char << 3); // shift these 5 bits to the left before logical addition
472 // with LFR_BP1_F0[i*9+8]
488 LFR_BP1[i*9+8] = LFR_BP1[i*9+8] | (tmp_u_char << 3); // shift these 5 bits to the left before logical addition
489 // with LFR_BP1[i*9+8]
473 490 #ifdef DEBUG_TCH
474 printf("LFR_BP1_F0[i*9+8] for VPHI exponent + significand : %u\n",LFR_BP1_F0[i*9+8]);
475 printf("LFR_BP1_F0[i*9+3] for VPHI sign + PSDB 'exponent' : %u\n",LFR_BP1_F0[i*9+3]);
491 printf("LFR_BP1[i*9+8] for VPHI exponent + significand : %u\n",LFR_BP1[i*9+8]);
492 printf("LFR_BP1[i*9+3] for VPHI sign + PSDB 'exponent' : %u\n",LFR_BP1[i*9+3]);
476 493 #endif
477 494 pt_u_char = (unsigned char*) &n_cross_e_scal_b_im; // Affect an unsigned char pointer with the adress of n_cross_e_scal_b_im
478 495 #ifdef LSB_FIRST_TCH
479 496 pt_u_char[3] = pt_u_char[3] & 0x7f; // Make n_cross_e_scal_b_im be positive in any case: |ImaSX|
480 497 #endif
481 498 #ifdef MSB_FIRST_TCH
482 499 pt_u_char[0] = pt_u_char[0] & 0x7f; // Make n_cross_e_scal_b_im be positive in any case: |ImaSX|
483 500 #endif
484 501 tmp_u_char = (n_cross_e_scal_b_im > n_cross_e_scal_b_re) ? 0x40 : 0x00; // Determine the sector argument of SX. If |Im| > |Re| affect
485 502 // an unsigned 8-bit char with 01000000; otherwise with null.
486 LFR_BP1_F0[i*9+3] = LFR_BP1_F0[i*9+3] | tmp_u_char; // Record it as a sign bit at the 7th bit position (from the right
487 // to the left) of LFR_BP1_F0[i*9+3], by simple logical addition.
503 LFR_BP1[i*9+3] = LFR_BP1[i*9+3] | tmp_u_char; // Record it as a sign bit at the 7th bit position (from the right
504 // to the left) of LFR_BP1[i*9+3], by simple logical addition.
488 505 #ifdef DEBUG_TCH
489 506 printf("|n_cross_e_scal_b_im| : %16.8e\n",n_cross_e_scal_b_im);
490 507 printf("|n_cross_e_scal_b_im|/bx_bx_star/2: %16.8e\n",n_cross_e_scal_b_im/bx_bx_star/2);
491 508 printf("ArgNEBX sign : %u\n",tmp_u_char);
492 printf("LFR_BP1_F0[i*9+3] for VPHI & ArgNEBX signs + PSDB 'exponent' : %u\n",LFR_BP1_F0[i*9+3]);
509 printf("LFR_BP1[i*9+3] for VPHI & ArgNEBX signs + PSDB 'exponent' : %u\n",LFR_BP1[i*9+3]);
493 510 #endif
494 511 }
495 512 }
496 513
497 void BP2_set(){
514 void BP2_set( float * compressed_spec_mat, unsigned char nb_bins_compressed_spec_mat, unsigned char * LFR_BP2 )
515 {
498 516 int i, exponent;
499 517 float aux, significand, cross_re, cross_im;
500 518 signed char nbitexp, nbitsig, expmin, expmax; // 8 bits
501 519 short int rangesig; // 16 bits
502 520 unsigned short int autocor, tmp_u_short_int; // 16 bits
503 521 unsigned short int *pt_u_short_int; // pointer on unsigned 16-bit words
504 522
505 523 #ifdef DEBUG_TCH
506 printf("Number of bins: %d\n", NB_BINS_COMPRESSED_MATRIX_f0);
524 printf("Number of bins: %d\n", nb_bins_compressed_spec_mat);
507 525 printf("BP2 : \n");
508 526 #endif
509 527
510 528 // For floating point data to be recorded on 16-bit words :
511 529 nbitexp = 6; // number of bits for the exponent
512 530 nbitsig = 16 - nbitexp; // number of bits for the significand
513 531 rangesig = (1 << nbitsig)-1; // == 2^nbitsig - 1
514 532 expmax = 32;
515 533 expmin = expmax - (1 << nbitexp) + 1;
516 534
517 535 #ifdef DEBUG_TCH
518 536 printf("nbitexp : %d, nbitsig : %d, rangesig : %d\n", nbitexp, nbitsig, rangesig);
519 537 printf("expmin : %d, expmax : %d\n", expmin, expmax);
520 538 #endif
521 539
522 for(i = 0; i<NB_BINS_COMPRESSED_MATRIX_f0; i++){
540 for(i = 0; i<nb_bins_compressed_spec_mat; i++){
523 541 //==============================================
524 542 // BP2 normalized cross correlations == PA_LFR_SC_BP2_CROSS_F0 == 10 * (8+8) bits
525 543 // == PA_LFR_SC_BP2_CROSS_RE_0_F0 == 8 bits
526 544 // == PA_LFR_SC_BP2_CROSS_IM_0_F0 == 8 bits
527 545 // == PA_LFR_SC_BP2_CROSS_RE_1_F0 == 8 bits
528 546 // == PA_LFR_SC_BP2_CROSS_IM_1_F0 == 8 bits
529 547 // == PA_LFR_SC_BP2_CROSS_RE_2_F0 == 8 bits
530 548 // == PA_LFR_SC_BP2_CROSS_IM_2_F0 == 8 bits
531 549 // == PA_LFR_SC_BP2_CROSS_RE_3_F0 == 8 bits
532 550 // == PA_LFR_SC_BP2_CROSS_IM_3_F0 == 8 bits
533 551 // == PA_LFR_SC_BP2_CROSS_RE_4_F0 == 8 bits
534 552 // == PA_LFR_SC_BP2_CROSS_IM_4_F0 == 8 bits
535 553 // == PA_LFR_SC_BP2_CROSS_RE_5_F0 == 8 bits
536 554 // == PA_LFR_SC_BP2_CROSS_IM_5_F0 == 8 bits
537 555 // == PA_LFR_SC_BP2_CROSS_RE_6_F0 == 8 bits
538 556 // == PA_LFR_SC_BP2_CROSS_IM_6_F0 == 8 bits
539 557 // == PA_LFR_SC_BP2_CROSS_RE_7_F0 == 8 bits
540 558 // == PA_LFR_SC_BP2_CROSS_IM_7_F0 == 8 bits
541 559 // == PA_LFR_SC_BP2_CROSS_RE_8_F0 == 8 bits
542 560 // == PA_LFR_SC_BP2_CROSS_IM_8_F0 == 8 bits
543 561 // == PA_LFR_SC_BP2_CROSS_RE_9_F0 == 8 bits
544 562 // == PA_LFR_SC_BP2_CROSS_IM_9_F0 == 8 bits
545 563 // S12
546 aux = sqrt(compressed_spectral_matrix_f0[i*30]*compressed_spectral_matrix_f0[i*30+9]);
547 cross_re = compressed_spectral_matrix_f0[i*30+1] / aux;
548 cross_im = compressed_spectral_matrix_f0[i*30+2] / aux;
549 LFR_BP2_F0[i*30+10] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
550 LFR_BP2_F0[i*30+20] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
564 aux = sqrt(compressed_spec_mat[i*30]*compressed_spec_mat[i*30+9]);
565 cross_re = compressed_spec_mat[i*30+1] / aux;
566 cross_im = compressed_spec_mat[i*30+2] / aux;
567 LFR_BP2[i*30+10] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
568 LFR_BP2[i*30+20] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
551 569 #ifdef DEBUG_TCH
552 printf("LFR_BP2_F0[i*30+10] for cross12_re (%16.8e) : %.3u\n",cross_re, LFR_BP2_F0[i*30+10]);
553 printf("LFR_BP2_F0[i*30+20] for cross12_im (%16.8e) : %.3u\n",cross_im, LFR_BP2_F0[i*30+20]);
570 printf("LFR_BP2[i*30+10] for cross12_re (%16.8e) : %.3u\n",cross_re, LFR_BP2[i*30+10]);
571 printf("LFR_BP2[i*30+20] for cross12_im (%16.8e) : %.3u\n",cross_im, LFR_BP2[i*30+20]);
554 572 #endif
555 573 // S13
556 aux = sqrt(compressed_spectral_matrix_f0[i*30]*compressed_spectral_matrix_f0[i*30+16]);
557 cross_re = compressed_spectral_matrix_f0[i*30+3] / aux;
558 cross_im = compressed_spectral_matrix_f0[i*30+4] / aux;
559 LFR_BP2_F0[i*30+11] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
560 LFR_BP2_F0[i*30+21] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
574 aux = sqrt(compressed_spec_mat[i*30]*compressed_spec_mat[i*30+16]);
575 cross_re = compressed_spec_mat[i*30+3] / aux;
576 cross_im = compressed_spec_mat[i*30+4] / aux;
577 LFR_BP2[i*30+11] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
578 LFR_BP2[i*30+21] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
561 579 #ifdef DEBUG_TCH
562 printf("LFR_BP2_F0[i*30+11] for cross13_re (%16.8e) : %.3u\n",cross_re, LFR_BP2_F0[i*30+11]);
563 printf("LFR_BP2_F0[i*30+21] for cross13_im (%16.8e) : %.3u\n",cross_im, LFR_BP2_F0[i*30+21]);
580 printf("LFR_BP2[i*30+11] for cross13_re (%16.8e) : %.3u\n",cross_re, LFR_BP2[i*30+11]);
581 printf("LFR_BP2[i*30+21] for cross13_im (%16.8e) : %.3u\n",cross_im, LFR_BP2[i*30+21]);
564 582 #endif
565 583 // S14
566 aux = sqrt(compressed_spectral_matrix_f0[i*30]*compressed_spectral_matrix_f0[i*30+21]);
567 cross_re = compressed_spectral_matrix_f0[i*30+5] / aux;
568 cross_im = compressed_spectral_matrix_f0[i*30+6] / aux;
569 LFR_BP2_F0[i*30+12] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
570 LFR_BP2_F0[i*30+22] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
584 aux = sqrt(compressed_spec_mat[i*30]*compressed_spec_mat[i*30+21]);
585 cross_re = compressed_spec_mat[i*30+5] / aux;
586 cross_im = compressed_spec_mat[i*30+6] / aux;
587 LFR_BP2[i*30+12] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
588 LFR_BP2[i*30+22] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
571 589 #ifdef DEBUG_TCH
572 printf("LFR_BP2_F0[i*30+12] for cross14_re (%16.8e) : %.3u\n",cross_re, LFR_BP2_F0[i*30+12]);
573 printf("LFR_BP2_F0[i*30+22] for cross14_im (%16.8e) : %.3u\n",cross_im, LFR_BP2_F0[i*30+22]);
590 printf("LFR_BP2[i*30+12] for cross14_re (%16.8e) : %.3u\n",cross_re, LFR_BP2[i*30+12]);
591 printf("LFR_BP2[i*30+22] for cross14_im (%16.8e) : %.3u\n",cross_im, LFR_BP2[i*30+22]);
574 592 #endif
575 593 // S15
576 aux = sqrt(compressed_spectral_matrix_f0[i*30]*compressed_spectral_matrix_f0[i*30+24]);
577 cross_re = compressed_spectral_matrix_f0[i*30+7] / aux;
578 cross_im = compressed_spectral_matrix_f0[i*30+8] / aux;
579 LFR_BP2_F0[i*30+13] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
580 LFR_BP2_F0[i*30+23] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
594 aux = sqrt(compressed_spec_mat[i*30]*compressed_spec_mat[i*30+24]);
595 cross_re = compressed_spec_mat[i*30+7] / aux;
596 cross_im = compressed_spec_mat[i*30+8] / aux;
597 LFR_BP2[i*30+13] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
598 LFR_BP2[i*30+23] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
581 599 #ifdef DEBUG_TCH
582 printf("LFR_BP2_F0[i*30+13] for cross15_re (%16.8e) : %.3u\n",cross_re, LFR_BP2_F0[i*30+13]);
583 printf("LFR_BP2_F0[i*30+23] for cross15_im (%16.8e) : %.3u\n",cross_im, LFR_BP2_F0[i*30+23]);
600 printf("LFR_BP2[i*30+13] for cross15_re (%16.8e) : %.3u\n",cross_re, LFR_BP2[i*30+13]);
601 printf("LFR_BP2[i*30+23] for cross15_im (%16.8e) : %.3u\n",cross_im, LFR_BP2[i*30+23]);
584 602 #endif
585 603 // S23
586 aux = sqrt(compressed_spectral_matrix_f0[i*30+9]*compressed_spectral_matrix_f0[i*30+16]);
587 cross_re = compressed_spectral_matrix_f0[i*30+10] / aux;
588 cross_im = compressed_spectral_matrix_f0[i*30+11] / aux;
589 LFR_BP2_F0[i*30+14] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
590 LFR_BP2_F0[i*30+24] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
604 aux = sqrt(compressed_spec_mat[i*30+9]*compressed_spec_mat[i*30+16]);
605 cross_re = compressed_spec_mat[i*30+10] / aux;
606 cross_im = compressed_spec_mat[i*30+11] / aux;
607 LFR_BP2[i*30+14] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
608 LFR_BP2[i*30+24] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
591 609 #ifdef DEBUG_TCH
592 printf("LFR_BP2_F0[i*30+14] for cross23_re (%16.8e) : %.3u\n",cross_re, LFR_BP2_F0[i*30+14]);
593 printf("LFR_BP2_F0[i*30+24] for cross23_im (%16.8e) : %.3u\n",cross_im, LFR_BP2_F0[i*30+24]);
610 printf("LFR_BP2[i*30+14] for cross23_re (%16.8e) : %.3u\n",cross_re, LFR_BP2[i*30+14]);
611 printf("LFR_BP2[i*30+24] for cross23_im (%16.8e) : %.3u\n",cross_im, LFR_BP2[i*30+24]);
594 612 #endif
595 613 // S24
596 aux = sqrt(compressed_spectral_matrix_f0[i*30+9]*compressed_spectral_matrix_f0[i*30+21]);
597 cross_re = compressed_spectral_matrix_f0[i*30+12] / aux;
598 cross_im = compressed_spectral_matrix_f0[i*30+13] / aux;
599 LFR_BP2_F0[i*30+15] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
600 LFR_BP2_F0[i*30+25] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
614 aux = sqrt(compressed_spec_mat[i*30+9]*compressed_spec_mat[i*30+21]);
615 cross_re = compressed_spec_mat[i*30+12] / aux;
616 cross_im = compressed_spec_mat[i*30+13] / aux;
617 LFR_BP2[i*30+15] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
618 LFR_BP2[i*30+25] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
601 619 #ifdef DEBUG_TCH
602 printf("LFR_BP2_F0[i*30+15] for cross24_re (%16.8e) : %.3u\n",cross_re, LFR_BP2_F0[i*30+15]);
603 printf("LFR_BP2_F0[i*30+25] for cross24_im (%16.8e) : %.3u\n",cross_im, LFR_BP2_F0[i*30+25]);
620 printf("LFR_BP2[i*30+15] for cross24_re (%16.8e) : %.3u\n",cross_re, LFR_BP2[i*30+15]);
621 printf("LFR_BP2[i*30+25] for cross24_im (%16.8e) : %.3u\n",cross_im, LFR_BP2[i*30+25]);
604 622 #endif
605 623 // S25
606 aux = sqrt(compressed_spectral_matrix_f0[i*30+9]*compressed_spectral_matrix_f0[i*30+24]);
607 cross_re = compressed_spectral_matrix_f0[i*30+14] / aux;
608 cross_im = compressed_spectral_matrix_f0[i*30+15] / aux;
609 LFR_BP2_F0[i*30+16] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
610 LFR_BP2_F0[i*30+26] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
624 aux = sqrt(compressed_spec_mat[i*30+9]*compressed_spec_mat[i*30+24]);
625 cross_re = compressed_spec_mat[i*30+14] / aux;
626 cross_im = compressed_spec_mat[i*30+15] / aux;
627 LFR_BP2[i*30+16] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
628 LFR_BP2[i*30+26] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
611 629 #ifdef DEBUG_TCH
612 printf("LFR_BP2_F0[i*30+16] for cross25_re (%16.8e) : %.3u\n",cross_re, LFR_BP2_F0[i*30+16]);
613 printf("LFR_BP2_F0[i*30+26] for cross25_im (%16.8e) : %.3u\n",cross_im, LFR_BP2_F0[i*30+26]);
630 printf("LFR_BP2[i*30+16] for cross25_re (%16.8e) : %.3u\n",cross_re, LFR_BP2[i*30+16]);
631 printf("LFR_BP2[i*30+26] for cross25_im (%16.8e) : %.3u\n",cross_im, LFR_BP2[i*30+26]);
614 632 #endif
615 633 // S34
616 aux = sqrt(compressed_spectral_matrix_f0[i*30+16]*compressed_spectral_matrix_f0[i*30+21]);
617 cross_re = compressed_spectral_matrix_f0[i*30+17] / aux;
618 cross_im = compressed_spectral_matrix_f0[i*30+18] / aux;
619 LFR_BP2_F0[i*30+17] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
620 LFR_BP2_F0[i*30+27] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
634 aux = sqrt(compressed_spec_mat[i*30+16]*compressed_spec_mat[i*30+21]);
635 cross_re = compressed_spec_mat[i*30+17] / aux;
636 cross_im = compressed_spec_mat[i*30+18] / aux;
637 LFR_BP2[i*30+17] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
638 LFR_BP2[i*30+27] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
621 639 #ifdef DEBUG_TCH
622 printf("LFR_BP2_F0[i*30+17] for cross34_re (%16.8e) : %.3u\n",cross_re, LFR_BP2_F0[i*30+17]);
623 printf("LFR_BP2_F0[i*30+27] for cross34_im (%16.8e) : %.3u\n",cross_im, LFR_BP2_F0[i*30+27]);
640 printf("LFR_BP2[i*30+17] for cross34_re (%16.8e) : %.3u\n",cross_re, LFR_BP2[i*30+17]);
641 printf("LFR_BP2[i*30+27] for cross34_im (%16.8e) : %.3u\n",cross_im, LFR_BP2[i*30+27]);
624 642 #endif
625 643 // S35
626 aux = sqrt(compressed_spectral_matrix_f0[i*30+16]*compressed_spectral_matrix_f0[i*30+24]);
627 cross_re = compressed_spectral_matrix_f0[i*30+19] / aux;
628 cross_im = compressed_spectral_matrix_f0[i*30+20] / aux;
629 LFR_BP2_F0[i*30+18] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
630 LFR_BP2_F0[i*30+28] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
644 aux = sqrt(compressed_spec_mat[i*30+16]*compressed_spec_mat[i*30+24]);
645 cross_re = compressed_spec_mat[i*30+19] / aux;
646 cross_im = compressed_spec_mat[i*30+20] / aux;
647 LFR_BP2[i*30+18] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
648 LFR_BP2[i*30+28] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
631 649 #ifdef DEBUG_TCH
632 printf("LFR_BP2_F0[i*30+18] for cross35_re (%16.8e) : %.3u\n",cross_re, LFR_BP2_F0[i*30+18]);
633 printf("LFR_BP2_F0[i*30+28] for cross35_im (%16.8e) : %.3u\n",cross_im, LFR_BP2_F0[i*30+28]);
650 printf("LFR_BP2[i*30+18] for cross35_re (%16.8e) : %.3u\n",cross_re, LFR_BP2[i*30+18]);
651 printf("LFR_BP2[i*30+28] for cross35_im (%16.8e) : %.3u\n",cross_im, LFR_BP2[i*30+28]);
634 652 #endif
635 653 // S45
636 aux = sqrt(compressed_spectral_matrix_f0[i*30+21]*compressed_spectral_matrix_f0[i*30+24]);
637 cross_re = compressed_spectral_matrix_f0[i*30+22] / aux;
638 cross_im = compressed_spectral_matrix_f0[i*30+23] / aux;
639 LFR_BP2_F0[i*30+19] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
640 LFR_BP2_F0[i*30+29] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
654 aux = sqrt(compressed_spec_mat[i*30+21]*compressed_spec_mat[i*30+24]);
655 cross_re = compressed_spec_mat[i*30+22] / aux;
656 cross_im = compressed_spec_mat[i*30+23] / aux;
657 LFR_BP2[i*30+19] = (unsigned char) (cross_re*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
658 LFR_BP2[i*30+29] = (unsigned char) (cross_im*127.5 + 128); // shift and cast into a 8-bit unsigned char (0, ..., 255) with rounding
641 659 #ifdef DEBUG_TCH
642 printf("LFR_BP2_F0[i*30+19] for cross45_re (%16.8e) : %.3u\n",cross_re, LFR_BP2_F0[i*30+19]);
643 printf("LFR_BP2_F0[i*30+29] for cross45_im (%16.8e) : %.3u\n",cross_im, LFR_BP2_F0[i*30+29]);
660 printf("LFR_BP2[i*30+19] for cross45_re (%16.8e) : %.3u\n",cross_re, LFR_BP2[i*30+19]);
661 printf("LFR_BP2[i*30+29] for cross45_im (%16.8e) : %.3u\n",cross_im, LFR_BP2[i*30+29]);
644 662 #endif
645 663 //==============================================
646 664 // BP2 auto correlations == PA_LFR_SC_BP2_AUTO_F0 == 5*16 bits = 5*[6 bits (exponent) + 10 bits (significand)]
647 665 // == PA_LFR_SC_BP2_AUTO_A0_F0 == 16 bits
648 666 // == PA_LFR_SC_BP2_AUTO_A1_F0 == 16 bits
649 667 // == PA_LFR_SC_BP2_AUTO_A2_F0 == 16 bits
650 668 // == PA_LFR_SC_BP2_AUTO_A3_F0 == 16 bits
651 669 // == PA_LFR_SC_BP2_AUTO_A4_F0 == 16 bits
652 670 // S11
653 significand = frexpf(compressed_spectral_matrix_f0[i*30], &exponent); // 0.5 <= significand < 1
671 significand = frexpf(compressed_spec_mat[i*30], &exponent); // 0.5 <= significand < 1
654 672 // S11 = significand * 2^exponent
655 673 #ifdef DEBUG_TCH
656 printf("S11 : %16.8e\n",compressed_spectral_matrix_f0[i*30]);
674 printf("S11 : %16.8e\n",compressed_spec_mat[i*30]);
657 675 printf("significand : %16.8e\n",significand);
658 676 printf("exponent : %d\n" ,exponent);
659 677 #endif
660 678 if (exponent < expmin) { // value should be >= 0.5 * 2^expmin
661 679 exponent = expmin;
662 680 significand = 0.5; // min value that can be recorded
663 681 }
664 682 if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)
665 683 exponent = expmax;
666 684 significand = 1.0; // max value that can be recorded
667 685 }
668 686 if (significand == 0) {// in that case exponent == 0 too
669 687 exponent = expmin;
670 688 significand = 0.5; // min value that can be recorded
671 689 }
672 690
673 691 autocor = (unsigned short int) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding
674 692 // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)
675 693 tmp_u_short_int = (unsigned short int) (exponent-expmin); // Shift and cast into a 16-bit unsigned int
676 694 // where just the first nbitexp bits are used (0, ..., 2^nbitexp-1)
677 pt_u_short_int = (unsigned short int*) &LFR_BP2_F0[i*30+0]; // Affect an unsigned short int pointer with the
695 pt_u_short_int = (unsigned short int*) &LFR_BP2[i*30+0]; // Affect an unsigned short int pointer with the
678 696 // adress where the 16-bit word result will be stored
679 697 *pt_u_short_int = autocor | (tmp_u_short_int << nbitsig); // Put the exponent bits (nbitexp) next to the
680 698 // left place of the significand bits (nbitsig), making
681 699 // the 16-bit word to be recorded, and record it using the pointer
682 700 #ifdef DEBUG_TCH
683 701 printf("autocor for S11 significand : %u\n",autocor );
684 702 printf("tmp_u_char for S11 exponent : %u\n",tmp_u_short_int );
685 703 printf("*pt_u_short_int for S11 exponent + significand : %.3d or %x\n",*pt_u_short_int, *pt_u_short_int);
686 printf("LFR_BP2_F0[i*30+1] : %u or %x\n",LFR_BP2_F0[i*30+1], LFR_BP2_F0[i*30+1]);
687 printf("LFR_BP2_F0[i*30+0] : %u or %x\n",LFR_BP2_F0[i*30+0], LFR_BP2_F0[i*30+0]);
704 printf("LFR_BP2[i*30+1] : %u or %x\n",LFR_BP2[i*30+1], LFR_BP2[i*30+1]);
705 printf("LFR_BP2[i*30+0] : %u or %x\n",LFR_BP2[i*30+0], LFR_BP2[i*30+0]);
688 706 #endif
689 707 // S22
690 significand = frexpf(compressed_spectral_matrix_f0[i*30+9], &exponent); // 0.5 <= significand < 1
708 significand = frexpf(compressed_spec_mat[i*30+9], &exponent); // 0.5 <= significand < 1
691 709 // S22 = significand * 2^exponent
692 710 #ifdef DEBUG_TCH
693 printf("S22 : %16.8e\n",compressed_spectral_matrix_f0[i*30+9]);
711 printf("S22 : %16.8e\n",compressed_spec_mat[i*30+9]);
694 712 printf("significand : %16.8e\n",significand);
695 713 printf("exponent : %d\n" ,exponent);
696 714 #endif
697 715 if (exponent < expmin) { // value should be >= 0.5 * 2^expmin
698 716 exponent = expmin;
699 717 significand = 0.5; // min value that can be recorded
700 718 }
701 719 if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)
702 720 exponent = expmax;
703 721 significand = 1.0; // max value that can be recorded
704 722 }
705 723 if (significand == 0) {// in that case exponent == 0 too
706 724 exponent = expmin;
707 725 significand = 0.5; // min value that can be recorded
708 726 }
709 727
710 728 autocor = (unsigned short int) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding
711 729 // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)
712 730 tmp_u_short_int = (unsigned short int) (exponent-expmin); // Shift and cast into a 16-bit unsigned int
713 731 // where just the first nbitexp bits are used (0, ..., 2^nbitexp-1)
714 pt_u_short_int = (unsigned short int*) &LFR_BP2_F0[i*30+2]; // Affect an unsigned short int pointer with the
732 pt_u_short_int = (unsigned short int*) &LFR_BP2[i*30+2]; // Affect an unsigned short int pointer with the
715 733 // adress where the 16-bit word result will be stored
716 734 *pt_u_short_int = autocor | (tmp_u_short_int << nbitsig); // Put the exponent bits (nbitexp) next to the
717 735 // left place of the significand bits (nbitsig), making
718 736 // the 16-bit word to be recorded, and record it using the pointer
719 737 #ifdef DEBUG_TCH
720 738 printf("autocor for S22 significand : %d\n",autocor );
721 739 printf("tmp_u_char for S22 exponent : %d\n",tmp_u_short_int );
722 740 printf("*pt_u_short_int for S22 exponent + significand : %.3d or %x\n",*pt_u_short_int, *pt_u_short_int);
723 printf("LFR_BP2_F0[i*30+3] : %.3d or %x\n",LFR_BP2_F0[i*30+3], LFR_BP2_F0[i*30+3]);
724 printf("LFR_BP2_F0[i*30+2] : %.3d or %x\n",LFR_BP2_F0[i*30+2], LFR_BP2_F0[i*30+2]);
741 printf("LFR_BP2[i*30+3] : %.3d or %x\n",LFR_BP2[i*30+3], LFR_BP2[i*30+3]);
742 printf("LFR_BP2[i*30+2] : %.3d or %x\n",LFR_BP2[i*30+2], LFR_BP2[i*30+2]);
725 743 #endif
726 744 // S33
727 significand = frexpf(compressed_spectral_matrix_f0[i*30+16], &exponent); // 0.5 <= significand < 1
745 significand = frexpf(compressed_spec_mat[i*30+16], &exponent); // 0.5 <= significand < 1
728 746 // S33 = significand * 2^exponent
729 747 #ifdef DEBUG_TCH
730 printf("S33 : %16.8e\n",compressed_spectral_matrix_f0[i*30+16]);
748 printf("S33 : %16.8e\n",compressed_spec_mat[i*30+16]);
731 749 printf("significand : %16.8e\n",significand);
732 750 printf("exponent : %d\n" ,exponent);
733 751 #endif
734 752 if (exponent < expmin) { // value should be >= 0.5 * 2^expmin
735 753 exponent = expmin;
736 754 significand = 0.5; // min value that can be recorded
737 755 }
738 756 if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)
739 757 exponent = expmax;
740 758 significand = 1.0; // max value that can be recorded
741 759 }
742 760 if (significand == 0) {// in that case exponent == 0 too
743 761 exponent = expmin;
744 762 significand = 0.5; // min value that can be recorded
745 763 }
746 764
747 765 autocor = (unsigned short int) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding
748 766 // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)
749 767 tmp_u_short_int = (unsigned short int) (exponent-expmin); // Shift and cast into a 16-bit unsigned int
750 768 // where just the first nbitexp bits are used (0, ..., 2^nbitexp-1)
751 pt_u_short_int = (unsigned short int*) &LFR_BP2_F0[i*30+4]; // Affect an unsigned short int pointer with the
769 pt_u_short_int = (unsigned short int*) &LFR_BP2[i*30+4]; // Affect an unsigned short int pointer with the
752 770 // adress where the 16-bit word result will be stored
753 771 *pt_u_short_int = autocor | (tmp_u_short_int << nbitsig); // Put the exponent bits (nbitexp) next to the
754 772 // left place of the significand bits (nbitsig), making
755 773 // the 16-bit word to be recorded, and record it using the pointer
756 774 #ifdef DEBUG_TCH
757 775 printf("autocor for S33 significand : %d\n",autocor );
758 776 printf("tmp_u_char for S33 exponent : %d\n",tmp_u_short_int );
759 777 printf("*pt_u_short_int for S33 exponent + significand : %.3d or %x\n",*pt_u_short_int, *pt_u_short_int);
760 printf("LFR_BP2_F0[i*30+5] : %.3d or %x\n",LFR_BP2_F0[i*30+5], LFR_BP2_F0[i*30+5]);
761 printf("LFR_BP2_F0[i*30+4] : %.3d or %x\n",LFR_BP2_F0[i*30+4], LFR_BP2_F0[i*30+4]);
778 printf("LFR_BP2[i*30+5] : %.3d or %x\n",LFR_BP2[i*30+5], LFR_BP2[i*30+5]);
779 printf("LFR_BP2[i*30+4] : %.3d or %x\n",LFR_BP2[i*30+4], LFR_BP2[i*30+4]);
762 780 #endif
763 781 // S44
764 significand = frexpf(compressed_spectral_matrix_f0[i*30+21], &exponent); // 0.5 <= significand < 1
782 significand = frexpf(compressed_spec_mat[i*30+21], &exponent); // 0.5 <= significand < 1
765 783 // S44 = significand * 2^exponent
766 784 #ifdef DEBUG_TCH
767 printf("S44 : %16.8e\n",compressed_spectral_matrix_f0[i*30+21]);
785 printf("S44 : %16.8e\n",compressed_spec_mat[i*30+21]);
768 786 printf("significand : %16.8e\n",significand);
769 787 printf("exponent : %d\n" ,exponent);
770 788 #endif
771 789
772 790 if (exponent < expmin) { // value should be >= 0.5 * 2^expmin
773 791 exponent = expmin;
774 792 significand = 0.5; // min value that can be recorded
775 793 }
776 794 if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)
777 795 exponent = expmax;
778 796 significand = 1.0; // max value that can be recorded
779 797 }
780 798 if (significand == 0) {// in that case exponent == 0 too
781 799 exponent = expmin;
782 800 significand = 0.5; // min value that can be recorded
783 801 }
784 802
785 803 autocor = (unsigned short int) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding
786 804 // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)
787 805 tmp_u_short_int = (unsigned short int) (exponent-expmin); // Shift and cast into a 16-bit unsigned int
788 806 // where just the first nbitexp bits are used (0, ..., 2^nbitexp-1)
789 pt_u_short_int = (unsigned short int*) &LFR_BP2_F0[i*30+6]; // Affect an unsigned short int pointer with the
807 pt_u_short_int = (unsigned short int*) &LFR_BP2[i*30+6]; // Affect an unsigned short int pointer with the
790 808 // adress where the 16-bit word result will be stored
791 809 *pt_u_short_int = autocor | (tmp_u_short_int << nbitsig); // Put the exponent bits (nbitexp) next to the
792 810 // left place of the significand bits (nbitsig), making
793 811 // the 16-bit word to be recorded, and record it using the pointer
794 812 #ifdef DEBUG_TCH
795 813 printf("autocor for S44 significand : %d\n",autocor );
796 814 printf("tmp_u_char for S44 exponent : %d\n",tmp_u_short_int );
797 815 printf("*pt_u_short_int for S44 exponent + significand : %.3d or %x\n",*pt_u_short_int, *pt_u_short_int);
798 printf("LFR_BP2_F0[i*30+7] : %.3d or %x\n",LFR_BP2_F0[i*30+7], LFR_BP2_F0[i*30+7]);
799 printf("LFR_BP2_F0[i*30+6] : %.3d or %x\n",LFR_BP2_F0[i*30+6], LFR_BP2_F0[i*30+6]);
816 printf("LFR_BP2[i*30+7] : %.3d or %x\n",LFR_BP2[i*30+7], LFR_BP2[i*30+7]);
817 printf("LFR_BP2[i*30+6] : %.3d or %x\n",LFR_BP2[i*30+6], LFR_BP2[i*30+6]);
800 818 #endif
801 819 // S55
802 significand = frexpf(compressed_spectral_matrix_f0[i*30+24], &exponent); // 0.5 <= significand < 1
820 significand = frexpf(compressed_spec_mat[i*30+24], &exponent); // 0.5 <= significand < 1
803 821 // S55 = significand * 2^exponent
804 822 #ifdef DEBUG_TCH
805 printf("S55 : %16.8e\n",compressed_spectral_matrix_f0[i*30+24]);
823 printf("S55 : %16.8e\n",compressed_spec_mat[i*30+24]);
806 824 printf("significand : %16.8e\n",significand);
807 825 printf("exponent : %d\n" ,exponent);
808 826 #endif
809 827 if (exponent < expmin) { // value should be >= 0.5 * 2^expmin
810 828 exponent = expmin;
811 829 significand = 0.5; // min value that can be recorded
812 830 }
813 831 if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)
814 832 exponent = expmax;
815 833 significand = 1.0; // max value that can be recorded
816 834 }
817 835 if (significand == 0) {// in that case exponent == 0 too
818 836 exponent = expmin;
819 837 significand = 0.5; // min value that can be recorded
820 838 }
821 839
822 840 autocor = (unsigned short int) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding
823 841 // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)
824 842 tmp_u_short_int = (unsigned short int) (exponent-expmin); // Shift and cast into a 16-bit unsigned int
825 843 // where just the first nbitexp bits are used (0, ..., 2^nbitexp-1)
826 pt_u_short_int = (unsigned short int*) &LFR_BP2_F0[i*30+8]; // Affect an unsigned short int pointer with the
844 pt_u_short_int = (unsigned short int*) &LFR_BP2[i*30+8]; // Affect an unsigned short int pointer with the
827 845 // adress where the 16-bit word result will be stored
828 846 *pt_u_short_int = autocor | (tmp_u_short_int << nbitsig); // Put the exponent bits (nbitexp) next to the
829 847 // left place of the significand bits (nbitsig), making
830 848 // the 16-bit word to be recorded, and record it using the pointer
831 849 #ifdef DEBUG_TCH
832 850 printf("autocor for S55 significand : %d\n",autocor );
833 851 printf("tmp_u_char for S55 exponent : %d\n",tmp_u_short_int );
834 852 printf("*pt_u_short_int for S55 exponent + significand : %.3d or %x\n",*pt_u_short_int, *pt_u_short_int);
835 printf("LFR_BP2_F0[i*30+9] : %.3d or %x\n",LFR_BP2_F0[i*30+9], LFR_BP2_F0[i*30+9]);
836 printf("LFR_BP2_F0[i*30+8] : %.3d or %x\n",LFR_BP2_F0[i*30+8], LFR_BP2_F0[i*30+8]);
853 printf("LFR_BP2[i*30+9] : %.3d or %x\n",LFR_BP2[i*30+9], LFR_BP2[i*30+9]);
854 printf("LFR_BP2[i*30+8] : %.3d or %x\n",LFR_BP2[i*30+8], LFR_BP2[i*30+8]);
837 855 #endif
838 856 }
839 857 }
840 858
Show file before | Show file after | Hide comments
@@ -1,32 +1,13
1 1 // In the frame of RPW LFR Sofware ICD Issue1 Rev8 (05/07/2013)
2 2 // version 1: 31/07/2013
3 3
4 4 #ifndef BASIC_PARAMETERS_H_INCLUDED
5 5 #define BASIC_PARAMETERS_H_INCLUDED
6 6
7 #define LPP_SPECTRAL_MATRIX_CTRL 0x80000700
8 #define LPP_SPECTRAL_MATRIX_1 0x80000704
9 #define LPP_SPECTRAL_MATRIX_2 0x80000708
7 #define NB_BINS_COMPRESSED_MATRIX_f0 1
10 8
11 #define NB_BINS_SPECTRAL_MATRIX 128
12 #define NB_VALUES_PER_SPECTRAL_MATRIX 25
13 #define TOTAL_SIZE_SPECTRAL_MATRIX NB_BINS_SPECTRAL_MATRIX * NB_VALUES_PER_SPECTRAL_MATRIX
14 #define NB_BINS_COMPRESSED_MATRIX_f0 1
15 #define SIZE_COMPRESSED_SPECTRAL_MATRIX_f1 13
16 #define SIZE_COMPRESSED_SPECTRAL_MATRIX_f2 12
17 #define TOTAL_SIZE_COMPRESSED_MATRIX_f0 NB_BINS_COMPRESSED_MATRIX_f0 * NB_VALUES_PER_SPECTRAL_MATRIX
18 #define NB_AVERAGE_NORMAL_f0 4
19
20 volatile int spectral_matrix_f0_a[TOTAL_SIZE_SPECTRAL_MATRIX];
21 volatile int spectral_matrix_f0_b[TOTAL_SIZE_SPECTRAL_MATRIX];
22 int averaged_spectral_matrix_f0[TOTAL_SIZE_SPECTRAL_MATRIX];
23
24 extern float compressed_spectral_matrix_f0[];
25
26 extern unsigned char LFR_BP1_F0[];
27 extern unsigned char LFR_BP2_F0[];
28
29 void BP1_set();
30 void BP2_set();
9 void init_k_f0( void );
10 void BP1_set(float * compressed_spec_mat, unsigned char nb_bins_compressed_spec_mat, unsigned char * LFR_BP1);
11 void BP2_set(float * compressed_spec_mat, unsigned char nb_bins_compressed_spec_mat, unsigned char * LFR_BP2);
31 12
32 13 #endif // BASIC_PARAMETERS_H_INCLUDED
General Comments 0
You need to be logged in to leave comments. Login now