Skip to content

GitLab

Commit bc550f69 authored by OP's avatar OP
Browse files

Ajout d'options pour calculer les cross-correlations VV et HV dans visiavg

parent 10062417
Hide whitespace changes
Inline Side-by-side
Showing with 281 additions and 5 deletions
p4autils.cc
View file @ bc550f69
......@@ -194,6 +194,17 @@ std::vector<sa_size_t> P4AVisiNumEncoder::getAllVCrossCor()
}
return P4AVisiNumEncoder::Convert2VisiNumber(nms);
}
std::vector<sa_size_t> P4AVisiNumEncoder::getAllHVCrossCor()
{
const char * vcxnm[16]={"1H-1V","1H-2V","1H-3V","1H-4V","2H-1V","2H-2V","2H-3V","2H-4V","3H-1V","3H-2V","3H-3V","3H-4V","4H-1V","4H-2V","4H-3V","4H-4V"};
std::vector<std::string> nms;
for(size_t i=0; i<16; i++) {
nms.push_back(std::string(vcxnm[i]));
}
return P4AVisiNumEncoder::Convert2VisiNumber(nms);
}
//---------------------------------------------------------------------
//------------------------ classe P4AnaParams ------------------------
//---------------------------------------------------------------------
......@@ -203,7 +214,7 @@ P4AnaParams::P4AnaParams()
gain_gnu_file_(""), gain_var_file_(""), gain_blind_filt_file_(""), phase_corr_param_file_(""),
outfile_("p4a.ppf"), fitsoutfile_(""), prtlev_(0), prtmodulo_(1),
Nmean_(10), fgTFM_(false), TFMtimebin_(30), TFMfreqbin_(4), fgdtable_(false),
fgRA_(false), RAcenter_(12.), RAwidth_(24.), RAresol_mn_(1.)
fgRA_(false), RAcenter_(12.), RAwidth_(24.), RAresol_mn_(1.), doVV_(false),doHV_(false)
{
}
......@@ -217,7 +228,7 @@ P4AnaParams::P4AnaParams(P4AnaParams const & a)
outfile_(a.outfile_), fitsoutfile_(a.fitsoutfile_), prtlev_(a.prtlev_), prtmodulo_(a.prtmodulo_),
Nmean_(a.Nmean_), fgTFM_(a.fgTFM_), TFMtimebin_(a.TFMtimebin_), TFMfreqbin_(a.TFMfreqbin_),
fgdtable_(a.fgdtable_), fbands_(a.fbands_),
fgRA_(a.fgRA_), RAcenter_(a.RAcenter_), RAwidth_(a.RAwidth_), RAresol_mn_(a.RAresol_mn_),
fgRA_(a.fgRA_), RAcenter_(a.RAcenter_), RAwidth_(a.RAwidth_), RAresol_mn_(a.RAresol_mn_),doVV_(a.doVV_),doHV_(a.doHV_),
lastargs_(a.lastargs_)
{
}
......@@ -234,7 +245,7 @@ P4AnaParams& P4AnaParams::operator = (P4AnaParams const & a)
prtlev_=a.prtlev_; prtmodulo_=a.prtmodulo_;
Nmean_=a.Nmean_; fgTFM_=a.fgTFM_; TFMfreqbin_=a.TFMfreqbin_; TFMtimebin_=a.TFMtimebin_;
fgdtable_=a.fgdtable_; fbands_=a.fbands_;
fgRA_=a.fgRA_; RAcenter_=a.RAcenter_; RAwidth_=a.RAwidth_; RAresol_mn_=a.RAresol_mn_;
fgRA_=a.fgRA_; RAcenter_=a.RAcenter_; RAwidth_=a.RAwidth_; RAresol_mn_=a.RAresol_mn_;doVV_=a.doVV_;doHV_=a.doHV_;
lastargs_=a.lastargs_;
return *this;
}
......@@ -275,6 +286,14 @@ int P4AnaParams::DecodeArgs(int narg, const char* arg[])
else if (fbo=="-reorderfreq") {
fgreorderfreq_=true; arg++; narg--; lastargs_.clear();
}
else if (fbo=="-dovv") {
cout << "P4AnaParams::DecodeArgs : VV cross correlations demanded "<<endl;
doVV_=true; arg++; narg--; lastargs_.clear();
}
else if (fbo=="-dohv") {
cout << "P4AnaParams::DecodeArgs : HV cross correlations demanded - may need lof of memory ! "<<endl;
doHV_=true; arg++; narg--; lastargs_.clear();
}
else if (fbo=="-gnu") {
if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl; return -1; }
gain_gnu_file_=arg[2]; arg+=2; narg-=2; lastargs_.clear();
......
p4autils.h
View file @ bc550f69
......@@ -135,6 +135,8 @@ public:
static std::vector<sa_size_t> getAllHCrossCor();
//! return a vector with visibility number (VisiMatrix rows) for all V-V type cross correlations ( 6 cx-corr)
static std::vector<sa_size_t> getAllVCrossCor();
//! return a vector with visibility number (VisiMatrix rows) for all H-V type cross correlations ( 16 cx-corr)
static std::vector<sa_size_t> getAllHVCrossCor();
//! prints a list of all visibility names and corresponding row numbers
ostream& PrintAll(ostream & os) const;
......@@ -213,6 +215,9 @@ public:
double RAcenter_, RAwidth_, RAresol_mn_;
// les arguments en fin de ligne de commande
vector<string> lastargs_;
bool doVV_; // true : also computes VxV cross correlations tf maps (visiavg)
bool doHV_; // true : also computes HxV cross correlations tf maps (visiavg)
};
inline ostream& operator << (ostream& os, P4AnaParams const& a)
......
visiavg.cc
View file @ bc550f69
......@@ -91,12 +91,15 @@ int main(int narg, const char* arg[])
P4AVisiNumEncoder visiencod;
vector<sa_size_t> KVAC = visiencod.getAllAutoCor();
vector<sa_size_t> KVCXHH = visiencod.getAllHCrossCor();
vector<sa_size_t> KVCXVV = visiencod.getAllVCrossCor();
vector<sa_size_t> KVCXHV = visiencod.getAllHVCrossCor();
cout << " List of AutoCorrelation rows:"<<endl;
for(size_t k=0; k<KVAC.size(); k++) {
cout << "KVAC["<<k<<"]="<<KVAC[k]<<" ->"<<visiencod.Convert2VisiName(KVAC[k])<<endl;
}
cout << " List of HH X-cor rows:"<<endl;
for(size_t k=0; k<KVCXHH.size(); k++) {
cout << "KVCXHH["<<k<<"]="<<KVCXHH[k]<<" ->"<<visiencod.Convert2VisiName(KVCXHH[k])<<endl;
......@@ -129,23 +132,49 @@ int main(int narg, const char* arg[])
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;
// for VV if needed
TMatrix< complex<r_4> > cxsum_vv;
// for sums of real and imag parts
TMatrix< r_4 > cxsum_vv_sq_rp;
TMatrix< r_4 > cxsum_vv_sq_ip;
// for HV if needed
TMatrix< complex<r_4> > cxsum_hv;
// for sums of real and imag parts
TMatrix< r_4 > cxsum_hv_sq_rp;
TMatrix< r_4 > cxsum_hv_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;
//----- 6 V-V cross-cor TimeFrequency maps
vector< TArray< complex<r_4> > > vtfm_vv;
//----- 6 V-V cross-cor TimeFrequency maps for the variances of real and imag parts
vector< TArray< r_4 > > vtfm_vv_rp_sq;
vector< TArray< r_4 > > vtfm_vv_ip_sq;
//----- 16 H-V cross-cor TimeFrequency maps
vector< TArray< complex<r_4> > > vtfm_hv;
//----- 16 H-V cross-cor TimeFrequency maps for the variances of real and imag parts
vector< TArray< r_4 > > vtfm_hv_rp_sq;
vector< TArray< r_4 > > vtfm_hv_ip_sq;
//----- 8 auto-corr TimeFrequency maps
vector< TArray< r_4 > > vtfmac;
//----- 8 auto-corr TimeFrequency variance maps
......@@ -164,12 +193,29 @@ int main(int narg, const char* arg[])
vismtx = wreader.getAverageVisMtx(cfdate);
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());
// VV if needed
if (params.doVV_){
cxsum_vv.SetSize(6, vismtx.NCols());
cxsum_vv_sq_rp.SetSize(6, vismtx.NCols());
cxsum_vv_sq_ip.SetSize(6, vismtx.NCols());
}
// HV if needed
if (params.doHV_){
cxsum_hv.SetSize(16, vismtx.NCols());
cxsum_hv_sq_rp.SetSize(16, vismtx.NCols());
cxsum_hv_sq_ip.SetSize(16, vismtx.NCols());
}
tfmSX=wreader.getTotalNbWindows()/deltaIavg;
tfmSY=vismtx.NCols()/TFMfbin;
//allocating 8 Auto-Corr time-frequency maps
......@@ -184,6 +230,26 @@ int main(int narg, const char* arg[])
cout << "and the 2x6 for squares of real & imaginary 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) );
// VV if needed
if (params.doVV_){
//allocating 6 Cross-Corr V-V time-frequency maps
cout<<"visiavg/Info: allocating V-V cross-cor Time-Frequency maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
for(int k=0; k<6; k++) vtfm_vv.push_back( TArray< complex<r_4> >(tfmSX, tfmSY) );
cout << "and the 2x6 for squares of real & imaginary parts " << endl;
for(int k=0; k<6; k++) vtfm_vv_rp_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
for(int k=0; k<6; k++) vtfm_vv_ip_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
}
// HV if needed
if (params.doHV_){
//allocating 16 Cross-Corr H-V time-frequency maps
cout<<"visiavg/Info: allocating 16 H-V cross-cor Time-Frequency maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
for(int k=0; k<16; k++) vtfm_hv.push_back( TArray< complex<r_4> >(tfmSX, tfmSY) );
cout << "and the 2x16 for squares of real & imaginary parts " << endl;
for(int k=0; k<16; k++) vtfm_hv_rp_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
for(int k=0; k<16; k++) vtfm_hv_ip_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
}
// recupere le jour de depart @ 0h
datestart = TimeStamp(cfdate.DaysPart(),0.);
......@@ -198,6 +264,21 @@ int main(int narg, const char* arg[])
cxsum = complex<r_4>(0.,0.);
cxsum_sq_rp = 0.;
cxsum_sq_ip = 0.;
// VV if needed
if (params.doVV_){
cxsum_vv = complex<r_4>(0.,0.);
cxsum_vv_sq_rp = 0.;
cxsum_vv_sq_ip = 0.;
}
// HV if needed
if (params.doHV_){
cxsum_hv = complex<r_4>(0.,0.);
cxsum_hv_sq_rp = 0.;
cxsum_hv_sq_ip = 0.;
}
}
// sum (integration) along the time axis
......@@ -217,6 +298,33 @@ int main(int narg, const char* arg[])
cxsum_sq_ip.Row(k) += tmp.MulElt(tmp,tmp) ;
}
if (params.doVV_){
for(size_t k=0; k<KVCXVV.size(); k++){
cxsum_vv.Row(k) += vismtx.Row(KVCXVV[k]); // les cross-correlations
TVector<r_4> tmp = real(vismtx.Row(KVCXVV[k]));
cxsum_vv_sq_rp.Row(k) += tmp.MulElt(tmp,tmp) ;
tmp = imag(vismtx.Row(KVCXVV[k]));
cxsum_vv_sq_ip.Row(k) += tmp.MulElt(tmp,tmp) ;
}
}
if (params.doHV_){
for(size_t k=0; k<KVCXHV.size(); k++){
cxsum_hv.Row(k) += vismtx.Row(KVCXHV[k]); // les cross-correlations HV
TVector<r_4> tmp = real(vismtx.Row(KVCXHV[k]));
cxsum_hv_sq_rp.Row(k) += tmp.MulElt(tmp,tmp) ;
tmp = imag(vismtx.Row(KVCXHV[k]));
cxsum_hv_sq_ip.Row(k) += tmp.MulElt(tmp,tmp) ;
}
}
I++; // incrementing DeltaTime counter
......@@ -254,7 +362,47 @@ int main(int narg, const char* arg[])
}
}
} //----- end of loop over the 6 Xcor
if (params.doVV_){ // Option VV
//---- On s'occupe des 6 cross-correlations 1V-2V ... 3V-4V
for(size_t k=0; k<KVCXVV.size(); k++) { // loop over the 6 Xcor
TVector< complex<r_4> > vcx = cxsum_vv.Row(k);
TVector<r_4> vcxprsq = cxsum_vv_sq_rp.Row(k);
TVector<r_4> vcxpisq = cxsum_vv_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_vv[k];
TArray< r_4 > & tfmapsqpr = vtfm_vv_rp_sq[k];
TArray< r_4 > & tfmapsqpi = vtfm_vv_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 VV
} // end VV option
if (params.doHV_){ // Option HV
//---- On s'occupe des 16 cross-correlations 1H-1V ... 4H-4V
for(size_t k=0; k<KVCXHV.size(); k++) { // loop over the 16 Xcor
TVector< complex<r_4> > vcx = cxsum_hv.Row(k);
TVector<r_4> vcxprsq = cxsum_hv_sq_rp.Row(k);
TVector<r_4> vcxpisq = cxsum_hv_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_hv[k];
TArray< r_4 > & tfmapsqpr = vtfm_hv_rp_sq[k];
TArray< r_4 > & tfmapsqpi = vtfm_hv_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 16 Xcor HV
} // end HV option
double tdif = cfdate.TimeDifferenceSeconds(cfdate,dateobs)/2.;
timevec(TFMtmidx) = dateobs.TimeDifferenceSeconds(dateobs.ShiftSeconds (tdif ),datestart); // centre du bin
ravec(TFMtmidx) = P4Coords::RAFromTimeTU(dateobs.ShiftSeconds (tdif ));
......@@ -346,6 +494,110 @@ int main(int narg, const char* arg[])
}
}
// ========================== end saving HH XCorr
//===========================
if (params.doVV_){ // Option VV
// --- renormalizing and saving V-V Cross-Corr time-frequency maps
cout<<" visiavg/Info: Saving 6 V-V cross-corr time-frequency maps to PPF file "<<outfile<<endl;
const char* tfmVV_names[6]={"TFM_1V2V", "TFM_1V3V", "TFM_1V4V", "TFM_2V3V", "TFM_2V4V", "TFM_3V4V"};
const char* vrtfmVV_names[6]={"RVARTFM_1V2V", "RVARTFM_1V3V", "RVARTFM_1V4V", "RVARTFM_2V3V", "RVARTFM_2V4V", "RVARTFM_3V4V"};
const char* vitfmVV_names[6]={"IVARTFM_1V2V", "IVARTFM_1V3V", "IVARTFM_1V4V", "IVARTFM_2V3V", "IVARTFM_2V4V", "IVARTFM_3V4V"};
for(int k=0; k<6; k++) { // loop over the 6 Xcor
TArray< complex<r_4> > & tfmap = vtfm_vv[k];
TArray< r_4 > & tfmap_sqpr = vtfm_vv_rp_sq[k];
TArray< r_4 > & tfmap_sqpi = vtfm_vv_ip_sq[k];
tfmap *= complex<r_4>((r_4)(1./((double)deltaIavg*(double)TFMfbin)), 0.);
potfm << PPFNameTag(tfmVV_names[k]) << tfmap;
if (fos != NULL) {
ext_names.push_back(string(tfmVV_names[k])+"_real");
(*fos)<< real(tfmap);
ext_names.push_back(string(tfmVV_names[k])+"_imag");
(*fos)<< imag(tfmap);
}
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(vrtfmVV_names[k]) << tfmap_sqpr;
potfm << PPFNameTag(vitfmVV_names[k]) << tfmap_sqpi;
if (fos != NULL) {
ext_names.push_back(vrtfmVV_names[k]);
(*fos)<< tfmap_sqpr;
ext_names.push_back(vitfmVV_names[k]);
(*fos)<< tfmap_sqpi;
}
}
} // end option VV
//===========================
if (params.doHV_){ // Option HV
// --- renormalizing and saving H-V Cross-Corr time-frequency maps
cout<<" visiavg/Info: Saving 16 H-V cross-corr time-frequency maps to PPF file "<<outfile<<endl;
const char* tfmHV_names[16]={"TFM_1H1V", "TFM_1H2V", "TFM_1H3V", "TFM_1H4V",
"TFM_2H1V", "TFM_2H2V", "TFM_2H3V", "TFM_2H4V",
"TFM_3H1V", "TFM_3H2V", "TFM_3H3V", "TFM_3H4V",
"TFM_4H1V", "TFM_4H2V", "TFM_4H3V", "TFM_4H4V"
};
const char* vrtfmHV_names[16]={"RVARTFM_1H1V", "RVARTFM_1H2V", "RVARTFM_1H3V", "RVARTFM_1H4V",
"RVARTFM_2H1V", "RVARTFM_2H2V", "RVARTFM_2H3V", "RVARTFM_2H4V",
"RVARTFM_3H1V", "RVARTFM_3H2V", "RVARTFM_3H3V", "RVARTFM_3H4V",
"RVARTFM_4H1V", "RVARTFM_4H2V", "RVARTFM_4H3V", "RVARTFM_4H4V"
};
const char* vitfmHV_names[16]={"IVARTFM_1H1V", "IVARTFM_1H2V", "IVARTFM_1H3V", "IVARTFM_1H4V",
"IVARTFM_2H1V", "IVARTFM_2H2V", "IVARTFM_2H3V", "IVARTFM_2H4V",
"IVARTFM_3H1V", "IVARTFM_3H2V", "IVARTFM_3H3V", "IVARTFM_3H4V",
"IVARTFM_4H1V", "IVARTFM_4H2V", "IVARTFM_4H3V", "IVARTFM_4H4V"
};
for(int k=0; k<16; k++) { // loop over the 6 Xcor
TArray< complex<r_4> > & tfmap = vtfm_hv[k];
TArray< r_4 > & tfmap_sqpr = vtfm_hv_rp_sq[k];
TArray< r_4 > & tfmap_sqpi = vtfm_hv_ip_sq[k];
tfmap *= complex<r_4>((r_4)(1./((double)deltaIavg*(double)TFMfbin)), 0.);
potfm << PPFNameTag(tfmHV_names[k]) << tfmap;
if (fos != NULL) {
ext_names.push_back(string(tfmHV_names[k])+"_real");
(*fos)<< real(tfmap);
ext_names.push_back(string(tfmHV_names[k])+"_imag");
(*fos)<< imag(tfmap);
}
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(vrtfmHV_names[k]) << tfmap_sqpr;
potfm << PPFNameTag(vitfmHV_names[k]) << tfmap_sqpi;
if (fos != NULL) {
ext_names.push_back(vrtfmHV_names[k]);
(*fos)<< tfmap_sqpr;
ext_names.push_back(vitfmHV_names[k]);
(*fos)<< tfmap_sqpi;
}
}
} // end option HV
potfm << PPFNameTag("TimeVec") << timevec ;
potfm << PPFNameTag("RAVec") << ravec ;
......
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