Commit a143e301 authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

(JEC/RZ) 4/7/16 subtract offset to visi matrix and normalize by NPASUM

parent cf642bc6
......@@ -113,14 +113,16 @@ int main(int narg, char* arg[])
VisiP4Reader vreader(paths, Imin,Imax,Istep,true);
vreader.setPrintLevel(prtlev);
bool fgok=true;
TMatrix< complex<r_4> > vismtx, meanmtx; //vismtx: visib matrice, meanmtx: mean FFT coef
TMatrix< complex<r_4> > vismtx_mean, meanmtx_mean;
//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;
TimeStamp dateobs;
double mttag;
int cnt=0;
while (fgok) {
if(fgrdmean) {
fgok=vreader.ReadNext(vismtx, meanmtx, dateobs, mttag);
vismtx_center = vreader.SubtractOffset(vismtx, meanmtx);
} else {
fgok=vreader.ReadNext(vismtx, dateobs, mttag);
}
......@@ -130,6 +132,8 @@ int main(int narg, char* arg[])
if(fgrdmean) {
if (cnt==0) meanmtx_mean = meanmtx;
else meanmtx_mean += meanmtx;
if (cnt==0) vismtx_center_mean = vismtx_center;
else vismtx_center_mean += vismtx_center;
}
cnt++;
if (cnt%nmod==0) cout<<"rdvisip4/Info file read count="<<cnt<<" / nfiles="<<nfiles<<endl;
......@@ -147,6 +151,10 @@ int main(int narg, char* arg[])
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;
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;
}
}
......
......@@ -98,14 +98,18 @@ bool VisiP4Reader::ReadNext_P(TMatrix< complex<r_4> >* vismtx, TMatrix< complex<
} else {
pin>>vmtx;
}
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(fgwithmean && (g == 0)) { //remplir la matrice des moyennes des coeff FFT
mmtxp->SubMatrix(Range::all(), Range(cvf, cvl)) = mmtx;
npaq = mmtx.Info()["NPAQSUM"];
mmtxp->Info()["NPAQSUM"] = npaq;
}
mttag=vmtx.Info()["MeanTT"]; mttag/=(125.e3);
vmtt.push_back(mttag);
......@@ -149,6 +153,32 @@ void VisiP4Reader::ReorderFreqs(TMatrix< complex<r_4> > & svismtx, TMatrix< comp
// Voir Remarque ci-dessus Z(N/2).real -> Z(0).image , cad real(fmax)
dvismtx.Column(0)=svismtx.Column(NFREQ-1);
dvismtx.Info() = svismtx.Info();
return;
}
/* --Methode-- */
TMatrix< complex<r_4> > VisiP4Reader::SubtractOffset(TMatrix< complex<r_4> > & vismtx, TMatrix< complex<r_4> > & meanmtx)
{
r_4 npaqV = vismtx.Info()["NPAQSUM"];
r_4 npaqM = meanmtx.Info()["NPAQSUM"];
r_4 npaqM2 = npaqM*npaqM;
TMatrix< complex<r_4> > resmtx(vismtx.NRows(), vismtx.NCols());
for(sa_size_t ir=0; ir<resmtx.NRows(); ir++){
sa_size_t iCh = chanpairs(ir,0);
sa_size_t jCh = chanpairs(ir,1);
for(sa_size_t ic=0; ic<resmtx.NCols(); ic++){
resmtx(ir,ic) = vismtx(ir,ic)/complex<r_4>(npaqV,0.)
- meanmtx(iCh,ic)*conj(meanmtx(jCh,ic))/complex<r_4>(npaqM2, 0.);
}//loop cols : freq
}//loop rows : visi
return resmtx;
}
......@@ -33,6 +33,8 @@ public:
void ReorderFreqs(TMatrix< complex<r_4> > & svismtx, TMatrix< complex<r_4> > & vismtx);
TMatrix< complex<r_4> > SubtractOffset(TMatrix< complex<r_4> > & vismtx, TMatrix< complex<r_4> > & meanmtx);
inline void setNChan(sa_size_t nchan=8)
{ NCHAN=nchan; }
......
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