# HG changeset patch # User chust # Date 2014-05-16 07:49:56 # Node ID 163519b5356a8bcd9197f504c31ad34b6233de38 # Parent f66dcaa5b47ba879af0690acdbce09f970609c41 version 1.4 déclaration des fonctions static inline (BP2 + BP1 ok cf V1.3) diff --git a/basic_parameters.c b/basic_parameters.c --- a/basic_parameters.c +++ b/basic_parameters.c @@ -5,47 +5,6 @@ // version 1.3: 02/05/2014 #include "basic_parameters.h" -#include -#include -#include - -#define K44_PE 0 -#define K55_PE 1 -#define K45_PE_RE 2 -#define K45_PE_IM 3 - -#define K14_SX_RE 4 -#define K14_SX_IM 5 -#define K15_SX_RE 6 -#define K15_SX_IM 7 -#define K24_SX_RE 8 -#define K24_SX_IM 9 -#define K25_SX_RE 10 -#define K25_SX_IM 11 -#define K34_SX_RE 12 -#define K34_SX_IM 13 -#define K35_SX_RE 14 -#define K35_SX_IM 15 - -#define K24_NY_RE 16 -#define K24_NY_IM 17 -#define K25_NY_RE 18 -#define K25_NY_IM 19 -#define K34_NY_RE 20 -#define K34_NY_IM 21 -#define K35_NY_RE 22 -#define K35_NY_IM 23 - -#define K24_NZ_RE 24 -#define K24_NZ_IM 25 -#define K25_NZ_RE 26 -#define K25_NZ_IM 27 -#define K34_NZ_RE 28 -#define K34_NZ_IM 29 -#define K35_NZ_RE 30 -#define K35_NZ_IM 31 - -float k_f0[NB_BINS_COMPRESSED_MATRIX_f0][32]; void init_k_f0( void ) { @@ -90,842 +49,3 @@ void init_k_f0( void ) } } -float alpha_M = 45 * (3.1415927/180); - -void BP1_set( float * compressed_spec_mat, uint8_t nb_bins_compressed_spec_mat, uint8_t * lfr_bp1 ){ - float PSDB; // 32-bit floating point - float PSDE; - float tmp; - float NVEC_V0; - float NVEC_V1; - float NVEC_V2; - float aux; - float tr_SB_SB; - float e_cross_b_re; - float e_cross_b_im; - float n_cross_e_scal_b_re; - float n_cross_e_scal_b_im; - float ny; - float nz; - float bx_bx_star; - float vphi; - float significand; - int exponent; // 32-bit signed integer - uint8_t nbitexp; // 8-bit unsigned integer - uint8_t nbitsig; - uint8_t tmp_uint8; - uint8_t *pt_uint8; // pointer on unsigned 8-bit integer - int8_t expmin; // 8-bit signed integer - int8_t expmax; - uint16_t rangesig; // 16-bit unsigned integer - uint16_t psd; - uint16_t exp; - uint16_t tmp_uint16; - uint16_t i; - - init_k_f0(); - -#ifdef DEBUG_TCH - printf("BP1 : \n"); - printf("Number of bins: %d\n", nb_bins_compressed_spec_mat); -#endif - - // initialization for managing the exponents of the floating point data: - nbitexp = 5; // number of bits for the exponent - expmax = 30; // maximum value of the exponent - expmin = expmax - (1 << nbitexp) + 1; // accordingly the minimum exponent value - // for floating point data to be recorded on 12-bit words: - nbitsig = 12 - nbitexp; // number of bits for the significand - rangesig = (1 << nbitsig)-1; // == 2^nbitsig - 1 - -#ifdef DEBUG_TCH - printf("nbitexp : %d, expmax : %d, expmin : %d\n", nbitexp, expmax, expmin); - printf("nbitsig : %d, rangesig : %d\n", nbitsig, rangesig); -#endif - - for(i=0; i= 0.5 * 2^expmin - exponent = expmin; - significand = 0.5; // min value that can be recorded - } - if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1) - exponent = expmax; - significand = 1.0; // max value that can be recorded - } - if (significand == 0) { // in that case exponent == 0 too - exponent = expmin; - significand = 0.5; // min value that can be recorded - } - - psd = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding - // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1) - exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just - // the first nbitexp bits are used (0, ..., 2^nbitexp-1) - tmp_uint16 = psd | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the - // left place of the significand bits (nbitsig), - // making the 16-bit word to be recorded - pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16 -#ifdef LSB_FIRST_TCH - lfr_bp1[i*NB_BYTES_BP1+2] = pt_uint8[0]; // Record LSB of tmp_uint16 - lfr_bp1[i*NB_BYTES_BP1+3] = pt_uint8[1]; // Record MSB of tmp_uint16 -#endif -#ifdef MSB_FIRST_TCH - lfr_bp1[i*NB_BYTES_BP1+2] = pt_uint8[1]; // Record LSB of tmp_uint16 - lfr_bp1[i*NB_BYTES_BP1+3] = pt_uint8[0]; // Record MSB of tmp_uint16 -#endif -#ifdef DEBUG_TCH - printf("\nBin number: %d\n", i); - printf("PSDB / 3 : %16.8e\n",PSDB/3); - printf("significand : %16.8e\n",significand); - printf("exponent : %d\n" ,exponent); - printf("psd for PSDB significand : %d\n",psd); - printf("exp for PSDB exponent : %d\n",exp); - printf("pt_uint8[1] for PSDB exponent + significand: %.3d or %.2x\n",pt_uint8[1], pt_uint8[1]); - printf("pt_uint8[0] for PSDB exponent + significand: %.3d or %.2x\n",pt_uint8[0], pt_uint8[0]); - printf("lfr_bp1[i*NB_BYTES_BP1+3] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+3], lfr_bp1[i*NB_BYTES_BP1+3]); - printf("lfr_bp1[i*NB_BYTES_BP1+2] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+2], lfr_bp1[i*NB_BYTES_BP1+2]); -#endif - //============================================== - // BP1 PSDE == PA_LFR_SC_BP1_PE_F0 == 12 bits = 5 bits (exponent) + 7 bits (significand) - PSDE = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21] * k_f0[i][K44_PE] // S44 - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24] * k_f0[i][K55_PE] // S55 - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+22] * k_f0[i][K45_PE_RE] // S45 Re - - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+23] * k_f0[i][K45_PE_IM]; // S45 Im - - significand = frexpf(PSDE/2, &exponent); // 0.5 <= significand < 1 - // PSDE/2 = significand * 2^exponent - // the division by 2 is to ensure that max value <= 2^30 - // should be reconsidered by taking into account the k-coefficients ... - - if (exponent < expmin) { // value should be >= 0.5 * 2^expmin - exponent = expmin; - significand = 0.5; // min value that can be recorded - } - if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1) - exponent = expmax; - significand = 1.0; // max value that can be recorded - } - if (significand == 0) {// in that case exponent == 0 too - exponent = expmin; - significand = 0.5; // min value that can be recorded - } - - psd = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding - // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1) - exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just - // the first nbitexp bits are used (0, ..., 2^nbitexp-1) - tmp_uint16 = psd | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the - // left place of the significand bits (nbitsig), - // making the 16-bit word to be recorded - pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16 -#ifdef LSB_FIRST_TCH - lfr_bp1[i*NB_BYTES_BP1+0] = pt_uint8[0]; // Record LSB of tmp_uint16 - lfr_bp1[i*NB_BYTES_BP1+1] = pt_uint8[1]; // Record MSB of tmp_uint16 -#endif -#ifdef MSB_FIRST_TCH - lfr_bp1[i*NB_BYTES_BP1+0] = pt_uint8[1]; // Record LSB of tmp_uint16 - lfr_bp1[i*NB_BYTES_BP1+1] = pt_uint8[0]; // Record MSB of tmp_uint16 -#endif -#ifdef DEBUG_TCH - printf("Bin number: %d\n", i); - printf("PSDE/2 : %16.8e\n",PSDE/2); - printf("significand : %16.8e\n",significand); - printf("exponent : %d\n" ,exponent); - printf("psd for PSDE significand : %d\n",psd); - printf("exp for PSDE exponent : %d\n",exp); - printf("pt_uint8[1] for PSDE exponent + significand: %.3d or %.2x\n",pt_uint8[1], pt_uint8[1]); - printf("pt_uint8[0] for PSDE exponent + significand: %.3d or %.2x\n",pt_uint8[0], pt_uint8[0]); - printf("lfr_bp1[i*NB_BYTES_BP1+1] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+1], lfr_bp1[i*NB_BYTES_BP1+1]); - printf("lfr_bp1[i*NB_BYTES_BP1+0] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+0], lfr_bp1[i*NB_BYTES_BP1+0]); -#endif - //============================================================================== - // BP1 normal wave vector == PA_LFR_SC_BP1_NVEC_V0_F0 == 8 bits - // == PA_LFR_SC_BP1_NVEC_V1_F0 == 8 bits - // == PA_LFR_SC_BP1_NVEC_V2_F0 == 1 sign bit - tmp = sqrt( compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] //Im S12 - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] //Im S13 - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11] //Im S23 - ); - NVEC_V0 = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11]/ tmp; // S23 Im => n1 - NVEC_V1 = -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] / tmp; // S13 Im => n2 - NVEC_V2 = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] / tmp; // S12 Im => n3 - - lfr_bp1[i*NB_BYTES_BP1+4] = (uint8_t) (NVEC_V0*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding - lfr_bp1[i*NB_BYTES_BP1+5] = (uint8_t) (NVEC_V1*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding - pt_uint8 = (uint8_t*) &NVEC_V2; // Affect an uint8_t pointer with the adress of NVEC_V2 -#ifdef LSB_FIRST_TCH - lfr_bp1[i*NB_BYTES_BP1+6] = pt_uint8[3] & 0x80; // Extract the sign bit of NVEC_V2 (32-bit float, sign bit in the 4th octet:PC convention) - // Record it at the 8th bit position (from the right to the left) of lfr_bp1[i*NB_BYTES_BP1+6] -#endif -#ifdef MSB_FIRST_TCH - lfr_bp1[i*NB_BYTES_BP1+6] = pt_uint8[0] & 0x80; // Extract the sign bit of NVEC_V2 (32-bit float, sign bit in the 0th octet:SPARC convention) - // Record it at the 8th bit position (from the right to the left) of lfr_bp1[i*NB_BYTES_BP1+6] -#endif -#ifdef DEBUG_TCH - printf("NVEC_V0 : %16.8e\n",NVEC_V0); - printf("NVEC_V1 : %16.8e\n",NVEC_V1); - printf("NVEC_V2 : %16.8e\n",NVEC_V2); - printf("lfr_bp1[i*NB_BYTES_BP1+4] for NVEC_V0 : %u\n",lfr_bp1[i*NB_BYTES_BP1+4]); - printf("lfr_bp1[i*NB_BYTES_BP1+5] for NVEC_V1 : %u\n",lfr_bp1[i*NB_BYTES_BP1+5]); - printf("lfr_bp1[i*NB_BYTES_BP1+6] for NVEC_V2 : %u\n",lfr_bp1[i*NB_BYTES_BP1+6]); -#endif - //======================================================= - // BP1 ellipticity == PA_LFR_SC_BP1_ELLIP_F0 == 4 bits - aux = 2*tmp / PSDB; // Compute the ellipticity - - tmp_uint8 = (uint8_t) (aux*15 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding - // where just the first 4 bits are used (0, ..., 15) - lfr_bp1[i*NB_BYTES_BP1+6] = lfr_bp1[i*NB_BYTES_BP1+6] | (tmp_uint8 << 3); // Put these 4 bits next to the right place - // of the sign bit of NVEC_V2 (recorded - // previously in lfr_bp1[i*NB_BYTES_BP1+6]) -#ifdef DEBUG_TCH - printf("ellipticity : %16.8e\n",aux); - printf("tmp_uint8 for ellipticity : %u\n",tmp_uint8); - printf("lfr_bp1[i*NB_BYTES_BP1+6] for NVEC_V2 + ellipticity : %u\n",lfr_bp1[i*NB_BYTES_BP1+6]); -#endif - //============================================================== - // BP1 degree of polarization == PA_LFR_SC_BP1_DOP_F0 == 3 bits - tr_SB_SB = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX] - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9] - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16] - + 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+1] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+1] - + 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] - + 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+3] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+3] - + 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] - + 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10] - + 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11]; - aux = PSDB*PSDB; - tmp = ( 3*tr_SB_SB - aux ) / ( 2 * aux ); // Compute the degree of polarisation - - tmp_uint8 = (uint8_t) (tmp*7 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding - // where just the first 3 bits are used (0, ..., 7) - lfr_bp1[i*NB_BYTES_BP1+6] = lfr_bp1[i*NB_BYTES_BP1+6] | tmp_uint8; // Record these 3 bits at the 3 first bit positions - // (from the right to the left) of lfr_bp1[i*NB_BYTES_BP1+6] -#ifdef DEBUG_TCH - printf("DOP : %16.8e\n",tmp); - printf("tmp_uint8 for DOP : %u\n",tmp_uint8); - printf("lfr_bp1[i*NB_BYTES_BP1+6] for NVEC_V2 + ellipticity + DOP : %u\n",lfr_bp1[i*NB_BYTES_BP1+6]); -#endif - //======================================================================================= - // BP1 X_SO-component of the Poynting flux == PA_LFR_SC_BP1_SX_F0 == 8 (+ 2) bits - // = 5 bits (exponent) + 3 bits (significand) - // + 1 sign bit + 1 argument bit (two sectors) - e_cross_b_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_f0[i][K34_SX_RE] //S34 Re - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_f0[i][K35_SX_RE] //S35 Re - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+5] *k_f0[i][K14_SX_RE] //S14 Re - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+7] *k_f0[i][K15_SX_RE] //S15 Re - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_f0[i][K24_SX_RE] //S24 Re - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_f0[i][K25_SX_RE] //S25 Re - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_f0[i][K34_SX_IM] //S34 Im - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_f0[i][K35_SX_IM] //S35 Im - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+6] *k_f0[i][K14_SX_IM] //S14 Im - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+8] *k_f0[i][K15_SX_IM] //S15 Im - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_f0[i][K24_SX_IM] //S24 Im - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_f0[i][K25_SX_IM]; //S25 Im - // Im(S_ji) = -Im(S_ij) - // k_ji = k_ij - e_cross_b_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_f0[i][K34_SX_IM] //S34 Re - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_f0[i][K35_SX_IM] //S35 Re - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+5] *k_f0[i][K14_SX_IM] //S14 Re - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+7] *k_f0[i][K15_SX_IM] //S15 Re - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_f0[i][K24_SX_IM] //S24 Re - + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_f0[i][K25_SX_IM] //S25 Re - - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_f0[i][K34_SX_RE] //S34 Im - - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_f0[i][K35_SX_RE] //S35 Im - - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+6] *k_f0[i][K14_SX_RE] //S14 Im - - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+8] *k_f0[i][K15_SX_RE] //S15 Im - - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_f0[i][K24_SX_RE] //S24 Im - - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_f0[i][K25_SX_RE]; //S25 Im -#ifdef DEBUG_TCH - printf("ReaSX / 2 : %16.8e\n",e_cross_b_re/2); -#endif - pt_uint8 = (uint8_t*) &e_cross_b_re; // Affect an uint8_t pointer with the adress of e_cross_b_re -#ifdef LSB_FIRST_TCH - lfr_bp1[i*NB_BYTES_BP1+1] = lfr_bp1[i*NB_BYTES_BP1+1] | (pt_uint8[3] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 4th octet:PC convention) - // Record it at the 8th bit position (from the right to the left) - // of lfr_bp1[i*NB_BYTES_BP1+1] - pt_uint8[3] = (pt_uint8[3] & 0x7f); // Make e_cross_b_re be positive in any case: |ReaSX| -#endif -#ifdef MSB_FIRST_TCH - lfr_bp1[i*NB_BYTES_BP1+1] = lfr_bp1[i*NB_BYTES_BP1+1] | (pt_uint8[0] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 0th octet:SPARC convention) - // Record it at the 8th bit position (from the right to the left) - // of lfr_bp1[i*NB_BYTES_BP1+1] - pt_uint8[0] = (pt_uint8[0] & 0x7f); // Make e_cross_b_re be positive in any case: |ReaSX| -#endif - significand = frexpf(e_cross_b_re/2, &exponent); // 0.5 <= significand < 1 - // ReaSX/2 = significand * 2^exponent - // The division by 2 is to ensure that max value <= 2^30 (rough estimate) - // Should be reconsidered by taking into account the k-coefficients ... - - if (exponent < expmin) { // value should be >= 0.5 * 2^expmin - exponent = expmin; - significand = 0.5; // min value that can be recorded - } - if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1) - exponent = expmax; - significand = 1.0; // max value that can be recorded - } - if (significand == 0) { // in that case exponent == 0 too - exponent = expmin; - significand = 0.5; // min value that can be recorded - } - - lfr_bp1[i*NB_BYTES_BP1+7] = (uint8_t) ((significand*2-1)*7 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding - // where just the first 3 bits are used (0, ..., 7) - tmp_uint8 = (uint8_t) (exponent-expmin); // Shift and cast into a 8-bit uint8_t where - // just the first 5 bits are used (0, ..., 2^5-1) -#ifdef DEBUG_TCH - printf("|ReaSX| / 2 : %16.8e\n",e_cross_b_re/2); - printf("significand : %16.8e\n",significand); - printf("exponent : %d\n" ,exponent); - printf("lfr_bp1[i*NB_BYTES_BP1+7] for ReaSX significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+7]); - printf("tmp_uint8 for ReaSX exponent : %d\n",tmp_uint8); -#endif - lfr_bp1[i*NB_BYTES_BP1+7] = lfr_bp1[i*NB_BYTES_BP1+7] | (tmp_uint8 << 3); // Shift these 5 bits to the left before logical addition - // with lfr_bp1[i*NB_BYTES_BP1+7] -#ifdef DEBUG_TCH - printf("lfr_bp1[i*NB_BYTES_BP1+7] for ReaSX exponent + significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+7]); - printf("lfr_bp1[i*NB_BYTES_BP1+1] for ReaSX sign + PSDE 'exponent' : %u\n",lfr_bp1[i*NB_BYTES_BP1+1]); - printf("ImaSX / 2 : %16.8e\n",e_cross_b_im/2); -#endif - pt_uint8 = (uint8_t*) &e_cross_b_im; // Affect an uint8_t pointer with the adress of e_cross_b_im -#ifdef LSB_FIRST_TCH - pt_uint8[3] = pt_uint8[3] & 0x7f; // Make e_cross_b_im be positive in any case: |ImaSX| -#endif -#ifdef MSB_FIRST_TCH - pt_uint8[0] = pt_uint8[0] & 0x7f; // Make e_cross_b_im be positive in any case: |ImaSX| -#endif - tmp_uint8 = (e_cross_b_im > e_cross_b_re) ? 0x40 : 0x00; // Determine the sector argument of SX. If |Im| > |Re| affect - // an unsigned 8-bit char with 01000000; otherwise with null. - lfr_bp1[i*NB_BYTES_BP1+1] = lfr_bp1[i*NB_BYTES_BP1+1] | tmp_uint8; // Record it as a sign bit at the 7th bit position (from the right - // to the left) of lfr_bp1[i*NB_BYTES_BP1+1], by simple logical addition. -#ifdef DEBUG_TCH - printf("|ImaSX| / 2 : %16.8e\n",e_cross_b_im/2); - printf("ArgSX sign : %u\n",tmp_uint8); - printf("lfr_bp1[i*NB_BYTES_BP1+1] for ReaSX & ArgSX signs + PSDE 'exponent' : %u\n",lfr_bp1[i*NB_BYTES_BP1+1]); -#endif - //====================================================================== - // BP1 phase velocity estimator == PA_LFR_SC_BP1_VPHI_F0 == 8 (+ 2) bits - // = 5 bits (exponent) + 3 bits (significand) - // + 1 sign bit + 1 argument bit (two sectors) - ny = sin(alpha_M)*NVEC_V1 + cos(alpha_M)*NVEC_V2; - nz = NVEC_V0; - bx_bx_star = cos(alpha_M)*cos(alpha_M)*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9] // S22 Re - + sin(alpha_M)*sin(alpha_M)*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16] // S33 Re - - 2*sin(alpha_M)*cos(alpha_M)*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10]; // S23 Re - - n_cross_e_scal_b_re = ny * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_f0[i][K24_NY_RE] //S24 Re - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_f0[i][K25_NY_RE] //S25 Re - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_f0[i][K34_NY_RE] //S34 Re - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_f0[i][K35_NY_RE] //S35 Re - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_f0[i][K24_NY_IM] //S24 Im - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_f0[i][K25_NY_IM] //S25 Im - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_f0[i][K34_NY_IM] //S34 Im - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_f0[i][K35_NY_IM]) //S35 Im - + nz * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_f0[i][K24_NZ_RE] //S24 Re - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_f0[i][K25_NZ_RE] //S25 Re - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_f0[i][K34_NZ_RE] //S34 Re - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_f0[i][K35_NZ_RE] //S35 Re - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_f0[i][K24_NZ_IM] //S24 Im - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_f0[i][K25_NZ_IM] //S25 Im - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_f0[i][K34_NZ_IM] //S34 Im - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_f0[i][K35_NZ_IM]);//S35 Im - // Im(S_ji) = -Im(S_ij) - // k_ji = k_ij - n_cross_e_scal_b_im = ny * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_f0[i][K24_NY_IM] //S24 Re - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_f0[i][K25_NY_IM] //S25 Re - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_f0[i][K34_NY_IM] //S34 Re - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_f0[i][K35_NY_IM] //S35 Re - -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_f0[i][K24_NY_RE] //S24 Im - -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_f0[i][K25_NY_RE] //S25 Im - -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_f0[i][K34_NY_RE] //S34 Im - -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_f0[i][K35_NY_RE]) //S35 Im - + nz * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_f0[i][K24_NZ_IM] //S24 Re - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_f0[i][K25_NZ_IM] //S25 Re - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_f0[i][K34_NZ_IM] //S34 Re - +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_f0[i][K35_NZ_IM] //S35 Re - -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_f0[i][K24_NZ_RE] //S24 Im - -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_f0[i][K25_NZ_RE] //S25 Im - -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_f0[i][K34_NZ_RE] //S34 Im - -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_f0[i][K35_NZ_RE]);//S35 Im -#ifdef DEBUG_TCH - printf("n_cross_e_scal_b_re : %16.8e\n",n_cross_e_scal_b_re); - printf("n_cross_e_scal_b_im : %16.8e\n",n_cross_e_scal_b_im); -#endif - // vphi = n_cross_e_scal_b_re / bx_bx_star => sign(VPHI) = sign(n_cross_e_scal_b_re) - pt_uint8 = (uint8_t*) &n_cross_e_scal_b_re; // Affect an uint8_t pointer with the adress of n_cross_e_scal_b_re -#ifdef LSB_FIRST_TCH - lfr_bp1[i*NB_BYTES_BP1+3] = lfr_bp1[i*NB_BYTES_BP1+3] | (pt_uint8[3] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 4th octet:PC convention) - // Record it at the 8th bit position (from the right to the left) - // of lfr_bp1[i*NB_BYTES_BP1+3] - pt_uint8[3] = (pt_uint8[3] & 0x7f); // Make n_cross_e_scal_b_re be positive in any case: |n_cross_e_scal_b_re| -#endif -#ifdef MSB_FIRST_TCH - lfr_bp1[i*NB_BYTES_BP1+3] = lfr_bp1[i*NB_BYTES_BP1+3] | (pt_uint8[0] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 0th octet:SPARC convention) - // Record it at the 8th bit position (from the right to the left) - // of lfr_bp1[i*NB_BYTES_BP1+3] - pt_uint8[0] = (pt_uint8[0] & 0x7f); // Make n_cross_e_scal_b_re be positive in any case: |n_cross_e_scal_b_re| -#endif - vphi = n_cross_e_scal_b_re / bx_bx_star; // Compute |VPHI| - - significand = frexpf(vphi/2, &exponent); // 0.5 <= significand < 1 - // vphi/2 = significand * 2^exponent - // The division by 2 is to ensure that max value <= 2^30 (rough estimate) - // Should be reconsidered by taking into account the k-coefficients ... - - if (exponent < expmin) { // value should be >= 0.5 * 2^expmin - exponent = expmin; - significand = 0.5; // min value that can be recorded - } - if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1) - exponent = expmax; - significand = 1.0; // max value that can be recorded - } - if (significand == 0) {// in that case exponent == 0 too - exponent = expmin; - significand = 0.5; // min value that can be recorded - } -#ifdef DEBUG_TCH - printf("|VPHI| / 2 : %16.8e\n",vphi/2); - printf("significand : %16.8e\n",significand); - printf("exponent : %d\n" ,exponent); -#endif - lfr_bp1[i*NB_BYTES_BP1+8] = (uint8_t) ((significand*2-1)*7 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding - // where just the first 3 bits are used (0, ..., 7) - tmp_uint8 = (uint8_t) (exponent-expmin); // Shift and cast into a 8-bit uint8_t where - // just the first 5 bits are used (0, ..., 2^5-1) -#ifdef DEBUG_TCH - printf("lfr_bp1[i*NB_BYTES_BP1+8] for VPHI significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+8]); - printf("tmp_uint8 for VPHI exponent : %d\n",tmp_uint8); -#endif - lfr_bp1[i*NB_BYTES_BP1+8] = lfr_bp1[i*NB_BYTES_BP1+8] | (tmp_uint8 << 3); // shift these 5 bits to the left before logical addition - // with lfr_bp1[i*NB_BYTES_BP1+8] -#ifdef DEBUG_TCH - printf("lfr_bp1[i*NB_BYTES_BP1+8] for VPHI exponent + significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+8]); - printf("lfr_bp1[i*NB_BYTES_BP1+3] for VPHI sign + PSDB 'exponent' : %u\n",lfr_bp1[i*NB_BYTES_BP1+3]); -#endif - pt_uint8 = (uint8_t*) &n_cross_e_scal_b_im; // Affect an uint8_t pointer with the adress of n_cross_e_scal_b_im -#ifdef LSB_FIRST_TCH - pt_uint8[3] = pt_uint8[3] & 0x7f; // Make n_cross_e_scal_b_im be positive in any case: |ImaSX| -#endif -#ifdef MSB_FIRST_TCH - pt_uint8[0] = pt_uint8[0] & 0x7f; // Make n_cross_e_scal_b_im be positive in any case: |ImaSX| -#endif - tmp_uint8 = (n_cross_e_scal_b_im > n_cross_e_scal_b_re) ? 0x40 : 0x00; // Determine the sector argument of SX. If |Im| > |Re| affect - // an unsigned 8-bit char with 01000000; otherwise with null. - lfr_bp1[i*NB_BYTES_BP1+3] = lfr_bp1[i*NB_BYTES_BP1+3] | tmp_uint8; // Record it as a sign bit at the 7th bit position (from the right - // to the left) of lfr_bp1[i*NB_BYTES_BP1+3], by simple logical addition. -#ifdef DEBUG_TCH - printf("|n_cross_e_scal_b_im| : %16.8e\n",n_cross_e_scal_b_im); - printf("|n_cross_e_scal_b_im|/bx_bx_star/2: %16.8e\n",n_cross_e_scal_b_im/bx_bx_star/2); - printf("ArgNEBX sign : %u\n",tmp_uint8); - printf("lfr_bp1[i*NB_BYTES_BP1+3] for VPHI & ArgNEBX signs + PSDB 'exponent' : %u\n",lfr_bp1[i*NB_BYTES_BP1+3]); -#endif - } -} - -void BP2_set( float * compressed_spec_mat, uint8_t nb_bins_compressed_spec_mat, uint8_t * lfr_bp2 ) -{ - float cross_re; // 32-bit floating point - float cross_im; - float aux; - float significand; - int exponent; // 32-bit signed integer - uint8_t nbitexp; // 8-bit unsigned integer - uint8_t nbitsig; - uint8_t *pt_uint8; // pointer on unsigned 8-bit integer - int8_t expmin; // 8-bit signed integer - int8_t expmax; - uint16_t rangesig; // 16-bit unsigned integer - uint16_t autocor; - uint16_t exp; - uint16_t tmp_uint16; - uint16_t i; - -#ifdef DEBUG_TCH - printf("BP2 : \n"); - printf("Number of bins: %d\n", nb_bins_compressed_spec_mat); -#endif - - // For floating point data to be recorded on 16-bit words : - nbitexp = 6; // number of bits for the exponent - nbitsig = 16 - nbitexp; // number of bits for the significand - rangesig = (1 << nbitsig)-1; // == 2^nbitsig - 1 - expmax = 32; - expmin = expmax - (1 << nbitexp) + 1; - -#ifdef DEBUG_TCH - printf("nbitexp : %d, nbitsig : %d, rangesig : %d\n", nbitexp, nbitsig, rangesig); - printf("expmin : %d, expmax : %d\n", expmin, expmax); -#endif - - for(i = 0; i= 0.5 * 2^expmin - exponent = expmin; - significand = 0.5; // min value that can be recorded - } - if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1) - exponent = expmax; - significand = 1.0; // max value that can be recorded - } - if (significand == 0) { // in that case exponent == 0 too - exponent = expmin; - significand = 0.5; // min value that can be recorded - } - - autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding - // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1) - exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just - // the first nbitexp bits are used (0, ..., 2^nbitexp-1) - tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the - // left place of the significand bits (nbitsig), - // making the 16-bit word to be recorded - pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16 -#ifdef LSB_FIRST_TCH - lfr_bp2[i*NB_BYTES_BP2+0] = pt_uint8[0]; // Record LSB of tmp_uint16 - lfr_bp2[i*NB_BYTES_BP2+1] = pt_uint8[1]; // Record MSB of tmp_uint16 -#endif -#ifdef MSB_FIRST_TCH - lfr_bp2[i*NB_BYTES_BP2+0] = pt_uint8[1]; // Record LSB of tmp_uint16 - lfr_bp2[i*NB_BYTES_BP2+1] = pt_uint8[0]; // Record MSB of tmp_uint16 -#endif -#ifdef DEBUG_TCH - printf("autocor for S11 significand : %u\n",autocor); - printf("exp for S11 exponent : %u\n",exp); - printf("pt_uint8[1] for S11 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]); - printf("pt_uint8[0] for S11 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]); - printf("lfr_bp2[i*NB_BYTES_BP2+1] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+1], lfr_bp2[i*NB_BYTES_BP2+1]); - printf("lfr_bp2[i*NB_BYTES_BP2+0] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+0], lfr_bp2[i*NB_BYTES_BP2+0]); -#endif - // S22 - significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9], &exponent); // 0.5 <= significand < 1 - // S22 = significand * 2^exponent -#ifdef DEBUG_TCH - printf("S22 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]); - printf("significand : %16.8e\n",significand); - printf("exponent : %d\n" ,exponent); -#endif - if (exponent < expmin) { // value should be >= 0.5 * 2^expmin - exponent = expmin; - significand = 0.5; // min value that can be recorded - } - if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1) - exponent = expmax; - significand = 1.0; // max value that can be recorded - } - if (significand == 0) { // in that case exponent == 0 too - exponent = expmin; - significand = 0.5; // min value that can be recorded - } - - autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding - // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1) - exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just - // the first nbitexp bits are used (0, ..., 2^nbitexp-1) - tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the - // left place of the significand bits (nbitsig), - // making the 16-bit word to be recorded - pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16 -#ifdef LSB_FIRST_TCH - lfr_bp2[i*NB_BYTES_BP2+2] = pt_uint8[0]; // Record LSB of tmp_uint16 - lfr_bp2[i*NB_BYTES_BP2+3] = pt_uint8[1]; // Record MSB of tmp_uint16 -#endif -#ifdef MSB_FIRST_TCH - lfr_bp2[i*NB_BYTES_BP2+2] = pt_uint8[1]; // Record LSB of tmp_uint16 - lfr_bp2[i*NB_BYTES_BP2+3] = pt_uint8[0]; // Record MSB of tmp_uint16 -#endif -#ifdef DEBUG_TCH - printf("autocor for S22 significand : %u\n",autocor); - printf("exp for S11 exponent : %u\n",exp); - printf("pt_uint8[1] for S22 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]); - printf("pt_uint8[0] for S22 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]); - printf("lfr_bp2[i*NB_BYTES_BP2+3] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+3], lfr_bp2[i*NB_BYTES_BP2+3]); - printf("lfr_bp2[i*NB_BYTES_BP2+2] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+2], lfr_bp2[i*NB_BYTES_BP2+2]); -#endif - // S33 - significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16], &exponent); // 0.5 <= significand < 1 - // S33 = significand * 2^exponent -#ifdef DEBUG_TCH - printf("S33 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]); - printf("significand : %16.8e\n",significand); - printf("exponent : %d\n" ,exponent); -#endif - if (exponent < expmin) { // value should be >= 0.5 * 2^expmin - exponent = expmin; - significand = 0.5; // min value that can be recorded - } - if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1) - exponent = expmax; - significand = 1.0; // max value that can be recorded - } - if (significand == 0) { // in that case exponent == 0 too - exponent = expmin; - significand = 0.5; // min value that can be recorded - } - - autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding - // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1) - exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just - // the first nbitexp bits are used (0, ..., 2^nbitexp-1) - tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the - // left place of the significand bits (nbitsig), - // making the 16-bit word to be recorded - pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16 -#ifdef LSB_FIRST_TCH - lfr_bp2[i*NB_BYTES_BP2+4] = pt_uint8[0]; // Record LSB of tmp_uint16 - lfr_bp2[i*NB_BYTES_BP2+5] = pt_uint8[1]; // Record MSB of tmp_uint16 -#endif -#ifdef MSB_FIRST_TCH - lfr_bp2[i*NB_BYTES_BP2+4] = pt_uint8[1]; // Record LSB of tmp_uint16 - lfr_bp2[i*NB_BYTES_BP2+5] = pt_uint8[0]; // Record MSB of tmp_uint16 -#endif -#ifdef DEBUG_TCH - printf("autocor for S33 significand : %u\n",autocor); - printf("exp for S33 exponent : %u\n",exp); - printf("pt_uint8[1] for S33 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]); - printf("pt_uint8[0] for S33 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]); - printf("lfr_bp2[i*NB_BYTES_BP2+5] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+5], lfr_bp2[i*NB_BYTES_BP2+5]); - printf("lfr_bp2[i*NB_BYTES_BP2+4] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+4], lfr_bp2[i*NB_BYTES_BP2+4]); -#endif - // S44 - significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21], &exponent); // 0.5 <= significand < 1 - // S44 = significand * 2^exponent -#ifdef DEBUG_TCH - printf("S44 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]); - printf("significand : %16.8e\n",significand); - printf("exponent : %d\n" ,exponent); -#endif - - if (exponent < expmin) { // value should be >= 0.5 * 2^expmin - exponent = expmin; - significand = 0.5; // min value that can be recorded - } - if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1) - exponent = expmax; - significand = 1.0; // max value that can be recorded - } - if (significand == 0) { // in that case exponent == 0 too - exponent = expmin; - significand = 0.5; // min value that can be recorded - } - - autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding - // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1) - exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just - // the first nbitexp bits are used (0, ..., 2^nbitexp-1) - tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the - // left place of the significand bits (nbitsig), - // making the 16-bit word to be recorded - pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16 -#ifdef LSB_FIRST_TCH - lfr_bp2[i*NB_BYTES_BP2+6] = pt_uint8[0]; // Record LSB of tmp_uint16 - lfr_bp2[i*NB_BYTES_BP2+7] = pt_uint8[1]; // Record MSB of tmp_uint16 -#endif -#ifdef MSB_FIRST_TCH - lfr_bp2[i*NB_BYTES_BP2+6] = pt_uint8[1]; // Record LSB of tmp_uint16 - lfr_bp2[i*NB_BYTES_BP2+7] = pt_uint8[0]; // Record MSB of tmp_uint16 -#endif -#ifdef DEBUG_TCH - printf("autocor for S44 significand : %u\n",autocor); - printf("exp for S44 exponent : %u\n",exp); - printf("pt_uint8[1] for S44 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]); - printf("pt_uint8[0] for S44 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]); - printf("lfr_bp2[i*NB_BYTES_BP2+7] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+7], lfr_bp2[i*NB_BYTES_BP2+7]); - printf("lfr_bp2[i*NB_BYTES_BP2+6] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+6], lfr_bp2[i*NB_BYTES_BP2+6]); -#endif - // S55 - significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24], &exponent); // 0.5 <= significand < 1 - // S55 = significand * 2^exponent -#ifdef DEBUG_TCH - printf("S55 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]); - printf("significand : %16.8e\n",significand); - printf("exponent : %d\n" ,exponent); -#endif - if (exponent < expmin) { // value should be >= 0.5 * 2^expmin - exponent = expmin; - significand = 0.5; // min value that can be recorded - } - if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1) - exponent = expmax; - significand = 1.0; // max value that can be recorded - } - if (significand == 0) { // in that case exponent == 0 too - exponent = expmin; - significand = 0.5; // min value that can be recorded - } - - autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding - // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1) - exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just - // the first nbitexp bits are used (0, ..., 2^nbitexp-1) - tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the - // left place of the significand bits (nbitsig), - // making the 16-bit word to be recorded - pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16 -#ifdef LSB_FIRST_TCH - lfr_bp2[i*NB_BYTES_BP2+8] = pt_uint8[0]; // Record LSB of tmp_uint16 - lfr_bp2[i*NB_BYTES_BP2+9] = pt_uint8[1]; // Record MSB of tmp_uint16 -#endif -#ifdef MSB_FIRST_TCH - lfr_bp2[i*NB_BYTES_BP2+8] = pt_uint8[1]; // Record LSB of tmp_uint16 - lfr_bp2[i*NB_BYTES_BP2+9] = pt_uint8[0]; // Record MSB of tmp_uint16 -#endif -#ifdef DEBUG_TCH - printf("autocor for S55 significand : %u\n",autocor); - printf("exp for S55 exponent : %u\n",exp); - printf("pt_uint8[1] for S55 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]); - printf("pt_uint8[0] for S55 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]); - printf("lfr_bp2[i*NB_BYTES_BP2+9] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+9], lfr_bp2[i*NB_BYTES_BP2+9]); - printf("lfr_bp2[i*NB_BYTES_BP2+8] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+8], lfr_bp2[i*NB_BYTES_BP2+8]); -#endif - } -} - diff --git a/basic_parameters.h b/basic_parameters.h --- a/basic_parameters.h +++ b/basic_parameters.h @@ -7,14 +7,861 @@ #ifndef BASIC_PARAMETERS_H_INCLUDED #define BASIC_PARAMETERS_H_INCLUDED -#define NB_VALUES_PER_SPECTRAL_MATRIX 25 -#define NB_BINS_COMPRESSED_MATRIX_f0 11 +#include +#include +#include -#define NB_BYTES_BP1 9 -#define NB_BYTES_BP2 30 +#include "basic_parameters_params.h" + +float k_f0[NB_BINS_COMPRESSED_MATRIX_f0][32]; void init_k_f0( void ); -void BP1_set(float * compressed_spec_mat, unsigned char nb_bins_compressed_spec_mat, unsigned char * lfr_bp1); -void BP2_set(float * compressed_spec_mat, unsigned char nb_bins_compressed_spec_mat, unsigned char * lfr_bp2); +static inline void BP1_set(float * compressed_spec_mat, unsigned char nb_bins_compressed_spec_mat, unsigned char * lfr_bp1); +static inline void BP2_set(float * compressed_spec_mat, unsigned char nb_bins_compressed_spec_mat, unsigned char * lfr_bp2); + +//*********************************** +// STATIC INLINE FUNCTION DEFINITIONS + +void BP1_set( float * compressed_spec_mat, uint8_t nb_bins_compressed_spec_mat, uint8_t * lfr_bp1 ){ + float PSDB; // 32-bit floating point + float PSDE; + float tmp; + float NVEC_V0; + float NVEC_V1; + float NVEC_V2; + float aux; + float tr_SB_SB; + float e_cross_b_re; + float e_cross_b_im; + float n_cross_e_scal_b_re; + float n_cross_e_scal_b_im; + float ny; + float nz; + float bx_bx_star; + float vphi; + float significand; + int exponent; // 32-bit signed integer + float alpha_M; + + uint8_t nbitexp; // 8-bit unsigned integer + uint8_t nbitsig; + uint8_t tmp_uint8; + uint8_t *pt_uint8; // pointer on unsigned 8-bit integer + int8_t expmin; // 8-bit signed integer + int8_t expmax; + uint16_t rangesig; // 16-bit unsigned integer + uint16_t psd; + uint16_t exp; + uint16_t tmp_uint16; + uint16_t i; + + alpha_M = 45 * (3.1415927/180); + + init_k_f0(); + +#ifdef DEBUG_TCH + printf("BP1 : \n"); + printf("Number of bins: %d\n", nb_bins_compressed_spec_mat); +#endif + + // initialization for managing the exponents of the floating point data: + nbitexp = 5; // number of bits for the exponent + expmax = 30; // maximum value of the exponent + expmin = expmax - (1 << nbitexp) + 1; // accordingly the minimum exponent value + // for floating point data to be recorded on 12-bit words: + nbitsig = 12 - nbitexp; // number of bits for the significand + rangesig = (1 << nbitsig)-1; // == 2^nbitsig - 1 + +#ifdef DEBUG_TCH + printf("nbitexp : %d, expmax : %d, expmin : %d\n", nbitexp, expmax, expmin); + printf("nbitsig : %d, rangesig : %d\n", nbitsig, rangesig); +#endif + + for(i=0; i= 0.5 * 2^expmin + exponent = expmin; + significand = 0.5; // min value that can be recorded + } + if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1) + exponent = expmax; + significand = 1.0; // max value that can be recorded + } + if (significand == 0) { // in that case exponent == 0 too + exponent = expmin; + significand = 0.5; // min value that can be recorded + } + + psd = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding + // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1) + exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just + // the first nbitexp bits are used (0, ..., 2^nbitexp-1) + tmp_uint16 = psd | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the + // left place of the significand bits (nbitsig), + // making the 16-bit word to be recorded + pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16 +#ifdef LSB_FIRST_TCH + lfr_bp1[i*NB_BYTES_BP1+2] = pt_uint8[0]; // Record LSB of tmp_uint16 + lfr_bp1[i*NB_BYTES_BP1+3] = pt_uint8[1]; // Record MSB of tmp_uint16 +#endif +#ifdef MSB_FIRST_TCH + lfr_bp1[i*NB_BYTES_BP1+2] = pt_uint8[1]; // Record LSB of tmp_uint16 + lfr_bp1[i*NB_BYTES_BP1+3] = pt_uint8[0]; // Record MSB of tmp_uint16 +#endif +#ifdef DEBUG_TCH + printf("\nBin number: %d\n", i); + printf("PSDB / 3 : %16.8e\n",PSDB/3); + printf("significand : %16.8e\n",significand); + printf("exponent : %d\n" ,exponent); + printf("psd for PSDB significand : %d\n",psd); + printf("exp for PSDB exponent : %d\n",exp); + printf("pt_uint8[1] for PSDB exponent + significand: %.3d or %.2x\n",pt_uint8[1], pt_uint8[1]); + printf("pt_uint8[0] for PSDB exponent + significand: %.3d or %.2x\n",pt_uint8[0], pt_uint8[0]); + printf("lfr_bp1[i*NB_BYTES_BP1+3] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+3], lfr_bp1[i*NB_BYTES_BP1+3]); + printf("lfr_bp1[i*NB_BYTES_BP1+2] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+2], lfr_bp1[i*NB_BYTES_BP1+2]); +#endif + //============================================== + // BP1 PSDE == PA_LFR_SC_BP1_PE_F0 == 12 bits = 5 bits (exponent) + 7 bits (significand) + PSDE = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21] * k_f0[i][K44_PE] // S44 + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24] * k_f0[i][K55_PE] // S55 + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+22] * k_f0[i][K45_PE_RE] // S45 Re + - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+23] * k_f0[i][K45_PE_IM]; // S45 Im + + significand = frexpf(PSDE/2, &exponent); // 0.5 <= significand < 1 + // PSDE/2 = significand * 2^exponent + // the division by 2 is to ensure that max value <= 2^30 + // should be reconsidered by taking into account the k-coefficients ... + + if (exponent < expmin) { // value should be >= 0.5 * 2^expmin + exponent = expmin; + significand = 0.5; // min value that can be recorded + } + if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1) + exponent = expmax; + significand = 1.0; // max value that can be recorded + } + if (significand == 0) {// in that case exponent == 0 too + exponent = expmin; + significand = 0.5; // min value that can be recorded + } + + psd = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding + // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1) + exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just + // the first nbitexp bits are used (0, ..., 2^nbitexp-1) + tmp_uint16 = psd | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the + // left place of the significand bits (nbitsig), + // making the 16-bit word to be recorded + pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16 +#ifdef LSB_FIRST_TCH + lfr_bp1[i*NB_BYTES_BP1+0] = pt_uint8[0]; // Record LSB of tmp_uint16 + lfr_bp1[i*NB_BYTES_BP1+1] = pt_uint8[1]; // Record MSB of tmp_uint16 +#endif +#ifdef MSB_FIRST_TCH + lfr_bp1[i*NB_BYTES_BP1+0] = pt_uint8[1]; // Record LSB of tmp_uint16 + lfr_bp1[i*NB_BYTES_BP1+1] = pt_uint8[0]; // Record MSB of tmp_uint16 +#endif +#ifdef DEBUG_TCH + printf("Bin number: %d\n", i); + printf("PSDE/2 : %16.8e\n",PSDE/2); + printf("significand : %16.8e\n",significand); + printf("exponent : %d\n" ,exponent); + printf("psd for PSDE significand : %d\n",psd); + printf("exp for PSDE exponent : %d\n",exp); + printf("pt_uint8[1] for PSDE exponent + significand: %.3d or %.2x\n",pt_uint8[1], pt_uint8[1]); + printf("pt_uint8[0] for PSDE exponent + significand: %.3d or %.2x\n",pt_uint8[0], pt_uint8[0]); + printf("lfr_bp1[i*NB_BYTES_BP1+1] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+1], lfr_bp1[i*NB_BYTES_BP1+1]); + printf("lfr_bp1[i*NB_BYTES_BP1+0] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+0], lfr_bp1[i*NB_BYTES_BP1+0]); +#endif + //============================================================================== + // BP1 normal wave vector == PA_LFR_SC_BP1_NVEC_V0_F0 == 8 bits + // == PA_LFR_SC_BP1_NVEC_V1_F0 == 8 bits + // == PA_LFR_SC_BP1_NVEC_V2_F0 == 1 sign bit + tmp = sqrt( compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] //Im S12 + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] //Im S13 + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11] //Im S23 + ); + NVEC_V0 = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11]/ tmp; // S23 Im => n1 + NVEC_V1 = -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] / tmp; // S13 Im => n2 + NVEC_V2 = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] / tmp; // S12 Im => n3 + + lfr_bp1[i*NB_BYTES_BP1+4] = (uint8_t) (NVEC_V0*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding + lfr_bp1[i*NB_BYTES_BP1+5] = (uint8_t) (NVEC_V1*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding + pt_uint8 = (uint8_t*) &NVEC_V2; // Affect an uint8_t pointer with the adress of NVEC_V2 +#ifdef LSB_FIRST_TCH + lfr_bp1[i*NB_BYTES_BP1+6] = pt_uint8[3] & 0x80; // Extract the sign bit of NVEC_V2 (32-bit float, sign bit in the 4th octet:PC convention) + // Record it at the 8th bit position (from the right to the left) of lfr_bp1[i*NB_BYTES_BP1+6] +#endif +#ifdef MSB_FIRST_TCH + lfr_bp1[i*NB_BYTES_BP1+6] = pt_uint8[0] & 0x80; // Extract the sign bit of NVEC_V2 (32-bit float, sign bit in the 0th octet:SPARC convention) + // Record it at the 8th bit position (from the right to the left) of lfr_bp1[i*NB_BYTES_BP1+6] +#endif +#ifdef DEBUG_TCH + printf("NVEC_V0 : %16.8e\n",NVEC_V0); + printf("NVEC_V1 : %16.8e\n",NVEC_V1); + printf("NVEC_V2 : %16.8e\n",NVEC_V2); + printf("lfr_bp1[i*NB_BYTES_BP1+4] for NVEC_V0 : %u\n",lfr_bp1[i*NB_BYTES_BP1+4]); + printf("lfr_bp1[i*NB_BYTES_BP1+5] for NVEC_V1 : %u\n",lfr_bp1[i*NB_BYTES_BP1+5]); + printf("lfr_bp1[i*NB_BYTES_BP1+6] for NVEC_V2 : %u\n",lfr_bp1[i*NB_BYTES_BP1+6]); +#endif + //======================================================= + // BP1 ellipticity == PA_LFR_SC_BP1_ELLIP_F0 == 4 bits + aux = 2*tmp / PSDB; // Compute the ellipticity + + tmp_uint8 = (uint8_t) (aux*15 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding + // where just the first 4 bits are used (0, ..., 15) + lfr_bp1[i*NB_BYTES_BP1+6] = lfr_bp1[i*NB_BYTES_BP1+6] | (tmp_uint8 << 3); // Put these 4 bits next to the right place + // of the sign bit of NVEC_V2 (recorded + // previously in lfr_bp1[i*NB_BYTES_BP1+6]) +#ifdef DEBUG_TCH + printf("ellipticity : %16.8e\n",aux); + printf("tmp_uint8 for ellipticity : %u\n",tmp_uint8); + printf("lfr_bp1[i*NB_BYTES_BP1+6] for NVEC_V2 + ellipticity : %u\n",lfr_bp1[i*NB_BYTES_BP1+6]); +#endif + //============================================================== + // BP1 degree of polarization == PA_LFR_SC_BP1_DOP_F0 == 3 bits + tr_SB_SB = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX] + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9] + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16] + + 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+1] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+1] + + 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] + + 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+3] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+3] + + 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] + + 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10] + + 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11]; + aux = PSDB*PSDB; + tmp = ( 3*tr_SB_SB - aux ) / ( 2 * aux ); // Compute the degree of polarisation + + tmp_uint8 = (uint8_t) (tmp*7 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding + // where just the first 3 bits are used (0, ..., 7) + lfr_bp1[i*NB_BYTES_BP1+6] = lfr_bp1[i*NB_BYTES_BP1+6] | tmp_uint8; // Record these 3 bits at the 3 first bit positions + // (from the right to the left) of lfr_bp1[i*NB_BYTES_BP1+6] +#ifdef DEBUG_TCH + printf("DOP : %16.8e\n",tmp); + printf("tmp_uint8 for DOP : %u\n",tmp_uint8); + printf("lfr_bp1[i*NB_BYTES_BP1+6] for NVEC_V2 + ellipticity + DOP : %u\n",lfr_bp1[i*NB_BYTES_BP1+6]); +#endif + //======================================================================================= + // BP1 X_SO-component of the Poynting flux == PA_LFR_SC_BP1_SX_F0 == 8 (+ 2) bits + // = 5 bits (exponent) + 3 bits (significand) + // + 1 sign bit + 1 argument bit (two sectors) + e_cross_b_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_f0[i][K34_SX_RE] //S34 Re + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_f0[i][K35_SX_RE] //S35 Re + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+5] *k_f0[i][K14_SX_RE] //S14 Re + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+7] *k_f0[i][K15_SX_RE] //S15 Re + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_f0[i][K24_SX_RE] //S24 Re + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_f0[i][K25_SX_RE] //S25 Re + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_f0[i][K34_SX_IM] //S34 Im + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_f0[i][K35_SX_IM] //S35 Im + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+6] *k_f0[i][K14_SX_IM] //S14 Im + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+8] *k_f0[i][K15_SX_IM] //S15 Im + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_f0[i][K24_SX_IM] //S24 Im + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_f0[i][K25_SX_IM]; //S25 Im + // Im(S_ji) = -Im(S_ij) + // k_ji = k_ij + e_cross_b_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_f0[i][K34_SX_IM] //S34 Re + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_f0[i][K35_SX_IM] //S35 Re + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+5] *k_f0[i][K14_SX_IM] //S14 Re + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+7] *k_f0[i][K15_SX_IM] //S15 Re + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_f0[i][K24_SX_IM] //S24 Re + + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_f0[i][K25_SX_IM] //S25 Re + - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_f0[i][K34_SX_RE] //S34 Im + - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_f0[i][K35_SX_RE] //S35 Im + - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+6] *k_f0[i][K14_SX_RE] //S14 Im + - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+8] *k_f0[i][K15_SX_RE] //S15 Im + - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_f0[i][K24_SX_RE] //S24 Im + - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_f0[i][K25_SX_RE]; //S25 Im +#ifdef DEBUG_TCH + printf("ReaSX / 2 : %16.8e\n",e_cross_b_re/2); +#endif + pt_uint8 = (uint8_t*) &e_cross_b_re; // Affect an uint8_t pointer with the adress of e_cross_b_re +#ifdef LSB_FIRST_TCH + lfr_bp1[i*NB_BYTES_BP1+1] = lfr_bp1[i*NB_BYTES_BP1+1] | (pt_uint8[3] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 4th octet:PC convention) + // Record it at the 8th bit position (from the right to the left) + // of lfr_bp1[i*NB_BYTES_BP1+1] + pt_uint8[3] = (pt_uint8[3] & 0x7f); // Make e_cross_b_re be positive in any case: |ReaSX| +#endif +#ifdef MSB_FIRST_TCH + lfr_bp1[i*NB_BYTES_BP1+1] = lfr_bp1[i*NB_BYTES_BP1+1] | (pt_uint8[0] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 0th octet:SPARC convention) + // Record it at the 8th bit position (from the right to the left) + // of lfr_bp1[i*NB_BYTES_BP1+1] + pt_uint8[0] = (pt_uint8[0] & 0x7f); // Make e_cross_b_re be positive in any case: |ReaSX| +#endif + significand = frexpf(e_cross_b_re/2, &exponent); // 0.5 <= significand < 1 + // ReaSX/2 = significand * 2^exponent + // The division by 2 is to ensure that max value <= 2^30 (rough estimate) + // Should be reconsidered by taking into account the k-coefficients ... + + if (exponent < expmin) { // value should be >= 0.5 * 2^expmin + exponent = expmin; + significand = 0.5; // min value that can be recorded + } + if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1) + exponent = expmax; + significand = 1.0; // max value that can be recorded + } + if (significand == 0) { // in that case exponent == 0 too + exponent = expmin; + significand = 0.5; // min value that can be recorded + } + + lfr_bp1[i*NB_BYTES_BP1+7] = (uint8_t) ((significand*2-1)*7 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding + // where just the first 3 bits are used (0, ..., 7) + tmp_uint8 = (uint8_t) (exponent-expmin); // Shift and cast into a 8-bit uint8_t where + // just the first 5 bits are used (0, ..., 2^5-1) +#ifdef DEBUG_TCH + printf("|ReaSX| / 2 : %16.8e\n",e_cross_b_re/2); + printf("significand : %16.8e\n",significand); + printf("exponent : %d\n" ,exponent); + printf("lfr_bp1[i*NB_BYTES_BP1+7] for ReaSX significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+7]); + printf("tmp_uint8 for ReaSX exponent : %d\n",tmp_uint8); +#endif + lfr_bp1[i*NB_BYTES_BP1+7] = lfr_bp1[i*NB_BYTES_BP1+7] | (tmp_uint8 << 3); // Shift these 5 bits to the left before logical addition + // with lfr_bp1[i*NB_BYTES_BP1+7] +#ifdef DEBUG_TCH + printf("lfr_bp1[i*NB_BYTES_BP1+7] for ReaSX exponent + significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+7]); + printf("lfr_bp1[i*NB_BYTES_BP1+1] for ReaSX sign + PSDE 'exponent' : %u\n",lfr_bp1[i*NB_BYTES_BP1+1]); + printf("ImaSX / 2 : %16.8e\n",e_cross_b_im/2); +#endif + pt_uint8 = (uint8_t*) &e_cross_b_im; // Affect an uint8_t pointer with the adress of e_cross_b_im +#ifdef LSB_FIRST_TCH + pt_uint8[3] = pt_uint8[3] & 0x7f; // Make e_cross_b_im be positive in any case: |ImaSX| +#endif +#ifdef MSB_FIRST_TCH + pt_uint8[0] = pt_uint8[0] & 0x7f; // Make e_cross_b_im be positive in any case: |ImaSX| +#endif + tmp_uint8 = (e_cross_b_im > e_cross_b_re) ? 0x40 : 0x00; // Determine the sector argument of SX. If |Im| > |Re| affect + // an unsigned 8-bit char with 01000000; otherwise with null. + lfr_bp1[i*NB_BYTES_BP1+1] = lfr_bp1[i*NB_BYTES_BP1+1] | tmp_uint8; // Record it as a sign bit at the 7th bit position (from the right + // to the left) of lfr_bp1[i*NB_BYTES_BP1+1], by simple logical addition. +#ifdef DEBUG_TCH + printf("|ImaSX| / 2 : %16.8e\n",e_cross_b_im/2); + printf("ArgSX sign : %u\n",tmp_uint8); + printf("lfr_bp1[i*NB_BYTES_BP1+1] for ReaSX & ArgSX signs + PSDE 'exponent' : %u\n",lfr_bp1[i*NB_BYTES_BP1+1]); +#endif + //====================================================================== + // BP1 phase velocity estimator == PA_LFR_SC_BP1_VPHI_F0 == 8 (+ 2) bits + // = 5 bits (exponent) + 3 bits (significand) + // + 1 sign bit + 1 argument bit (two sectors) + ny = sin(alpha_M)*NVEC_V1 + cos(alpha_M)*NVEC_V2; + nz = NVEC_V0; + bx_bx_star = cos(alpha_M)*cos(alpha_M)*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9] // S22 Re + + sin(alpha_M)*sin(alpha_M)*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16] // S33 Re + - 2*sin(alpha_M)*cos(alpha_M)*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10]; // S23 Re + + n_cross_e_scal_b_re = ny * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_f0[i][K24_NY_RE] //S24 Re + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_f0[i][K25_NY_RE] //S25 Re + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_f0[i][K34_NY_RE] //S34 Re + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_f0[i][K35_NY_RE] //S35 Re + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_f0[i][K24_NY_IM] //S24 Im + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_f0[i][K25_NY_IM] //S25 Im + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_f0[i][K34_NY_IM] //S34 Im + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_f0[i][K35_NY_IM]) //S35 Im + + nz * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_f0[i][K24_NZ_RE] //S24 Re + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_f0[i][K25_NZ_RE] //S25 Re + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_f0[i][K34_NZ_RE] //S34 Re + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_f0[i][K35_NZ_RE] //S35 Re + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_f0[i][K24_NZ_IM] //S24 Im + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_f0[i][K25_NZ_IM] //S25 Im + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_f0[i][K34_NZ_IM] //S34 Im + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_f0[i][K35_NZ_IM]);//S35 Im + // Im(S_ji) = -Im(S_ij) + // k_ji = k_ij + n_cross_e_scal_b_im = ny * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_f0[i][K24_NY_IM] //S24 Re + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_f0[i][K25_NY_IM] //S25 Re + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_f0[i][K34_NY_IM] //S34 Re + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_f0[i][K35_NY_IM] //S35 Re + -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_f0[i][K24_NY_RE] //S24 Im + -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_f0[i][K25_NY_RE] //S25 Im + -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_f0[i][K34_NY_RE] //S34 Im + -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_f0[i][K35_NY_RE]) //S35 Im + + nz * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_f0[i][K24_NZ_IM] //S24 Re + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_f0[i][K25_NZ_IM] //S25 Re + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_f0[i][K34_NZ_IM] //S34 Re + +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_f0[i][K35_NZ_IM] //S35 Re + -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_f0[i][K24_NZ_RE] //S24 Im + -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_f0[i][K25_NZ_RE] //S25 Im + -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_f0[i][K34_NZ_RE] //S34 Im + -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_f0[i][K35_NZ_RE]);//S35 Im +#ifdef DEBUG_TCH + printf("n_cross_e_scal_b_re : %16.8e\n",n_cross_e_scal_b_re); + printf("n_cross_e_scal_b_im : %16.8e\n",n_cross_e_scal_b_im); +#endif + // vphi = n_cross_e_scal_b_re / bx_bx_star => sign(VPHI) = sign(n_cross_e_scal_b_re) + pt_uint8 = (uint8_t*) &n_cross_e_scal_b_re; // Affect an uint8_t pointer with the adress of n_cross_e_scal_b_re +#ifdef LSB_FIRST_TCH + lfr_bp1[i*NB_BYTES_BP1+3] = lfr_bp1[i*NB_BYTES_BP1+3] | (pt_uint8[3] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 4th octet:PC convention) + // Record it at the 8th bit position (from the right to the left) + // of lfr_bp1[i*NB_BYTES_BP1+3] + pt_uint8[3] = (pt_uint8[3] & 0x7f); // Make n_cross_e_scal_b_re be positive in any case: |n_cross_e_scal_b_re| +#endif +#ifdef MSB_FIRST_TCH + lfr_bp1[i*NB_BYTES_BP1+3] = lfr_bp1[i*NB_BYTES_BP1+3] | (pt_uint8[0] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 0th octet:SPARC convention) + // Record it at the 8th bit position (from the right to the left) + // of lfr_bp1[i*NB_BYTES_BP1+3] + pt_uint8[0] = (pt_uint8[0] & 0x7f); // Make n_cross_e_scal_b_re be positive in any case: |n_cross_e_scal_b_re| +#endif + vphi = n_cross_e_scal_b_re / bx_bx_star; // Compute |VPHI| + + significand = frexpf(vphi/2, &exponent); // 0.5 <= significand < 1 + // vphi/2 = significand * 2^exponent + // The division by 2 is to ensure that max value <= 2^30 (rough estimate) + // Should be reconsidered by taking into account the k-coefficients ... + + if (exponent < expmin) { // value should be >= 0.5 * 2^expmin + exponent = expmin; + significand = 0.5; // min value that can be recorded + } + if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1) + exponent = expmax; + significand = 1.0; // max value that can be recorded + } + if (significand == 0) {// in that case exponent == 0 too + exponent = expmin; + significand = 0.5; // min value that can be recorded + } +#ifdef DEBUG_TCH + printf("|VPHI| / 2 : %16.8e\n",vphi/2); + printf("significand : %16.8e\n",significand); + printf("exponent : %d\n" ,exponent); +#endif + lfr_bp1[i*NB_BYTES_BP1+8] = (uint8_t) ((significand*2-1)*7 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding + // where just the first 3 bits are used (0, ..., 7) + tmp_uint8 = (uint8_t) (exponent-expmin); // Shift and cast into a 8-bit uint8_t where + // just the first 5 bits are used (0, ..., 2^5-1) +#ifdef DEBUG_TCH + printf("lfr_bp1[i*NB_BYTES_BP1+8] for VPHI significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+8]); + printf("tmp_uint8 for VPHI exponent : %d\n",tmp_uint8); +#endif + lfr_bp1[i*NB_BYTES_BP1+8] = lfr_bp1[i*NB_BYTES_BP1+8] | (tmp_uint8 << 3); // shift these 5 bits to the left before logical addition + // with lfr_bp1[i*NB_BYTES_BP1+8] +#ifdef DEBUG_TCH + printf("lfr_bp1[i*NB_BYTES_BP1+8] for VPHI exponent + significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+8]); + printf("lfr_bp1[i*NB_BYTES_BP1+3] for VPHI sign + PSDB 'exponent' : %u\n",lfr_bp1[i*NB_BYTES_BP1+3]); +#endif + pt_uint8 = (uint8_t*) &n_cross_e_scal_b_im; // Affect an uint8_t pointer with the adress of n_cross_e_scal_b_im +#ifdef LSB_FIRST_TCH + pt_uint8[3] = pt_uint8[3] & 0x7f; // Make n_cross_e_scal_b_im be positive in any case: |ImaSX| +#endif +#ifdef MSB_FIRST_TCH + pt_uint8[0] = pt_uint8[0] & 0x7f; // Make n_cross_e_scal_b_im be positive in any case: |ImaSX| +#endif + tmp_uint8 = (n_cross_e_scal_b_im > n_cross_e_scal_b_re) ? 0x40 : 0x00; // Determine the sector argument of SX. If |Im| > |Re| affect + // an unsigned 8-bit char with 01000000; otherwise with null. + lfr_bp1[i*NB_BYTES_BP1+3] = lfr_bp1[i*NB_BYTES_BP1+3] | tmp_uint8; // Record it as a sign bit at the 7th bit position (from the right + // to the left) of lfr_bp1[i*NB_BYTES_BP1+3], by simple logical addition. +#ifdef DEBUG_TCH + printf("|n_cross_e_scal_b_im| : %16.8e\n",n_cross_e_scal_b_im); + printf("|n_cross_e_scal_b_im|/bx_bx_star/2: %16.8e\n",n_cross_e_scal_b_im/bx_bx_star/2); + printf("ArgNEBX sign : %u\n",tmp_uint8); + printf("lfr_bp1[i*NB_BYTES_BP1+3] for VPHI & ArgNEBX signs + PSDB 'exponent' : %u\n",lfr_bp1[i*NB_BYTES_BP1+3]); +#endif + } +} + +void BP2_set( float * compressed_spec_mat, uint8_t nb_bins_compressed_spec_mat, uint8_t * lfr_bp2 ) +{ + float cross_re; // 32-bit floating point + float cross_im; + float aux; + float significand; + int exponent; // 32-bit signed integer + uint8_t nbitexp; // 8-bit unsigned integer + uint8_t nbitsig; + uint8_t *pt_uint8; // pointer on unsigned 8-bit integer + int8_t expmin; // 8-bit signed integer + int8_t expmax; + uint16_t rangesig; // 16-bit unsigned integer + uint16_t autocor; + uint16_t exp; + uint16_t tmp_uint16; + uint16_t i; + +#ifdef DEBUG_TCH + printf("BP2 : \n"); + printf("Number of bins: %d\n", nb_bins_compressed_spec_mat); +#endif + + // For floating point data to be recorded on 16-bit words : + nbitexp = 6; // number of bits for the exponent + nbitsig = 16 - nbitexp; // number of bits for the significand + rangesig = (1 << nbitsig)-1; // == 2^nbitsig - 1 + expmax = 32; + expmin = expmax - (1 << nbitexp) + 1; + +#ifdef DEBUG_TCH + printf("nbitexp : %d, nbitsig : %d, rangesig : %d\n", nbitexp, nbitsig, rangesig); + printf("expmin : %d, expmax : %d\n", expmin, expmax); +#endif + + for(i = 0; i= 0.5 * 2^expmin + exponent = expmin; + significand = 0.5; // min value that can be recorded + } + if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1) + exponent = expmax; + significand = 1.0; // max value that can be recorded + } + if (significand == 0) { // in that case exponent == 0 too + exponent = expmin; + significand = 0.5; // min value that can be recorded + } + + autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding + // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1) + exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just + // the first nbitexp bits are used (0, ..., 2^nbitexp-1) + tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the + // left place of the significand bits (nbitsig), + // making the 16-bit word to be recorded + pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16 +#ifdef LSB_FIRST_TCH + lfr_bp2[i*NB_BYTES_BP2+0] = pt_uint8[0]; // Record LSB of tmp_uint16 + lfr_bp2[i*NB_BYTES_BP2+1] = pt_uint8[1]; // Record MSB of tmp_uint16 +#endif +#ifdef MSB_FIRST_TCH + lfr_bp2[i*NB_BYTES_BP2+0] = pt_uint8[1]; // Record LSB of tmp_uint16 + lfr_bp2[i*NB_BYTES_BP2+1] = pt_uint8[0]; // Record MSB of tmp_uint16 +#endif +#ifdef DEBUG_TCH + printf("autocor for S11 significand : %u\n",autocor); + printf("exp for S11 exponent : %u\n",exp); + printf("pt_uint8[1] for S11 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]); + printf("pt_uint8[0] for S11 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]); + printf("lfr_bp2[i*NB_BYTES_BP2+1] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+1], lfr_bp2[i*NB_BYTES_BP2+1]); + printf("lfr_bp2[i*NB_BYTES_BP2+0] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+0], lfr_bp2[i*NB_BYTES_BP2+0]); +#endif + // S22 + significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9], &exponent); // 0.5 <= significand < 1 + // S22 = significand * 2^exponent +#ifdef DEBUG_TCH + printf("S22 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]); + printf("significand : %16.8e\n",significand); + printf("exponent : %d\n" ,exponent); +#endif + if (exponent < expmin) { // value should be >= 0.5 * 2^expmin + exponent = expmin; + significand = 0.5; // min value that can be recorded + } + if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1) + exponent = expmax; + significand = 1.0; // max value that can be recorded + } + if (significand == 0) { // in that case exponent == 0 too + exponent = expmin; + significand = 0.5; // min value that can be recorded + } + + autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding + // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1) + exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just + // the first nbitexp bits are used (0, ..., 2^nbitexp-1) + tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the + // left place of the significand bits (nbitsig), + // making the 16-bit word to be recorded + pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16 +#ifdef LSB_FIRST_TCH + lfr_bp2[i*NB_BYTES_BP2+2] = pt_uint8[0]; // Record LSB of tmp_uint16 + lfr_bp2[i*NB_BYTES_BP2+3] = pt_uint8[1]; // Record MSB of tmp_uint16 +#endif +#ifdef MSB_FIRST_TCH + lfr_bp2[i*NB_BYTES_BP2+2] = pt_uint8[1]; // Record LSB of tmp_uint16 + lfr_bp2[i*NB_BYTES_BP2+3] = pt_uint8[0]; // Record MSB of tmp_uint16 +#endif +#ifdef DEBUG_TCH + printf("autocor for S22 significand : %u\n",autocor); + printf("exp for S11 exponent : %u\n",exp); + printf("pt_uint8[1] for S22 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]); + printf("pt_uint8[0] for S22 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]); + printf("lfr_bp2[i*NB_BYTES_BP2+3] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+3], lfr_bp2[i*NB_BYTES_BP2+3]); + printf("lfr_bp2[i*NB_BYTES_BP2+2] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+2], lfr_bp2[i*NB_BYTES_BP2+2]); +#endif + // S33 + significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16], &exponent); // 0.5 <= significand < 1 + // S33 = significand * 2^exponent +#ifdef DEBUG_TCH + printf("S33 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]); + printf("significand : %16.8e\n",significand); + printf("exponent : %d\n" ,exponent); +#endif + if (exponent < expmin) { // value should be >= 0.5 * 2^expmin + exponent = expmin; + significand = 0.5; // min value that can be recorded + } + if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1) + exponent = expmax; + significand = 1.0; // max value that can be recorded + } + if (significand == 0) { // in that case exponent == 0 too + exponent = expmin; + significand = 0.5; // min value that can be recorded + } + + autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding + // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1) + exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just + // the first nbitexp bits are used (0, ..., 2^nbitexp-1) + tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the + // left place of the significand bits (nbitsig), + // making the 16-bit word to be recorded + pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16 +#ifdef LSB_FIRST_TCH + lfr_bp2[i*NB_BYTES_BP2+4] = pt_uint8[0]; // Record LSB of tmp_uint16 + lfr_bp2[i*NB_BYTES_BP2+5] = pt_uint8[1]; // Record MSB of tmp_uint16 +#endif +#ifdef MSB_FIRST_TCH + lfr_bp2[i*NB_BYTES_BP2+4] = pt_uint8[1]; // Record LSB of tmp_uint16 + lfr_bp2[i*NB_BYTES_BP2+5] = pt_uint8[0]; // Record MSB of tmp_uint16 +#endif +#ifdef DEBUG_TCH + printf("autocor for S33 significand : %u\n",autocor); + printf("exp for S33 exponent : %u\n",exp); + printf("pt_uint8[1] for S33 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]); + printf("pt_uint8[0] for S33 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]); + printf("lfr_bp2[i*NB_BYTES_BP2+5] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+5], lfr_bp2[i*NB_BYTES_BP2+5]); + printf("lfr_bp2[i*NB_BYTES_BP2+4] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+4], lfr_bp2[i*NB_BYTES_BP2+4]); +#endif + // S44 + significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21], &exponent); // 0.5 <= significand < 1 + // S44 = significand * 2^exponent +#ifdef DEBUG_TCH + printf("S44 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]); + printf("significand : %16.8e\n",significand); + printf("exponent : %d\n" ,exponent); +#endif + + if (exponent < expmin) { // value should be >= 0.5 * 2^expmin + exponent = expmin; + significand = 0.5; // min value that can be recorded + } + if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1) + exponent = expmax; + significand = 1.0; // max value that can be recorded + } + if (significand == 0) { // in that case exponent == 0 too + exponent = expmin; + significand = 0.5; // min value that can be recorded + } + + autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding + // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1) + exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just + // the first nbitexp bits are used (0, ..., 2^nbitexp-1) + tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the + // left place of the significand bits (nbitsig), + // making the 16-bit word to be recorded + pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16 +#ifdef LSB_FIRST_TCH + lfr_bp2[i*NB_BYTES_BP2+6] = pt_uint8[0]; // Record LSB of tmp_uint16 + lfr_bp2[i*NB_BYTES_BP2+7] = pt_uint8[1]; // Record MSB of tmp_uint16 +#endif +#ifdef MSB_FIRST_TCH + lfr_bp2[i*NB_BYTES_BP2+6] = pt_uint8[1]; // Record LSB of tmp_uint16 + lfr_bp2[i*NB_BYTES_BP2+7] = pt_uint8[0]; // Record MSB of tmp_uint16 +#endif +#ifdef DEBUG_TCH + printf("autocor for S44 significand : %u\n",autocor); + printf("exp for S44 exponent : %u\n",exp); + printf("pt_uint8[1] for S44 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]); + printf("pt_uint8[0] for S44 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]); + printf("lfr_bp2[i*NB_BYTES_BP2+7] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+7], lfr_bp2[i*NB_BYTES_BP2+7]); + printf("lfr_bp2[i*NB_BYTES_BP2+6] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+6], lfr_bp2[i*NB_BYTES_BP2+6]); +#endif + // S55 + significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24], &exponent); // 0.5 <= significand < 1 + // S55 = significand * 2^exponent +#ifdef DEBUG_TCH + printf("S55 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]); + printf("significand : %16.8e\n",significand); + printf("exponent : %d\n" ,exponent); +#endif + if (exponent < expmin) { // value should be >= 0.5 * 2^expmin + exponent = expmin; + significand = 0.5; // min value that can be recorded + } + if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1) + exponent = expmax; + significand = 1.0; // max value that can be recorded + } + if (significand == 0) { // in that case exponent == 0 too + exponent = expmin; + significand = 0.5; // min value that can be recorded + } + + autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding + // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1) + exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just + // the first nbitexp bits are used (0, ..., 2^nbitexp-1) + tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the + // left place of the significand bits (nbitsig), + // making the 16-bit word to be recorded + pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16 +#ifdef LSB_FIRST_TCH + lfr_bp2[i*NB_BYTES_BP2+8] = pt_uint8[0]; // Record LSB of tmp_uint16 + lfr_bp2[i*NB_BYTES_BP2+9] = pt_uint8[1]; // Record MSB of tmp_uint16 +#endif +#ifdef MSB_FIRST_TCH + lfr_bp2[i*NB_BYTES_BP2+8] = pt_uint8[1]; // Record LSB of tmp_uint16 + lfr_bp2[i*NB_BYTES_BP2+9] = pt_uint8[0]; // Record MSB of tmp_uint16 +#endif +#ifdef DEBUG_TCH + printf("autocor for S55 significand : %u\n",autocor); + printf("exp for S55 exponent : %u\n",exp); + printf("pt_uint8[1] for S55 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]); + printf("pt_uint8[0] for S55 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]); + printf("lfr_bp2[i*NB_BYTES_BP2+9] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+9], lfr_bp2[i*NB_BYTES_BP2+9]); + printf("lfr_bp2[i*NB_BYTES_BP2+8] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+8], lfr_bp2[i*NB_BYTES_BP2+8]); +#endif + } +} + #endif // BASIC_PARAMETERS_H_INCLUDED diff --git a/basic_parameters_params.h b/basic_parameters_params.h new file mode 100755 --- /dev/null +++ b/basic_parameters_params.h @@ -0,0 +1,48 @@ +#ifndef BASIC_PARAMETERS_PARAMS_H +#define BASIC_PARAMETERS_PARAMS_H + +#define NB_VALUES_PER_SPECTRAL_MATRIX 25 +#define NB_BINS_COMPRESSED_MATRIX_f0 11 + +#define NB_BYTES_BP1 9 +#define NB_BYTES_BP2 30 + +//** +// K +#define K44_PE 0 +#define K55_PE 1 +#define K45_PE_RE 2 +#define K45_PE_IM 3 + +#define K14_SX_RE 4 +#define K14_SX_IM 5 +#define K15_SX_RE 6 +#define K15_SX_IM 7 +#define K24_SX_RE 8 +#define K24_SX_IM 9 +#define K25_SX_RE 10 +#define K25_SX_IM 11 +#define K34_SX_RE 12 +#define K34_SX_IM 13 +#define K35_SX_RE 14 +#define K35_SX_IM 15 + +#define K24_NY_RE 16 +#define K24_NY_IM 17 +#define K25_NY_RE 18 +#define K25_NY_IM 19 +#define K34_NY_RE 20 +#define K34_NY_IM 21 +#define K35_NY_RE 22 +#define K35_NY_IM 23 + +#define K24_NZ_RE 24 +#define K24_NZ_IM 25 +#define K25_NZ_RE 26 +#define K25_NZ_IM 27 +#define K34_NZ_RE 28 +#define K34_NZ_IM 29 +#define K35_NZ_RE 30 +#define K35_NZ_IM 31 + +#endif // BASIC_PARAMETERS_PARAMS_H diff --git a/tests3.pro b/tests3.pro --- a/tests3.pro +++ b/tests3.pro @@ -1,18 +1,19 @@ -TEMPLATE = app -CONFIG += console -CONFIG -= app_bundle -CONFIG -= qt - -DEFINES += DEBUG_TCH -#DEFINES += MSB_FIRST_TCH # SPARC convention -DEFINES += LSB_FIRST_TCH # PC convention - -SOURCES += main.c \ - basic_parameters.c \ - file_utilities.c - -HEADERS += \ - basic_parameters.h \ - file_utilities.h - - +TEMPLATE = app +CONFIG += console +CONFIG -= app_bundle +CONFIG -= qt + +DEFINES += DEBUG_TCH +#DEFINES += MSB_FIRST_TCH # SPARC convention +DEFINES += LSB_FIRST_TCH # PC convention + +SOURCES += main.c \ + basic_parameters.c \ + file_utilities.c + +HEADERS += \ + basic_parameters.h \ + file_utilities.h \ + basic_parameters_params.h + +