Skip to content

GitLab

Commit 6e5960df authored by perdereau's avatar perdereau
Browse files

select freq band in visiavg

parent 5ae93f0f
Hide whitespace changes
Inline Side-by-side
Showing with 79 additions and 36 deletions
p4autils.cc
View file @ 6e5960df
......@@ -357,11 +357,13 @@ int P4AnaParams::DecodeArgs(int narg, const char* arg[])
else if (fbo=="-freq") {
if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl; return -1; }
double f0,df;
cout << " selecting freqs "<<endl;
sscanf(arg[2],"%lg,%lg",&f0,&df); fbands_.push_back(P4FreqBand(f0,df)); fgdtable_=true;
arg+=2; narg-=2; lastargs_.clear();
}
else if (fbo=="-killfreq") {
if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl; return -1; }
cout << " killing freqs "<<endl;
double f0,df;
sscanf(arg[2],"%lg,%lg",&f0,&df); kill_fbands_.push_back(P4FreqBand(f0,df));
arg+=2; narg-=2; lastargs_.clear();
......@@ -598,8 +600,8 @@ ostream& P4AnaParams::Print(ostream& os) const
os<<" Right Ascension: center_="<<RAcenter_<<" Width="<<RAwidth_<<" (hours) Resolution="<<RAresol_mn_
<<"("<<RAcenter_-RAwidth_/2<<"<ra<"<<RAcenter_-RAwidth_/2<<")"<<endl;
if (fbands_.size()==0) os<<" No frequency band defined ..."<<endl;
else for(size_t i=0; i<fbands_.size(); i++) os<<"FreqBand["<<i<<"]: "<<fbands_[i]<<endl;
if (kill_fbands_.size()==0) {
else for(size_t i=0; i<fbands_.size(); i++) os<<"Selected FreqBand["<<i<<"]: "<<fbands_[i]<<endl;
if (kill_fbands_.size()>0) {
os<<" Frequency bands to mask/kill defined :"<<endl;
for(size_t i=0; i<fbands_.size(); i++) os<<"KillFreqBand["<<i<<"]: "<<fbands_[i]<<endl;
}
......
p4autils.h
View file @ 6e5960df
/* ----------------------------------------------------------
Projet BAORadio/PAON4 - (C) LAL/IRFU 2008-2017
Projet BAORadio/PAON4 - (C) LAL 2008-2017
Classes et fonctions utilitaires pour les programmes d'analyse PAON4
R. Ansari, C. Magneville Fevrier-Mars 2017
......
visiavg.cc
View file @ 6e5960df
......@@ -87,7 +87,8 @@ int main(int narg, const char* arg[])
cout <<"visiavg/Info: Path BAO5:"<<params.inpath5_<<" BAO6:"<<params.inpath6_<<"\n"
<<"fgreorderfreq="<<params.fgreorderfreq_<<"\n"
<<" DeltaIAvg="<<deltaIavg<<"\n"
<<"outfile="<<outfile<<" PrtLev="<<prtlev<<endl;
<< " fband # "<<params.fbands_.size()
<<" outfile="<<outfile<<" PrtLev="<<prtlev<<endl;
if (!FgTFM) {
cout<<" visiavg/parameter error : specify Time-Frequency map parameter with -tfm "<<endl;
......@@ -131,7 +132,7 @@ int main(int narg, const char* arg[])
bool fgok=true;
// vecteur de noms
vector <string> ext_names;
// un vecteur avec les temps
// un vecteur avec les temps & un avec RA
TVector< double > timevec(wreader.getTotalNbWindows()/deltaIavg);
TVector< double > ravec(wreader.getTotalNbWindows()/deltaIavg);
TMatrix< complex<r_4> > vismtx;
......@@ -190,6 +191,12 @@ int main(int narg, const char* arg[])
sa_size_t tfmSX, tfmSY;
long serialnum;
size_t rdidx, lastrdidx=0;
double freqmin = P4FreqBand::freqstart_;
double freqmax = P4FreqBand::freqend_;
sa_size_t ifrmin = 0 ;
sa_size_t ifrmax = P4FreqBand::nbnufft_ ;
while (fgok) {
//reads next visibility matrix window
fgok = wreader.Shift(rdidx, serialnum);
......@@ -203,39 +210,58 @@ int main(int narg, const char* arg[])
lastrdidx=rdidx;
vismtx = wreader.getAverageVisMtx(cfdate);
///cout << vismtx << endl;
sa_size_t NTFCols = vismtx.NCols(); // bin en freq a ecrire
if (cnt==0) { //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations
acsum.SetSize(8, vismtx.NCols());
tfmSX=wreader.getTotalNbWindows()/deltaIavg;
tfmSY=vismtx.NCols()/TFMfbin;
// default
if ( params.fbands_.size()>0 ){
sa_size_t XtfmSY = (sa_size_t) floor(params.fbands_[0].df_/params.fbands_[0].deltanufft_ /TFMfbin ) ;
cout << " taille finale en freq "<< XtfmSY << " bin "<< TFMfbin <<endl;
tfmSY=XtfmSY;
ifrmin = (sa_size_t) params.fbands_[0].jfmin_ ;
ifrmax = (sa_size_t) params.fbands_[0].jfmax_ ;
freqmin = params.fbands_[0].getP4Frequency(ifrmin);
freqmax = params.fbands_[0].getP4Frequency(ifrmax);
NTFCols = ifrmax - ifrmin +1 ; //## CHECK THIS
cout << " fmin,max "<< freqmin<<","<< freqmax<<" ident "<<ifrmin<<" "<< ifrmax << endl;
}
acsum.SetSize(8, NTFCols);
if (params.doSigma_) // for computing Sigma, if required
acsum_sq.SetSize(8, vismtx.NCols());
acsum_sq.SetSize(8, NTFCols);
cxsum.SetSize(6, vismtx.NCols());
cxsum.SetSize(6, NTFCols);
if (params.doSigma_) { // for computing Sigma, if required
cxsum_sq_rp.SetSize(6, vismtx.NCols());
cxsum_sq_ip.SetSize(6, vismtx.NCols());
cxsum_sq_rp.SetSize(6, NTFCols);
cxsum_sq_ip.SetSize(6, NTFCols);
}
// VV if needed
if (params.doVV_){
cxsum_vv.SetSize(6, vismtx.NCols());
cxsum_vv.SetSize(6, NTFCols);
if (params.doSigma_) { // for computing Sigma, if required
cxsum_vv_sq_rp.SetSize(6, vismtx.NCols());
cxsum_vv_sq_ip.SetSize(6, vismtx.NCols());
cxsum_vv_sq_rp.SetSize(6, NTFCols);
cxsum_vv_sq_ip.SetSize(6, NTFCols);
}
}
// HV if needed
if (params.doHV_){
cxsum_hv.SetSize(16, vismtx.NCols());
cxsum_hv.SetSize(16, NTFCols);
if (params.doSigma_) { // for computing Sigma, if required
cxsum_hv_sq_rp.SetSize(16, vismtx.NCols());
cxsum_hv_sq_ip.SetSize(16, vismtx.NCols());
cxsum_hv_sq_rp.SetSize(16, NTFCols);
cxsum_hv_sq_ip.SetSize(16, NTFCols);
}
}
tfmSX=wreader.getTotalNbWindows()/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) );
......@@ -312,31 +338,38 @@ int main(int narg, const char* arg[])
}// end if(I==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<KVAC.size(); k++){
//cout << vismtx.Row(KVAC[k]) << endl;
//cout << (vismtx.Row(KVAC[k])).SubMatrix(Range::all(),Range(ifrmin,ifrmax))<< endl;
//cout <<acsum.Row(k) << endl;
acsum.Row(k) += (vismtx.Row(KVAC[k])).SubMatrix( Range::all(),Range(ifrmin,ifrmax));
}
//cout << " end ac "<<endl;
if (params.doSigma_) { // for computing Sigma, if required
for(size_t k=0; k<KVAC.size(); k++) {
TVector<r_4> tmp = real(vismtx.Row(KVAC[k]));
TVector<r_4> tmp = real(vismtx.Row(KVAC[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax))) ; //(Range(ifrmin,ifrmax)));
acsum_sq.Row(k) += tmp.MulElt(tmp,tmp) ;
}
}
//cout << " end ac "<<endl;
for(size_t k=0; k<KVCXHH.size(); k++){
cxsum.Row(k) += vismtx.Row(KVCXHH[k]); // les cross-correlations
cxsum.Row(k) += vismtx.Row(KVCXHH[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)); // (Range(ifrmin,ifrmax)); // les cross-correlations
if (params.doSigma_) { // for computing Sigma, if required
TVector<r_4> tmp = real(vismtx.Row(KVCXHH[k]));
TVector<r_4> tmp = real(vismtx.Row(KVCXHH[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)));//(Range(ifrmin,ifrmax)));
cxsum_sq_rp.Row(k) += tmp.MulElt(tmp,tmp) ;
tmp = imag(vismtx.Row(KVCXHH[k]));
tmp = imag(vismtx.Row(KVCXHH[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)));//(Range(ifrmin,ifrmax)));
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
cxsum_vv.Row(k) += vismtx.Row(KVCXVV[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)) ;//(Range(ifrmin,ifrmax)); // les cross-correlations
if (params.doSigma_) { // for computing Sigma, if required
TVector<r_4> tmp = real(vismtx.Row(KVCXVV[k]));
TVector<r_4> tmp = real(vismtx.Row(KVCXVV[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)));//(Range(ifrmin,ifrmax)));
cxsum_vv_sq_rp.Row(k) += tmp.MulElt(tmp,tmp) ;
tmp = imag(vismtx.Row(KVCXVV[k]));
tmp = imag(vismtx.Row(KVCXVV[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)));//(Range(ifrmin,ifrmax)));
cxsum_vv_sq_ip.Row(k) += tmp.MulElt(tmp,tmp) ;
}
}
......@@ -344,17 +377,17 @@ int main(int narg, const char* arg[])
if (params.doHV_){
for(size_t k=0; k<KVCXHV.size(); k++){
cxsum_hv.Row(k) += vismtx.Row(KVCXHV[k]); // les cross-correlations HV
cxsum_hv.Row(k) += vismtx.Row(KVCXHV[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax));//(Range(ifrmin,ifrmax)); // les cross-correlations HV
if (params.doSigma_) { // for computing Sigma, if required
TVector<r_4> tmp = real(vismtx.Row(KVCXHV[k]));
TVector<r_4> tmp = real(vismtx.Row(KVCXHV[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)));//(Range(ifrmin,ifrmax)));
cxsum_hv_sq_rp.Row(k) += tmp.MulElt(tmp,tmp) ;
tmp = imag(vismtx.Row(KVCXHV[k]));
tmp = imag(vismtx.Row(KVCXHV[k]).SubMatrix(Range::all(),Range(ifrmin,ifrmax)));//(Range(ifrmin,ifrmax)));
cxsum_hv_sq_ip.Row(k) += tmp.MulElt(tmp,tmp) ;
}
}
}
//cout << " at i++ "<<endl;
I++; // incrementing DeltaTime counter
// we check that our time index did not go beyond the allocated array size (might not be necessary)
......@@ -649,15 +682,19 @@ int main(int narg, const char* arg[])
}
} // end option HV
P4FreqBand myp4fre;
if(params.fbands_.size()>0 )
myp4fre=params.fbands_[0];
TVector <double> lim_freq(2);
if (params.gain_gnu_file_.length()>0) {
P4gnuGain p4gnu( params.gain_gnu_file_ );
lim_freq(0) = p4gnu. minGoodF();
lim_freq(1) = p4gnu. maxGoodF();
}else{
lim_freq(0) = myp4fre.freqstart_;
lim_freq(1) = myp4fre.freqend_;
lim_freq(0) = myp4fre.getP4Frequency(params.fbands_[0].jfmin_);
lim_freq(1) = myp4fre.getP4Frequency(params.fbands_[0].jfmax_);
}
cout << " frequences limites "<< lim_freq(0) <<" ; "<< lim_freq(1)<<endl;
......@@ -668,11 +705,15 @@ int main(int narg, const char* arg[])
// resu.Update();
size_t szFreqs = myp4fre.getP4NbFreqChannels()/TFMfbin;
if (params.fbands_.size()>0 )
szFreqs = (params.fbands_[0].jfmax_ - params.fbands_[0].jfmin_)/TFMfbin;
TVector <double> avg_freqs( szFreqs );
TVector <double> avg_freqs( myp4fre.getP4NbFreqChannels()/TFMfbin);
double frbase = myp4fre.freqstart_ + myp4fre.getP4FreqResolution()/2. ;
for (int kf=0 ; kf< myp4fre.getP4NbFreqChannels()/TFMfbin ; kf++,frbase += myp4fre.getP4FreqResolution()*TFMfbin )
double frbase = myp4fre.getP4Frequency(params.fbands_[0].jfmin_) + myp4fre.getP4FreqResolution()/2. ;
for (int kf=0 ; kf< szFreqs ; kf++,frbase += myp4fre.getP4FreqResolution()*TFMfbin )
avg_freqs(kf) = frbase ;
potfm << PPFNameTag("FreqVec") << avg_freqs;
......
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