##// END OF EJS Templates
version 1.2 qui règle le problème d'alignement mémoire pour BP1 (pour BP2 plus tard ...)
version 1.2 qui règle le problème d'alignement mémoire pour BP1 (pour BP2 plus tard ...)

File last commit:

r6:be100cfb14bd default
r6:be100cfb14bd default
Show More
basic_parameters.c
881 lines | 57.8 KiB | text/x-c | CLexer
// In the frame of RPW LFR Sofware ICD Issue1 Rev8 (05/07/2013)
// version 1.0: 31/07/2013
// version 1.1: 02/04/2014
// version 1.2: 30/04/2014
#include "basic_parameters.h"
#include <math.h>
#include <stdio.h>
#include <stdint.h>
#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 )
{
uint16_t i; // 16 bits unsigned
for(i=0; i<NB_BINS_COMPRESSED_MATRIX_f0; i++){
k_f0[i][K44_PE] = 1;
k_f0[i][K55_PE] = 1;
k_f0[i][K45_PE_RE] = 1;
k_f0[i][K45_PE_IM] = 1;
k_f0[i][K14_SX_RE] = 1;
k_f0[i][K14_SX_IM] = 1;
k_f0[i][K15_SX_RE] = 1;
k_f0[i][K15_SX_IM] = 1;
k_f0[i][K24_SX_RE] = 1;
k_f0[i][K24_SX_IM] = 1;
k_f0[i][K25_SX_RE] = 1;
k_f0[i][K25_SX_IM] = 1;
k_f0[i][K34_SX_RE] = 1;
k_f0[i][K34_SX_IM] = 1;
k_f0[i][K35_SX_RE] = 1;
k_f0[i][K35_SX_IM] = 1;
k_f0[i][K24_NY_RE] = 1;
k_f0[i][K24_NY_IM] = 1;
k_f0[i][K25_NY_RE] = 1;
k_f0[i][K25_NY_IM] = 1;
k_f0[i][K34_NY_RE] = 1;
k_f0[i][K34_NY_IM] = 1;
k_f0[i][K35_NY_RE] = 1;
k_f0[i][K35_NY_IM] = 1;
k_f0[i][K24_NZ_RE] = 1;
k_f0[i][K24_NZ_IM] = 1;
k_f0[i][K25_NZ_RE] = 1;
k_f0[i][K25_NZ_IM] = 1;
k_f0[i][K34_NZ_RE] = 1;
k_f0[i][K34_NZ_IM] = 1;
k_f0[i][K35_NZ_RE] = 1;
k_f0[i][K35_NZ_IM] = 1;
}
}
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 ){
int exponent; // 32 bits signed
float PSDB; // 32 bits 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;
uint8_t nbitexp; // 8 bits unsigned
uint8_t nbitsig;
uint8_t tmp_uint8;
uint8_t *pt_uint8; // pointer on unsigned 8-bit bytes
int8_t expmin; // 8 bits signed
int8_t expmax;
int16_t rangesig; // 16 bits unsigned
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<nb_bins_compressed_spec_mat; i++){
//==============================================
// BP1 PSDB == PA_LFR_SC_BP1_PB_F0 == 12 bits = 5 bits (exponent) + 7 bits (significand)
PSDB = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX] // S11
+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9] // S22
+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]; // S33
significand = frexpf(PSDB/3, &exponent); // 0.5 <= significand < 1
// PSDB/3 = significand * 2^exponent
// the division by 3 is to ensure that max value <= 2^30
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+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 )
{
int i, exponent;
float aux, significand, cross_re, cross_im;
int8_t nbitexp, nbitsig, expmin, expmax; // 8 bits
int16_t rangesig; // 16 bits
uint16_t autocor, tmp_uint16; // 16 bits
uint16_t *pt_u_short_int; // pointer on unsigned 16-bit words
#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<nb_bins_compressed_spec_mat; i++){
//==============================================
// BP2 normalized cross correlations == PA_LFR_SC_BP2_CROSS_F0 == 10 * (8+8) bits
// == PA_LFR_SC_BP2_CROSS_RE_0_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_IM_0_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_RE_1_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_IM_1_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_RE_2_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_IM_2_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_RE_3_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_IM_3_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_RE_4_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_IM_4_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_RE_5_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_IM_5_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_RE_6_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_IM_6_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_RE_7_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_IM_7_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_RE_8_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_IM_8_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_RE_9_F0 == 8 bits
// == PA_LFR_SC_BP2_CROSS_IM_9_F0 == 8 bits
// S12
aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]);
cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+1] / aux;
cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] / aux;
lfr_bp2[i*NB_BYTES_BP2+10] = (uint8_t) (cross_re*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
lfr_bp2[i*NB_BYTES_BP2+20] = (uint8_t) (cross_im*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
#ifdef DEBUG_TCH
printf("\nBin number: %d\n", i);
printf("lfr_bp2[i*NB_BYTES_BP2+10] for cross12_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+10]);
printf("lfr_bp2[i*NB_BYTES_BP2+20] for cross12_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+20]);
#endif
// S13
aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]);
cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+3] / aux;
cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] / aux;
lfr_bp2[i*NB_BYTES_BP2+11] = (uint8_t) (cross_re*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
lfr_bp2[i*NB_BYTES_BP2+21] = (uint8_t) (cross_im*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
#ifdef DEBUG_TCH
printf("lfr_bp2[i*NB_BYTES_BP2+11] for cross13_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+11]);
printf("lfr_bp2[i*NB_BYTES_BP2+21] for cross13_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+21]);
#endif
// S14
aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]);
cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+5] / aux;
cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+6] / aux;
lfr_bp2[i*NB_BYTES_BP2+12] = (uint8_t) (cross_re*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
lfr_bp2[i*NB_BYTES_BP2+22] = (uint8_t) (cross_im*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
#ifdef DEBUG_TCH
printf("lfr_bp2[i*NB_BYTES_BP2+12] for cross14_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+12]);
printf("lfr_bp2[i*NB_BYTES_BP2+22] for cross14_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+22]);
#endif
// S15
aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);
cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+7] / aux;
cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+8] / aux;
lfr_bp2[i*NB_BYTES_BP2+13] = (uint8_t) (cross_re*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
lfr_bp2[i*NB_BYTES_BP2+23] = (uint8_t) (cross_im*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
#ifdef DEBUG_TCH
printf("lfr_bp2[i*NB_BYTES_BP2+13] for cross15_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+13]);
printf("lfr_bp2[i*NB_BYTES_BP2+23] for cross15_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+23]);
#endif
// S23
aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]);
cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10] / aux;
cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11] / aux;
lfr_bp2[i*NB_BYTES_BP2+14] = (uint8_t) (cross_re*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
lfr_bp2[i*NB_BYTES_BP2+24] = (uint8_t) (cross_im*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
#ifdef DEBUG_TCH
printf("lfr_bp2[i*NB_BYTES_BP2+14] for cross23_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+14]);
printf("lfr_bp2[i*NB_BYTES_BP2+24] for cross23_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+24]);
#endif
// S24
aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]);
cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12] / aux;
cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13] / aux;
lfr_bp2[i*NB_BYTES_BP2+15] = (uint8_t) (cross_re*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
lfr_bp2[i*NB_BYTES_BP2+25] = (uint8_t) (cross_im*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
#ifdef DEBUG_TCH
printf("lfr_bp2[i*NB_BYTES_BP2+15] for cross24_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+15]);
printf("lfr_bp2[i*NB_BYTES_BP2+25] for cross24_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+25]);
#endif
// S25
aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);
cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14] / aux;
cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15] / aux;
lfr_bp2[i*NB_BYTES_BP2+16] = (uint8_t) (cross_re*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
lfr_bp2[i*NB_BYTES_BP2+26] = (uint8_t) (cross_im*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
#ifdef DEBUG_TCH
printf("lfr_bp2[i*NB_BYTES_BP2+16] for cross25_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+16]);
printf("lfr_bp2[i*NB_BYTES_BP2+26] for cross25_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+26]);
#endif
// S34
aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]);
cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17] / aux;
cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18] / aux;
lfr_bp2[i*NB_BYTES_BP2+17] = (uint8_t) (cross_re*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
lfr_bp2[i*NB_BYTES_BP2+27] = (uint8_t) (cross_im*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
#ifdef DEBUG_TCH
printf("lfr_bp2[i*NB_BYTES_BP2+17] for cross34_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+17]);
printf("lfr_bp2[i*NB_BYTES_BP2+27] for cross34_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+27]);
#endif
// S35
aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);
cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19] / aux;
cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20] / aux;
lfr_bp2[i*NB_BYTES_BP2+18] = (uint8_t) (cross_re*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
lfr_bp2[i*NB_BYTES_BP2+28] = (uint8_t) (cross_im*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
#ifdef DEBUG_TCH
printf("lfr_bp2[i*NB_BYTES_BP2+18] for cross35_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+18]);
printf("lfr_bp2[i*NB_BYTES_BP2+28] for cross35_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+28]);
#endif
// S45
aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);
cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+22] / aux;
cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+23] / aux;
lfr_bp2[i*NB_BYTES_BP2+19] = (uint8_t) (cross_re*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
lfr_bp2[i*NB_BYTES_BP2+29] = (uint8_t) (cross_im*127.5 + 128); // shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
#ifdef DEBUG_TCH
printf("lfr_bp2[i*NB_BYTES_BP2+19] for cross45_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+19]);
printf("lfr_bp2[i*NB_BYTES_BP2+29] for cross45_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+29]);
#endif
//==============================================
// BP2 auto correlations == PA_LFR_SC_BP2_AUTO_F0 == 5*16 bits = 5*[6 bits (exponent) + 10 bits (significand)]
// == PA_LFR_SC_BP2_AUTO_A0_F0 == 16 bits
// == PA_LFR_SC_BP2_AUTO_A1_F0 == 16 bits
// == PA_LFR_SC_BP2_AUTO_A2_F0 == 16 bits
// == PA_LFR_SC_BP2_AUTO_A3_F0 == 16 bits
// == PA_LFR_SC_BP2_AUTO_A4_F0 == 16 bits
// S11
significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX], &exponent); // 0.5 <= significand < 1
// S11 = significand * 2^exponent
#ifdef DEBUG_TCH
printf("S11 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]);
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)
tmp_uint16 = (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)
pt_u_short_int = (uint16_t*) &lfr_bp2[i*NB_BYTES_BP2+0]; // Affect an uint16_t pointer with the
// adress where the 16-bit word result will be stored
*pt_u_short_int = autocor | (tmp_uint16 << nbitsig); // Put the exponent bits (nbitexp) next to the
// left place of the significand bits (nbitsig), making
// the 16-bit word to be recorded, and record it using the pointer
#ifdef DEBUG_TCH
printf("autocor for S11 significand : %u\n",autocor );
printf("tmp_uint8 for S11 exponent : %u\n",tmp_uint16 );
printf("*pt_u_short_int for S11 exponent + significand : %.3d or %x\n",*pt_u_short_int, *pt_u_short_int);
printf("lfr_bp2[i*NB_BYTES_BP2+1] : %u or %x\n",lfr_bp2[i*NB_BYTES_BP2+1], lfr_bp2[i*NB_BYTES_BP2+1]);
printf("lfr_bp2[i*NB_BYTES_BP2+0] : %u or %x\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)
tmp_uint16 = (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)
pt_u_short_int = (uint16_t*) &lfr_bp2[i*NB_BYTES_BP2+2]; // Affect an uint16_t pointer with the
// adress where the 16-bit word result will be stored
*pt_u_short_int = autocor | (tmp_uint16 << nbitsig); // Put the exponent bits (nbitexp) next to the
// left place of the significand bits (nbitsig), making
// the 16-bit word to be recorded, and record it using the pointer
#ifdef DEBUG_TCH
printf("autocor for S22 significand : %d\n",autocor );
printf("tmp_uint8 for S22 exponent : %d\n",tmp_uint16 );
printf("*pt_u_short_int for S22 exponent + significand : %.3d or %x\n",*pt_u_short_int, *pt_u_short_int);
printf("lfr_bp2[i*NB_BYTES_BP2+3] : %.3d or %x\n",lfr_bp2[i*NB_BYTES_BP2+3], lfr_bp2[i*NB_BYTES_BP2+3]);
printf("lfr_bp2[i*NB_BYTES_BP2+2] : %.3d or %x\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)
tmp_uint16 = (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)
pt_u_short_int = (uint16_t*) &lfr_bp2[i*NB_BYTES_BP2+4]; // Affect an uint16_t pointer with the
// adress where the 16-bit word result will be stored
*pt_u_short_int = autocor | (tmp_uint16 << nbitsig); // Put the exponent bits (nbitexp) next to the
// left place of the significand bits (nbitsig), making
// the 16-bit word to be recorded, and record it using the pointer
#ifdef DEBUG_TCH
printf("autocor for S33 significand : %d\n",autocor );
printf("tmp_uint8 for S33 exponent : %d\n",tmp_uint16 );
printf("*pt_u_short_int for S33 exponent + significand : %.3d or %x\n",*pt_u_short_int, *pt_u_short_int);
printf("lfr_bp2[i*NB_BYTES_BP2+5] : %.3d or %x\n",lfr_bp2[i*NB_BYTES_BP2+5], lfr_bp2[i*NB_BYTES_BP2+5]);
printf("lfr_bp2[i*NB_BYTES_BP2+4] : %.3d or %x\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)
tmp_uint16 = (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)
pt_u_short_int = (uint16_t*) &lfr_bp2[i*NB_BYTES_BP2+6]; // Affect an uint16_t pointer with the
// adress where the 16-bit word result will be stored
*pt_u_short_int = autocor | (tmp_uint16 << nbitsig); // Put the exponent bits (nbitexp) next to the
// left place of the significand bits (nbitsig), making
// the 16-bit word to be recorded, and record it using the pointer
#ifdef DEBUG_TCH
printf("autocor for S44 significand : %d\n",autocor );
printf("tmp_uint8 for S44 exponent : %d\n",tmp_uint16 );
printf("*pt_u_short_int for S44 exponent + significand : %.3d or %x\n",*pt_u_short_int, *pt_u_short_int);
printf("lfr_bp2[i*NB_BYTES_BP2+7] : %.3d or %x\n",lfr_bp2[i*NB_BYTES_BP2+7], lfr_bp2[i*NB_BYTES_BP2+7]);
printf("lfr_bp2[i*NB_BYTES_BP2+6] : %.3d or %x\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)
tmp_uint16 = (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)
pt_u_short_int = (uint16_t*) &lfr_bp2[i*NB_BYTES_BP2+8]; // Affect an uint16_t pointer with the
// adress where the 16-bit word result will be stored
*pt_u_short_int = autocor | (tmp_uint16 << nbitsig); // Put the exponent bits (nbitexp) next to the
// left place of the significand bits (nbitsig), making
// the 16-bit word to be recorded, and record it using the pointer
#ifdef DEBUG_TCH
printf("autocor for S55 significand : %d\n",autocor );
printf("tmp_uint8 for S55 exponent : %d\n",tmp_uint16 );
printf("*pt_u_short_int for S55 exponent + significand : %.3d or %x\n",*pt_u_short_int, *pt_u_short_int);
printf("lfr_bp2[i*NB_BYTES_BP2+9] : %.3d or %x\n",lfr_bp2[i*NB_BYTES_BP2+9], lfr_bp2[i*NB_BYTES_BP2+9]);
printf("lfr_bp2[i*NB_BYTES_BP2+8] : %.3d or %x\n",lfr_bp2[i*NB_BYTES_BP2+8], lfr_bp2[i*NB_BYTES_BP2+8]);
#endif
}
}