Commit 8b05e6cb authored by Reza Ansari's avatar Reza Ansari
Browse files

1/ Implementation lecture Cross-Frequency correlation matrix ds VisiP4ReaderBase

2/ Normalisation des matrices de visibilites par NPAQSUM lors de la lecture
3/ Ajout option lecture Cross-Frequency correlation et moyannage ds rdvisip4.cc
                           Reza 17/10/2017
parent beab735f
......@@ -44,8 +44,11 @@
int Usage(void)
{
cout<<"--- rdvisip4.cc : Read PPF files produced by mfacq time-frequency\n"<<endl;
cout<<"Usage: rdvisip4 [-arguments] [-meanfft]\n";
cout<<" -meanfft : try to read mean FFT coeff and subtract it from visib. \n";
cout<<"Usage: rdvisip4 [-options] [meanfft] [cxfreq] \n";
cout<<" meanfft : try to read mean FFT coeff and subtract it from visib. \n";
cout<<" cxfreq : try to read cross-frequency correlation matrix and compute mean of it \n";
cout<<" WARNING: Only one of the two options (meanfft , cxfreq) can be specified \n";
P4AnaParams::UsageOptions();
cout<<endl;
return 1;
......@@ -66,15 +69,19 @@ int main(int narg, const char* arg[])
int prtlev=params.prtlev_;
bool fgreorderfreq = params.fgreorderfreq_;
bool fgrdmean = false;
if( (params.lastargs_.size()>0) &&(params.lastargs_[0]==string("-m")) ) fgrdmean = true;
bool fgrdcxfreq = false;
for (size_t i=0; i<params.lastargs_.size(); i++) {
if ( params.lastargs_[i]==string("meanfft") ) fgrdmean = true;
if ( params.lastargs_[i]==string("cxfreq") ) fgrdcxfreq = true;
}
params.Print(cout);
cout<<"rdvisip4/Info: Path BAO5:"<<path5<<" BAO6:"<<path6<<"\n"
<<"fgreorderfreq="<<params.fgreorderfreq_<<"\n"
<<"Imin,max,step="<<Imin<<","<<Imax<<","<<Istep<<"\n"
<<"PrtLev="<<prtlev<<" OutFile:"<<outfile<<endl;
if(fgrdmean) cout<<"Trying to read also the FFT mean coef matrix"<<endl;
if(fgrdcxfreq) cout<<"Trying to read also the Cross-Frequency correlation "<<endl;
P4AVisiNumEncoder vencod;
cout << vencod;
// ---
......@@ -112,20 +119,30 @@ int main(int narg, const char* arg[])
//vismtx: visib matrice, meanmtx: mean FFT coef, vismtx_center: visib matrix corrigee de l'offset (normalisee avec NPAQSUM)
TMatrix< complex<r_4> > vismtx, meanmtx, vismtx_center;
TMatrix< complex<r_4> > vismtx_mean, meanmtx_mean, vismtx_center_mean;
TMatrix< complex<r_4> > cxfreqmtx , cxfreqmtx_mean;
TimeStamp dateobs;
double mttag;
int cnt=0;
while (fgok) {
if(fgrdmean) {
fgok=vreader.ReadNext(vismtx, meanmtx, dateobs, mttag);
if (fgrdcxfreq) {
fgok=vreader.ReadNextCxFreq(vismtx, cxfreqmtx, dateobs, mttag);
}
else if (fgrdmean) {
fgok=vreader.ReadNextMFFT(vismtx, meanmtx, dateobs, mttag);
vismtx_center = vreader.SubtractOffset(vismtx, meanmtx);
} else {
}
else {
fgok=vreader.ReadNext(vismtx, dateobs, mttag);
}
if (fgok) {
if (cnt==0) vismtx_mean=vismtx;
else vismtx_mean+=vismtx;
if(fgrdmean) {
if (fgrdcxfreq) {
if (cnt==0) cxfreqmtx_mean=cxfreqmtx;
else cxfreqmtx_mean+=cxfreqmtx;
}
if (fgrdmean) {
if (cnt==0) meanmtx_mean = meanmtx;
else meanmtx_mean += meanmtx;
if (cnt==0) vismtx_center_mean = vismtx_center;
......@@ -142,15 +159,20 @@ int main(int narg, const char* arg[])
P4gnuGain::computeSaveGain(vismtx_mean, gainfile);
cout<<" rdvisip4/Info: Saving vismtx_mean to PPF file "<<outfile<<endl;
POutPersist po(outfile);
po<<vismtx_mean;
po<<PPFNameTag("visimean")<<vismtx_mean;
if (fgrdcxfreq) {
cxfreqmtx_mean /= complex<r_4>((r_4)cnt, (r_4)0.);
cout<<" rdvisip4/Info: Saving cxfreqmtx_mean to PPF file "<<outfile<<endl;
po<<PPFNameTag("cxfreqmean")<<cxfreqmtx_mean;
}
if(fgrdmean){
meanmtx_mean /= complex<r_4>((r_4)cnt, (r_4)0.);
cout<<" rdvisip4/Info: Saving meanmtx_mean to PPF file "<<outfile<<endl;
po<<meanmtx_mean;
po<<PPFNameTag("fftmean")<<meanmtx_mean;
vismtx_center_mean /= complex<r_4>((r_4)cnt, (r_4)0.);
cout<<" rdvisip4/Info: Saving vismtx_center_mean to PPF file "<<outfile<<endl;
po<<vismtx_center_mean;
po<<PPFNameTag("visicentermean")<<vismtx_center_mean;
}
}
......
......@@ -107,7 +107,7 @@ void VisiP4ReaderBase::SelectAll()
/* --Methode-- */
bool VisiP4ReaderBase::Read_P(long sernum, TMatrix< complex<r_4> >* vismtx, TMatrix< complex<r_4> >* meanmtx,
TimeStamp& dateobs, double& mttag)
TMatrix< complex<r_4> >* freqcxmtx, TimeStamp& dateobs, double& mttag)
{
if (paths_.size() != NBAND) {
cout << " VisiP4ReaderBase::Read_P() / ERROR: paths_.size() != NBAND ..."<<endl;
......@@ -120,11 +120,14 @@ bool VisiP4ReaderBase::Read_P(long sernum, TMatrix< complex<r_4> >* vismtx, TMat
}
bool fgwithmean = false;
if(meanmtx) fgwithmean = true;
if (meanmtx) fgwithmean = true;
bool fgwithfreqcx = false;
if (freqcxmtx) fgwithfreqcx = true;
if ( (prtlevel>1)||((prtlevel>0)&&(sernum%20==0)) ) {
cout << "---- VisiP4ReaderBase::Read_P() - sernum= "<<sernum<<" (LastNum="<<ser_last_<<")"
<< (fgwithmean ? " Read also mean coef ":" ") <<endl;
<< (fgwithmean ? " + MeanFFTCoef ":" ")
<< (fgwithfreqcx ? " + FreqCorrVii ":" ") <<endl;
}
vector<double> vmtt;
......@@ -132,18 +135,24 @@ bool VisiP4ReaderBase::Read_P(long sernum, TMatrix< complex<r_4> >* vismtx, TMat
vismtx->ReSize(NVIS, NFREQ);
if(fgwithmean) meanmtx->ReSize(NCHAN, NFREQ);
if(fgwithfreqcx) freqcxmtx->ReSize(NCHAN, NFREQ);
TMatrix< complex<r_4> > svismtx, smeanmtx;
TMatrix< complex<r_4> > svismtx, smeanmtx, sfreqcxmtx;
TMatrix< complex<r_4> > * vmtxp = vismtx;
TMatrix< complex<r_4> > * mmtxp = meanmtx;
TMatrix< complex<r_4> > * fcxmtxp = freqcxmtx;
if (fgreorderfreq_) { // s'il faut reordonner les frequences, on a besoin d'une matrice tampon
svismtx.ReSize(NVIS, NFREQ);
vmtxp = &svismtx;
if(fgwithmean){
if (fgwithmean){
smeanmtx.ReSize(NCHAN, NFREQ);
mmtxp = &smeanmtx;
}
if (fgwithfreqcx){
sfreqcxmtx.ReSize(NCHAN, NFREQ);
fcxmtxp = &sfreqcxmtx;
}
}
TMatrix< complex<r_4> > & myvismtx = (*vmtxp);
// loop over the files created by each VisCalculator (each file contains a set of visibilities)
......@@ -160,22 +169,36 @@ bool VisiP4ReaderBase::Read_P(long sernum, TMatrix< complex<r_4> >* vismtx, TMat
cout << " VisiP4ReaderBase::Read_P()() - Readding VisMtx with SerialNum= "<<sernum
<<" freqRange="<<cvf<<"-"<<cvl;
{
TMatrix< complex<r_4> > vmtx, mmtx;
TMatrix< complex<r_4> > * pmmtx=(fgwithmean && (g == 0) )?&mmtx:NULL;
ReadVisMtx(b,g,sernum,vmtx,pmmtx);
TMatrix< complex<r_4> > vmtx, mmtx, fcxmtx;
TMatrix< complex<r_4> > * pmmtx=(fgwithmean && (g == 0) )?&mmtx:NULL;
TMatrix< complex<r_4> > * pfcxmtx=(fgwithfreqcx && (g == 0) )?&fcxmtx:NULL;
ReadVisMtx(b,g,sernum,vmtx,pmmtx,pfcxmtx);
int npaq = vmtx.Info()["NPAQSUM"];
string sdate=vmtx.Info()["DATEOBS"];
dateobs=TimeStamp(sdate);
if (prtlevel>1) cout << " DateObs="<<dateobs<<endl;
if (prtlevel>2) cout << " ... vmtx: "<<vmtx.InfoString()<<endl;
if (prtlevel>3) cout << " ... vmtx.Info(): \n"<<vmtx.Info()<<endl;
myvismtx.SubMatrix(Range(rvf,rvl), Range(cvf, cvl)) = vmtx;
myvismtx.Info()["NPAQSUM"] = npaq;
if (prtlevel>3) cout << " ... vmtx.Info(): \n"<<vmtx.Info()<<endl;
// Dividing by the number of paquets used to compute the visibility.
complex<r_4> zfacnorm((r_4)npaq, 0.);
vmtx /= zfacnorm;
myvismtx.SubMatrix(Range(rvf,rvl), Range(cvf, cvl)) = vmtx;
if ((g==0)&&(b==0)) myvismtx.Info() = vmtx.Info();
char nmpaqsum[48];
sprintf(nmpaqsum,"NPAQSUM_%d_%d",b,g);
myvismtx.Info()[nmpaqsum] = npaq;
if(fgwithmean && (g == 0)) { //remplir la matrice des moyennes des coeff FFT
mmtx /= zfacnorm;
mmtxp->SubMatrix(Range::all(), Range(cvf, cvl)) = mmtx;
mmtxp->Info() = mmtx.Info();
}
if(fgwithfreqcx && (g == 0)) { //remplir la matrice des moyennes des coeff FFT
fcxmtx /= zfacnorm;
fcxmtxp->SubMatrix(Range::all(), Range(cvf, cvl)) = fcxmtx;
npaq = mmtx.Info()["NPAQSUM"];
mmtxp->Info()["NPAQSUM"] = npaq;
fcxmtxp->Info() = fcxmtx.Info();
}
mttag=vmtx.Info()["MeanTT"]; mttag/=(125.e3);
vmtt.push_back(mttag);
......@@ -190,6 +213,8 @@ bool VisiP4ReaderBase::Read_P(long sernum, TMatrix< complex<r_4> >* vismtx, TMat
if (fgreorderfreq_){
ReorderFreqs(myvismtx, *vismtx);
if(fgwithmean) ReorderFreqs(*mmtxp, *meanmtx);
if(fgwithfreqcx) ReorderFreqs(*fcxmtxp, *freqcxmtx);
}
return true;
......@@ -259,16 +284,23 @@ VisiP4ReaderNoDB::VisiP4ReaderNoDB(vector<string> const& vpath, int first, int l
}
/* --Methode-- */
void VisiP4ReaderNoDB::ReadVisMtx(int band, int visgrp, long sernum, TMatrix< complex<r_4> > & vmtx, TMatrix< complex<r_4> > * pmmtx)
void VisiP4ReaderNoDB::ReadVisMtx(int band, int visgrp, long sernum, TMatrix< complex<r_4> > & vmtx, TMatrix< complex<r_4> > * pmmtx, TMatrix< complex<r_4> > * fxmtx)
{
char filenamebuff[2048];
sprintf(filenamebuff,"%s/vismtx_%d_%ld.ppf",paths_[band].c_str(), visgrp, sernum);
if (prtlevel>3)
cout << " VisiP4ReaderNoDB::ReadVisMtx() - opening file: "<<filenamebuff<<endl;
PInPersist pin(filenamebuff);
if( pmmtx ) {
if( pmmtx && fxmtx ) {
pin>>vmtx >> (*pmmtx) >> (*fxmtx);
}
else if( pmmtx ) {
pin>>vmtx >> (*pmmtx);
} else {
}
else if( fxmtx ) {
pin>>vmtx >> (*fxmtx);
}
else {
pin>>vmtx;
}
return;
......@@ -383,7 +415,7 @@ void VisiP4ReaderDB::Init()
}
/* --Methode-- */
void VisiP4ReaderDB::ReadVisMtx(int band, int visgrp, long sernum, TMatrix< complex<r_4> > & vmtx, TMatrix< complex<r_4> > * mmtxp)
void VisiP4ReaderDB::ReadVisMtx(int band, int visgrp, long sernum, TMatrix< complex<r_4> > & vmtx, TMatrix< complex<r_4> > * mmtxp, TMatrix< complex<r_4> > * fxmtx)
{
if ((sernum < vdb_[cur_vdb_index_].ser_start)||(sernum >= vdb_[cur_vdb_index_].ser_end)) {
if ((sernum < vdb_[cur_vdb_index_+1].ser_start)||(sernum >= vdb_[cur_vdb_index_+1].ser_end)) {
......@@ -403,12 +435,17 @@ void VisiP4ReaderDB::ReadVisMtx(int band, int visgrp, long sernum, TMatrix< comp
if (!fgopenok_) OpenFiles(cur_vdb_index_);
char tagbuff[128];
// -- reading the visibility matrix
sprintf(tagbuff,"vismtx_%d_%ld",visgrp,sernum);
(*cur_filep_[band]) >> PPFNameTag(tagbuff) >> vmtx;
if (mmtxp) {
if (mmtxp) { // reading the mean fourier coeff matrix if requested
sprintf(tagbuff,"meanfourcoeff_%d_%ld",visgrp,sernum);
(*cur_filep_[band]) >> PPFNameTag(tagbuff) >> (*mmtxp);
}
if (fxmtx) { // reading the frequency cross-corr fo V_ii if requested
sprintf(tagbuff,"cxfreq_%d_%ld",visgrp,sernum);
(*cur_filep_[band]) >> PPFNameTag(tagbuff) >> (*mmtxp);
}
return;
}
......
......@@ -59,16 +59,25 @@ public:
fills also the corresponding TimeStamp (TU) and mean hardware time-tag. return true if OK, false otherwise. */
bool ReadNext(TMatrix< complex<r_4> > & vismtx, TimeStamp& dateobs, double& mttag)
{
bool readok = Read_P(cur_serialnum_, &vismtx, NULL, dateobs, mttag);
bool readok = Read_P(cur_serialnum_, &vismtx, NULL, NULL, dateobs, mttag);
if (readok) cur_serialnum_+=ser_step_;
return readok;
}
/*! \brief Read the next visibility matrix and mean Fourier coeff matrix in the selected visi matrix number range.
fills also the corresponding TimeStamp (TU) and mean hardware time-tag. return true if OK, false otherwise. */
bool ReadNext(TMatrix< complex<r_4> > & vismtx, TMatrix< complex<r_4> > & meanmtx,
TimeStamp& dateobs, double& mttag)
bool ReadNextMFFT(TMatrix< complex<r_4> > & vismtx, TMatrix< complex<r_4> > & meanmtx,
TimeStamp& dateobs, double& mttag)
{
bool readok = Read_P(cur_serialnum_, &vismtx, &meanmtx, dateobs, mttag);
bool readok = Read_P(cur_serialnum_, &vismtx, &meanmtx, NULL, dateobs, mttag);
if (readok) cur_serialnum_+=ser_step_;
return readok;
}
/*! \brief Read the next visibility matrix and cross-frequency correlation matrix in the selected visi matrix number range.
fills also the corresponding TimeStamp (TU) and mean hardware time-tag. return true if OK, false otherwise. */
bool ReadNextCxFreq(TMatrix< complex<r_4> > & vismtx, TMatrix< complex<r_4> > & freqcx,
TimeStamp& dateobs, double& mttag)
{
bool readok = Read_P(cur_serialnum_, &vismtx, NULL, &freqcx, dateobs, mttag);
if (readok) cur_serialnum_+=ser_step_;
return readok;
}
......@@ -79,16 +88,26 @@ public:
\warning : sernum should be in the selected range */
bool Read(long sernum, TMatrix< complex<r_4> > & vismtx, TimeStamp& dateobs, double& mttag)
{
return Read_P(sernum, &vismtx, NULL, dateobs, mttag);
return Read_P(sernum, &vismtx, NULL, NULL, dateobs, mttag);
}
/*! \brief Read the next visibility matrix and mean Fourier coeff matrix identified by the serial number \b sernum.
fills also the corresponding TimeStamp (TU) and mean hardware time-tag. return true if OK, false otherwise.
\warning : sernum should be in the selected range */
bool Read(long sernum, TMatrix< complex<r_4> > & vismtx, TMatrix< complex<r_4> > & meanmtx,
bool ReadMFFT(long sernum, TMatrix< complex<r_4> > & vismtx, TMatrix< complex<r_4> > & meanmtx,
TimeStamp& dateobs, double& mttag)
{
return Read_P(sernum, &vismtx, &meanmtx, dateobs, mttag);
return Read_P(sernum, &vismtx, &meanmtx, NULL, dateobs, mttag);
}
/*! \brief Read the next visibility matrix and cross frequency correlation matrix identified by the serial number \b sernum.
fills also the corresponding TimeStamp (TU) and mean hardware time-tag. return true if OK, false otherwise.
\warning : sernum should be in the selected range */
bool ReadCxFreq(long sernum, TMatrix< complex<r_4> > & vismtx, TMatrix< complex<r_4> > & freqcx,
TimeStamp& dateobs, double& mttag)
{
return Read_P(sernum, &vismtx, NULL, &freqcx, dateobs, mttag);
}
//! perform the reordering of the different frequency components (for FFT firmware of the ADC board)
......@@ -127,10 +146,12 @@ public:
protected:
bool Read_P(long sernum, TMatrix< complex<r_4> >* vismtx, TMatrix< complex<r_4> >* meanmtx,
// vismtx : matrice de visibilite , meanmtx: pointeur matrice moyenne coeff de Fourier , freqcxmtx: pointeur matrice correlation en deux freq pour les V_ii
bool Read_P(long sernum, TMatrix< complex<r_4> >* vismtx, TMatrix< complex<r_4> >* meanmtx, TMatrix< complex<r_4> >* freqcxmtx,
TimeStamp& dateobs, double& mttag);
virtual void ReadVisMtx(int band, int visgrp, long sernum, TMatrix< complex<r_4> > & vmtx, TMatrix< complex<r_4> > * mmtx) = 0;
// vmtx : matrice de visibilite , mmtx pointeur matrice moyenne coeff de Fourier, fxmtx : pointeur matrice correlation entre deux freq pour les V_ii
virtual void ReadVisMtx(int band, int visgrp, long sernum, TMatrix< complex<r_4> > & vmtx, TMatrix< complex<r_4> > * mmtx,
TMatrix< complex<r_4> > * fxmtx) = 0;
vector<string> paths_; // path for visibility files produced for each BAO5,BAO6 (first,second frequency band)
long ser_first_, ser_last_, ser_step_; // VisMtx serial number I such as ser_first_<=I<=ser_last_ , one matrix out of ser_step_
......@@ -158,7 +179,8 @@ public:
VisiP4ReaderNoDB(vector<string> const& vpath, int first, int last, int step=1, bool reorderfreq=false);
protected:
virtual void ReadVisMtx(int band, int visgrp, long sernum, TMatrix< complex<r_4> > & vmtx, TMatrix< complex<r_4> > * mmtx);
virtual void ReadVisMtx(int band, int visgrp, long sernum, TMatrix< complex<r_4> > & vmtx, TMatrix< complex<r_4> > * mmtx,
TMatrix< complex<r_4> > * fxmtx);
};
//------------------------------------------------------------------------------------------------
......@@ -180,7 +202,8 @@ public:
protected:
void Init();
virtual void ReadVisMtx(int band, int visgrp, long sernum, TMatrix< complex<r_4> > & vmtx, TMatrix< complex<r_4> > * mmtx);
virtual void ReadVisMtx(int band, int visgrp, long sernum, TMatrix< complex<r_4> > & vmtx, TMatrix< complex<r_4> > * mmtx,
TMatrix< complex<r_4> > * fxmtx);
void OpenFiles(size_t vdbidx);
long db_first_ser_, db_end_ser_; // first and last Matrix sequence number in directory
......
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