Commit beab735f authored by OP's avatar OP
Browse files

add variance time-freq maps

parent 3e2fca74
......@@ -137,38 +137,62 @@ int main(int narg, const char* arg[])
TVector< double > timevec((Imax-Imin)/Istep/deltaIavg);
TMatrix< complex<r_4> > vismtx;
TMatrix< complex<r_4> > acsum;
TMatrix< r_4 > acsum_sq; // wil sum only the real part ^2
TMatrix< complex<r_4> > cxsum;
// for sums of real and imag parts
TMatrix< r_4 > cxsum_sq_rp;
TMatrix< r_4 > cxsum_sq_ip;
TimeStamp dateobs, cfdate,datestart;
TimeStamp dateorg(2015,1,1,12,0,0.); // Date origine 1 jan 2015
double mttag;
int cnt=0, cntnt=0, pcntnt=0;
int I=0;
// for
//----- 6 H-H cross-cor TimeFrequency maps
vector< TArray< complex<r_4> > > vtfm;
//----- 6 H-H cross-cor TimeFrequency maps for the variances of real and imag parts
vector< TArray< r_4 > > vtfm_rp_sq;
vector< TArray< r_4 > > vtfm_ip_sq;
//----- 8 auto-corr TimeFrequency maps
vector< TArray< r_4 > > vtfmac;
//----- 8 auto-corr TimeFrequency variance maps
vector< TArray< r_4 > > vtfmac_sq;
//---- for the time-freqency map filling
sa_size_t TFMtmidx=0;
sa_size_t tfmSX, tfmSY;
while (fgok) {
//reads next visimtx
fgok=vreader.ReadNext(vismtx, cfdate, mttag);
if (!fgok) break;
// Apply gain g(nu)
p4g.applyGain(vismtx);
if (cnt==0) { //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations
acsum.SetSize(8, vismtx.NCols());
acsum_sq.SetSize(8, vismtx.NCols());
cxsum.SetSize(6, vismtx.NCols());
cxsum_sq_rp.SetSize(6, vismtx.NCols());
cxsum_sq_ip.SetSize(6, vismtx.NCols());
tfmSX=(Imax-Imin)/Istep/deltaIavg;
tfmSY=vismtx.NCols()/TFMfbin;
//allocating 8 Auto-Corr time-frequency maps
cout<<"visiavg/Info: allocating 8 AutoCor Time-Frequency maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
for(int k=0; k<8; k++) vtfmac.push_back( TArray< r_4 >(tfmSX, tfmSY) );
cout <<" and 8 for the the variance maps "<<endl;
for(int k=0; k<8; k++) vtfmac_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
//allocating 6 Cross-Corr H-H time-frequency maps
cout<<"visiavg/Info: allocating H-H cross-cor Time-Frequency maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
for(int k=0; k<6; k++) vtfm.push_back( TArray< complex<r_4> >(tfmSX, tfmSY) );
cout << "and the 2x6 for squares of rela & imag parts " << endl;
for(int k=0; k<6; k++) vtfm_rp_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
for(int k=0; k<6; k++) vtfm_ip_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
// recupere le jour de depart @ 0h
datestart = TimeStamp(cfdate.DaysPart(),0.);
......@@ -179,32 +203,62 @@ int main(int narg, const char* arg[])
if (prtlev>0)
cout<<"visiavg/Info: dateobs="<<dateobs<<" SecondsPart()="<<dateobs.SecondsPart()<<endl;
acsum = complex<r_4>(0.,0.);
acsum_sq = 0.;
cxsum = complex<r_4>(0.,0.);
cxsum_sq_rp = 0.;
cxsum_sq_ip = 0.;
}
// sum (integration) along the time axis
for(size_t k=0; k<KVAC.size(); k++) acsum.Row(k) += vismtx.Row(KVAC[k]); // Les auto-correlations
for(size_t k=0; k<KVCXHH.size(); k++) cxsum.Row(k) += vismtx.Row(KVCXHH[k]); // les cross-correlations
for(size_t k=0; k<KVAC.size(); k++) {
TVector<r_4> tmp = real(vismtx.Row(KVAC[k]));
acsum_sq.Row(k) += tmp.MulElt(tmp,tmp) ;
}
for(size_t k=0; k<KVCXHH.size(); k++){
cxsum.Row(k) += vismtx.Row(KVCXHH[k]); // les cross-correlations
TVector<r_4> tmp = real(vismtx.Row(KVCXHH[k]));
cxsum_sq_rp.Row(k) += tmp.MulElt(tmp,tmp) ;
tmp = imag(vismtx.Row(KVCXHH[k]));
cxsum_sq_ip.Row(k) += tmp.MulElt(tmp,tmp) ;
}
I++; // incrementing DeltaTime counter
if (I==deltaIavg) {
//---- On s'occupe d'abord des autocorrelations P1 ... P8
for(size_t k=0; k<KVAC.size(); k++) { // Loop over the 8 auto-correlations
TVector<r_4> vac = real(acsum.Row(k));
TVector<r_4> vacsq = acsum_sq.Row(k);
if (TFMtmidx<tfmSX) { // we check that our time index did not go beyond the allocated array size (might not be necessary)
TArray< r_4 > & tfmap = vtfmac[k];
TArray< r_4 > & tfmap_sq = vtfmac_sq[k];
for(sa_size_t jy=0; jy<tfmSY; jy++) { // frequency binning
tfmap(TFMtmidx, jy) = vac( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
tfmap_sq(TFMtmidx, jy) = vacsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
}
}
}
} //----- end of loop over the 8 AutoCor
//---- On s'occupe des cross-correlations 1H-2H ... 3H-4H
//---- On s'occupe des 6 cross-correlations 1H-2H ... 3H-4H
for(size_t k=0; k<KVCXHH.size(); k++) { // loop over the 6 Xcor
TVector< complex<r_4> > vcx = cxsum.Row(k);
TVector<r_4> vcxprsq = cxsum_sq_rp.Row(k);
TVector<r_4> vcxpisq = cxsum_sq_ip.Row(k);
if (TFMtmidx<tfmSX) { // we check that our time index did not go beyond the allocated array size (might not be necessary)
TArray< complex<r_4> > & tfmap = vtfm[k];
TArray< r_4 > & tfmapsqpr = vtfm_rp_sq[k];
TArray< r_4 > & tfmapsqpi = vtfm_ip_sq[k];
for(sa_size_t jy=0; jy<tfmSY; jy++) {
tfmap(TFMtmidx, jy) = vcx( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
tfmapsqpr(TFMtmidx, jy) = vcxprsq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
tfmapsqpi(TFMtmidx, jy) = vcxpisq( Range(jy*TFMfbin, (jy+1)*TFMfbin-1) ).Sum();
}
}
} //----- end of loop over the 6 Xcor
......@@ -226,18 +280,42 @@ int main(int narg, const char* arg[])
// --- renormalizing and saving AutoCorr time-frequency maps
cout<<" visiavg/Info: Saving 8 AutoCorr time-frequency maps to PPF file "<<outfile<<endl;
const char* tfm_names[8]={"TFM_1H", "TFM_2H", "TFM_3H", "TFM_4H", "TFM_1V", "TFM_2V", "TFM_3V", "TFM_4V"};
const char* tfmsq_names[8]={"VARTFM_1H", "VARTFM_2H", "VARTFM_3H", "VARTFM_4H", "VARTFM_1V", "VARTFM_2V", "VARTFM_3V", "VARTFM_4V"};
for(int k=0; k<8; k++) { // loop over the 8 AutoCorr
TArray< r_4 > & tfmap = vtfmac[k];
tfmap *= (r_4)(1./((double)deltaIavg*(double)TFMfbin));
TArray< r_4 > & tfmapsq = vtfmac_sq[k];
tfmapsq *= (r_4)(1./((double)deltaIavg*(double)TFMfbin));
tfmapsq = tfmapsq - tfmap.MulElt(tfmap,tfmap) ;
potfm << PPFNameTag(tfm_names[k]) << tfmap;
potfm << PPFNameTag(tfmsq_names[k]) << tfmapsq;
}
// --- renormalizing and saving H-H Cross-Corr time-frequency maps
cout<<" visiavg/Info: Saving 6 H-H cross-corr time-frequency maps to PPF file "<<outfile<<endl;
const char* tfmCC_names[6]={"TFM_1H2H", "TFM_1H3H", "TFM_1H4H", "TFM_2H3H", "TFM_2H4H", "TFM_3H4H"};
const char* vrtfmCC_names[6]={"RVARTFM_1H2H", "RVARTFM_1H3H", "RVARTFM_1H4H", "RVARTFM_2H3H", "RVARTFM_2H4H", "RVARTFM_3H4H"};
const char* vitfmCC_names[6]={"IVARTFM_1H2H", "IVARTFM_1H3H", "IVARTFM_1H4H", "IVARTFM_2H3H", "IVARTFM_2H4H", "IVARTFM_3H4H"};
for(int k=0; k<6; k++) { // loop over the 6 Xcor
TArray< complex<r_4> > & tfmap = vtfm[k];
TArray< r_4 > & tfmap_sqpr = vtfm_rp_sq[k];
TArray< r_4 > & tfmap_sqpi = vtfm_ip_sq[k];
tfmap *= complex<r_4>((r_4)(1./((double)deltaIavg*(double)TFMfbin)), 0.);
TArray< r_4 > tfmapr = real(tfmap);
TArray< r_4 > tfmapi = imag(tfmap);
tfmap_sqpr *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
tfmap_sqpi *= ((r_4)(1./((double)deltaIavg*(double)TFMfbin)));
tfmapr = tfmapr.MulElt(tfmapr,tfmapr) ;
tfmap_sqpr -= tfmapr ;
tfmapi = tfmapi.MulElt(tfmapi,tfmapi) ;
tfmap_sqpi -= tfmapi;
potfm << PPFNameTag(tfmCC_names[k]) << tfmap;
potfm << PPFNameTag(vrtfmCC_names[k]) << tfmap_sqpr;
potfm << PPFNameTag(vitfmCC_names[k]) << tfmap_sqpi;
}
potfm << PPFNameTag("TimeVec") << timevec ;
// --- FIN sauvegarde cartes temps-frequence
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment