Commit 5636ca33 authored by Réza Ansari's avatar Réza Ansari
Browse files

1/ Debut d'implementation de la classe P4gnuGain pour calcul et application

du gain g(nu) reponse en frequence pour PAON-4
2/ Adaptation du makefile et programme rdvisip4.cc
                          Reza 16/09/2017
parent 467afa37
......@@ -5,23 +5,42 @@ include $(SOPHYABASE)/include/sophyamake.inc
OBJ = ./Objs/
EXE = ./Objs/
# List of include files of this package, and .o files to handle dependencies
MYINCLISTHERE = p4autils.h visip4reader.h p4gnugain.h
MYOLISTHERE = $(OBJ)/p4autils.o $(OBJ)/visip4reader.o $(OBJ)/p4gnugain.o
# Define our target list
all : rdvisip4 visi2ntac visi2dtacx visi2tmfreq p4conv2fits msvis2dt
clean :
rm -f $(EXE)/rdvisip4 $(EXE)/visi2ntac $(EXE)/visi2dtacx $(EXE)/visi2tmfreq $(EXE)/p4conv2fits $(EXE)/msvis2dt
rm -f $(OBJ)/rdvisip4.o $(OBJ)/visi2ntac.o $(OBJ)/visi2dtacx.o $(OBJ)/visi2tmfreq.o $(OBJ)/p4conv2fits.o $(OBJ)/msvis2dt.o
rm -f $(OBJ)/p4autils.o $(OBJ)/visip4reader.o
rm -f $(MYOLISTHERE)
###############################################################
### Compilation de .o
######
$(OBJ)/p4autils.o : p4autils.cc $(MYINCLISTHERE)
$(CXXCOMPILE) -o $(OBJ)/p4autils.o p4autils.cc
######
$(OBJ)/visip4reader.o : visip4reader.cc $(MYINCLISTHERE)
$(CXXCOMPILE) -o $(OBJ)/visip4reader.o visip4reader.cc
######
$(OBJ)/p4gnugain.o : p4gnugain.cc $(MYINCLISTHERE)
$(CXXCOMPILE) -o $(OBJ)/p4gnugain.o p4gnugain.cc
###############################################################
###### Compilation et link des executables
## programme de remplissage de DataTable avec cross-correlation fct du temps
visi2dtacx : $(EXE)/visi2dtacx
echo '---visi2dtacx made'
$(EXE)/visi2dtacx : $(OBJ)/visi2dtacx.o $(OBJ)/p4autils.o $(OBJ)/visip4reader.o
$(CXXLINK) -o $(EXE)/visi2dtacx $(OBJ)/visi2dtacx.o $(OBJ)/p4autils.o $(OBJ)/visip4reader.o $(SOPHYAEXTSLBLIST)
$(EXE)/visi2dtacx : $(OBJ)/visi2dtacx.o $(MYOLISTHERE)
$(CXXLINK) -o $(EXE)/visi2dtacx $(OBJ)/visi2dtacx.o $(MYOLISTHERE) $(SOPHYAEXTSLBLIST)
$(OBJ)/visi2dtacx.o : visi2dtacx.cc visip4reader.h
$(OBJ)/visi2dtacx.o : visi2dtacx.cc $(MYINCLISTHERE)
$(CXXCOMPILE) -o $(OBJ)/visi2dtacx.o visi2dtacx.cc
######
......@@ -29,10 +48,10 @@ $(OBJ)/visi2dtacx.o : visi2dtacx.cc visip4reader.h
visi2ntac : $(EXE)/visi2ntac
echo '---visi2ntac made'
$(EXE)/visi2ntac : $(OBJ)/visi2ntac.o $(OBJ)/visip4reader.o
$(CXXLINK) -o $(EXE)/visi2ntac $(OBJ)/visi2ntac.o $(OBJ)/visip4reader.o $(SOPHYAEXTSLBLIST)
$(EXE)/visi2ntac : $(OBJ)/visi2ntac.o $(MYOLISTHERE)
$(CXXLINK) -o $(EXE)/visi2ntac $(OBJ)/visi2ntac.o $(MYOLISTHERE) $(SOPHYAEXTSLBLIST)
$(OBJ)/visi2ntac.o : visi2ntac.cc visip4reader.h
$(OBJ)/visi2ntac.o : visi2ntac.cc $(MYINCLISTHERE)
$(CXXCOMPILE) -o $(OBJ)/visi2ntac.o visi2ntac.cc
######
......@@ -40,10 +59,10 @@ $(OBJ)/visi2ntac.o : visi2ntac.cc visip4reader.h
visi2tmfreq : $(EXE)/visi2tmfreq
echo '---visi2tmfreq made'
$(EXE)/visi2tmfreq : $(OBJ)/visi2tmfreq.o $(OBJ)/p4autils.o $(OBJ)/visip4reader.o
$(CXXLINK) -o $(EXE)/visi2tmfreq $(OBJ)/visi2tmfreq.o $(OBJ)/p4autils.o $(OBJ)/visip4reader.o $(SOPHYAEXTSLBLIST)
$(EXE)/visi2tmfreq : $(OBJ)/visi2tmfreq.o $(MYOLISTHERE)
$(CXXLINK) -o $(EXE)/visi2tmfreq $(OBJ)/visi2tmfreq.o $(MYOLISTHERE) $(SOPHYAEXTSLBLIST)
$(OBJ)/visi2tmfreq.o : visi2tmfreq.cc visip4reader.h p4autils.h
$(OBJ)/visi2tmfreq.o : visi2tmfreq.cc $(MYINCLISTHERE)
$(CXXCOMPILE) -o $(OBJ)/visi2tmfreq.o visi2tmfreq.cc
######
......@@ -51,10 +70,10 @@ $(OBJ)/visi2tmfreq.o : visi2tmfreq.cc visip4reader.h p4autils.h
rdvisip4 : $(EXE)/rdvisip4
echo '---rdvisip4 made'
$(EXE)/rdvisip4 : $(OBJ)/rdvisip4.o $(OBJ)/p4autils.o $(OBJ)/visip4reader.o
$(CXXLINK) -o $(EXE)/rdvisip4 $(OBJ)/rdvisip4.o $(OBJ)/p4autils.o $(OBJ)/visip4reader.o $(SOPHYAEXTSLBLIST)
$(EXE)/rdvisip4 : $(OBJ)/rdvisip4.o $(MYOLISTHERE)
$(CXXLINK) -o $(EXE)/rdvisip4 $(OBJ)/rdvisip4.o $(MYOLISTHERE) $(SOPHYAEXTSLBLIST)
$(OBJ)/rdvisip4.o : rdvisip4.cc visip4reader.h
$(OBJ)/rdvisip4.o : rdvisip4.cc $(MYINCLISTHERE)
$(CXXCOMPILE) -o $(OBJ)/rdvisip4.o rdvisip4.cc
######
......@@ -80,10 +99,3 @@ $(EXE)/msvis2dt : $(OBJ)/msvis2dt.o $(OBJ)/msvis2dt.o
$(OBJ)/msvis2dt.o : msvis2dt.cc
$(CXXCOMPILE) -o $(OBJ)/msvis2dt.o msvis2dt.cc
######
$(OBJ)/p4autils.o : p4autils.cc p4autils.h
$(CXXCOMPILE) -o $(OBJ)/p4autils.o p4autils.cc
######
$(OBJ)/visip4reader.o : visip4reader.cc visip4reader.h
$(CXXCOMPILE) -o $(OBJ)/visip4reader.o visip4reader.cc
/* ----------------------------------------------------------
Projet BAORadio/PAON4 - (C) LAL/IRFU 2008-2017
Classes et fonctions utilitaires pour les programmes d'analyse PAON4
R. Ansari, C. Magneville Automne 2017
---------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <string>
#include <string>
#include "pexceptions.h"
#include "array.h"
#include "fioarr.h"
#include "p4gnugain.h"
#include "p4autils.h"
#include "nbtrixx.h"
#include "cspline.h"
using namespace SOPHYA ;
using namespace std;
/*--Fonction--*/
void P4gnuGain::computeSaveGain(SOPHYA::TMatrix< complex<r_4> > & vismtx_mean, std::string const & gainfile)
{
sa_size_t KVAC_CK[8]={0,8,15,21,26,30,33,35};
vector<sa_size_t> KVAC = P4AVisiNumEncoder::getAllAutoCor();
if (KVAC.size() != 8) throw PError("P4gnuGain::computeSaveGain() KVAC.size() != 8");
for(size_t i=0; i<KVAC.size(); i++) {
if (KVAC[i] != KVAC_CK[i]) {
cout << "P4gnuGain::computeSaveGain()/BUG i="<<i<<" KVAC[i]="<<KVAC[i]<<" KVAC_CK[i]="<<KVAC_CK[i]<<endl;
throw PError("P4gnuGain::computeSaveGain() KVAC.size() != 8");
}
}
Matrix acs(8, vismtx_mean.NCols());
for(sa_size_t r=0; r<8; r++)
for(sa_size_t c=0; c<acs.NCols(); c++) acs(r,c)=vismtx_mean(KVAC[r],c).real();
int SZW=96; float frackeep=0.75;
cout << "P4gnuGain::computeSaveGain/Info: computing gains with SZW="<<SZW<<" frackeep="<<frackeep<<" ..."<<endl;
Matrix gains = P4gnuGain::computeGain2(acs, SZW, frackeep);
Vector gn = P4gnuGain::normalizeGain2(gains);
cout << "P4gnuGain::computeSaveGain/Info: saving gains, gn to "<<gainfile<<endl;
POutPersist po(gainfile);
po<<PPFNameTag("gains")<<gains;
po<<PPFNameTag("gn")<<gn;
return;
}
// ----- adaptation de la fonction de calcul de gain, ecrit par cmv pour l'analyse Amas (Aout 2013)
/* --- fonction de calcul de gain avec fit de spline qui flag tout seul
les mauvais points du spectre - calcule la moyenne par fenetre (en
utilisant une mediane) et fit de spline sur les un point / fenetre
base sur la fonction rawgain2 de cmv (en date du gain (fonction de normalisation)
des spectres ON et OFF separes
*/
/*--Fonction--*/
SOPHYA::Matrix P4gnuGain::computeGain2(SOPHYA::Matrix & acs, int SZW, float frackeep)
{
if(SZW%2) SZW++;
int nrows = acs.NRows(), ncols = acs.NCols();
double *val = new double[2*SZW+2];
Matrix gain(nrows,ncols); gain = 0.;
for(int r=0;r<nrows;r++) {
TVector<r_8> DVon(ncols); DVon = 0.;
for(int c=1; c<ncols-1; c++) DVon(c) = fabs(acs(r,c) - 0.5*(acs(r,c-1)+acs(r,c+1)));
vector<double> vsplx, vsply;
for(int c=0; c<ncols; c+=SZW/4) {
// on enleve d'abord les RFI tres fortes
int ng = 0;
for(int cc=c;cc<c+SZW&&cc<ncols;cc++) {val[ng] = DVon(cc); ng++;}
if(ng<2) continue;
TabSort(ng,val);
double dvcut = val[int(ng*frackeep)];
//
ng = 0;
double x = 0.;
for(int cc=c;cc<c+SZW&&cc<ncols;cc++) if(DVon(cc)<dvcut) {val[ng] = acs(r,cc); x+=cc; ng++;}
if(ng!=0) x /= (double)ng;
////cout<<"rawgain2: r="<<r<<" c="<<c<<" dvcut="<<dvcut<<" ng="<<ng<<" x="<<x<<" v="<<val[ng/2]<<endl;
if(ng<2) continue;
TabSort(ng,val);
vsplx.push_back(x);
vsply.push_back(val[ng/2]);
}
TVector<r_8> vvx(vsplx), vvy(vsply);
CSpline csp(vvx.Size(),vvx.Data(),vvy.Data());
csp.ComputeCSpline();
for(int c=0;c<ncols;c++) gain(r,c) = csp.CSplineInt(c);
}
delete [] val;
return gain;
}
/*--Fonction--*/
SOPHYA::Vector P4gnuGain::normalizeGain2(SOPHYA::Matrix & gain)
{
char buff[16];
Vector gnorm(gain.NRows());
for(sa_size_t r=0;r<gain.NRows();r++) {
double gn=gain.Row(r).Sum();
double fcg=1000./gn; // we decide arbitray to renormalize the gain to have the sum[gains(i)]=1000
gain.Row(r) *= fcg;
sprintf(buff,"GNORM%d",(int)r);
gain.Info()[buff]=gn;
gnorm(r)=gn/1000.;
size_t badcnt=0;
#define LOWGTHR 1.e-9
for(sa_size_t c=0; c<gain.NCols(); c++)
if (gain(r,c)<LOWGTHR) {
badcnt++; gain(r,c)=LOWGTHR;
}
cout << " normalizeGain2()/Info GNorm["<<r<<"]= "<<gn<<" (badcount="<<badcnt<<")"<<endl;
}
return gnorm;
}
/* ----------------------------------------------------------
Projet BAORadio/PAON4 - (C) LAL/IRFU 2008-2017
Classes et fonctions utilitaires pour les programmes d'analyse PAON4
R. Ansari, C. Magneville Automne 2017
---------------------------------------------------------- */
#ifndef P4GNUGAINS_SEEN
#define P4GNUGAINS_SEEN
#include "array.h"
//---------------------------------------------------------------------------------
//---- Gestion des gains g(nu) - reponse en frequence pour PAON-4
//! \brief class for managing the frequency dependent gain g(nu) for PAON-4
class P4gnuGain {
public:
/*! \brief compute the g(nu) frequency response gains for the 8 PAON-4 signals and saves
the resulting
\b vismtx_mean is a visibility matrix averaged over few minutes, from which the gains g(nu)
is determined. The 8 g_i(nu) gains are saved as a matrix, as well as overall nomalisation
to the PPF file \b gainfile.
*/
static void computeSaveGain(SOPHYA::TMatrix< complex<r_4> > & vismtx_mean, std::string const & gainfile);
protected:
// ----- adaptation de la fonction de calcul de gain, ecrit par cmv pour l'analyse Amas (Aout 2013)
static SOPHYA::Matrix computeGain2(SOPHYA::Matrix & acs, int SZW=96, float frackeep=0.75);
// ----- fonction pour normaliser les gains
static SOPHYA::Vector normalizeGain2(SOPHYA::Matrix & gain);
};
#endif
......@@ -8,7 +8,7 @@
Programme de lecture des fichiers matrices de
visibilites de PAON4 produits mfacq
R. Ansari, J.E. Campagne, C. Magneville - LAL/Irfu
V0 : Fev 2015
V0 : Fev 2015 - Updated Sep 2017
---------------------------------------------------------- */
// include standard c/c++
......@@ -39,15 +39,8 @@
// include lecteur de fichiers visibilites
#include "p4autils.h"
#include "visip4reader.h"
#include "p4gnugain.h"
// ---- compute gain for each channel from visibilies and save it to file
void computeSaveGain(TMatrix< complex<r_4> > & vismtx_mean, string const & gainfile);
// ----- adaptation de la fonction de calcul de gain, ecrit par cmv pour l'analyse Amas (Aout 2013)
Matrix computeGain2(Matrix & acs, int SZW=96,float frackeep=0.75);
// ----- fonction pour normaliser les gains
Vector normalizeGain2(Matrix& gain);
int Usage(void);
int Usage(void)
{
cout<<"--- rdvisip4.cc : Read PPF files produced by mfacq time-frequency\n"<<endl;
......@@ -130,7 +123,7 @@ int main(int narg, const char* arg[])
if (cnt>0) {
vismtx_mean /= complex<r_4>((r_4)cnt, (r_4)0.);
cout<<" rdvisip4/Info: computing gains ..."<<endl;
computeSaveGain(vismtx_mean, gainfile);
P4gnuGain::computeSaveGain(vismtx_mean, gainfile);
cout<<" rdvisip4/Info: Saving vismtx_mean to PPF file "<<outfile<<endl;
POutPersist po(outfile);
po<<vismtx_mean;
......@@ -168,94 +161,3 @@ int main(int narg, const char* arg[])
}
//--------------------- Fonctions pour le calcul de gain -----------------------
//----- include specifique pour computeGain2()
#include "nbtrixx.h"
#include "cspline.h"
/*--Fonction--*/
void computeSaveGain(TMatrix< complex<r_4> > & vismtx_mean, string const & gainfile)
{
sa_size_t KVAC[8]={0,8,15,21,26,30,33,35};
Matrix acs(8, vismtx_mean.NCols());
for(sa_size_t r=0; r<8; r++)
for(sa_size_t c=0; c<acs.NCols(); c++) acs(r,c)=vismtx_mean(KVAC[r],c).real();
int SZW=96; float frackeep=0.75;
cout << " computeSaveGain/Info: computing gains with SZW="<<SZW<<" frackeep="<<frackeep<<" ..."<<endl;
Matrix gains = computeGain2(acs, SZW, frackeep);
Vector gn = normalizeGain2(gains);
cout << " computeSaveGain/Info: saving gains, gn to "<<gainfile<<endl;
POutPersist po(gainfile);
po<<PPFNameTag("gains")<<gains;
po<<PPFNameTag("gn")<<gn;
return;
}
// ----- adaptation de la fonction de calcul de gain, ecrit par cmv pour l'analyse Amas (Aout 2013)
/* --- fonction de calcul de gain avec fit de spline qui flag tout seul
les mauvais points du spectre - calcule la moyenne par fenetre (en
utilisant une mediane) et fit de spline sur les un point / fenetre
base sur la fonction rawgain2 de cmv (en date du gain (fonction de normalisation)
des spectres ON et OFF separes
*/
/*--Fonction--*/
Matrix computeGain2(Matrix& acs, int SZW,float frackeep)
{
if(SZW%2) SZW++;
int nrows = acs.NRows(), ncols = acs.NCols();
double *val = new double[2*SZW+2];
Matrix gain(nrows,ncols); gain = 0.;
for(int r=0;r<nrows;r++) {
TVector<r_8> DVon(ncols); DVon = 0.;
for(int c=1; c<ncols-1; c++) DVon(c) = fabs(acs(r,c) - 0.5*(acs(r,c-1)+acs(r,c+1)));
vector<double> vsplx, vsply;
for(int c=0; c<ncols; c+=SZW/4) {
// on enleve d'abord les RFI tres fortes
int ng = 0;
for(int cc=c;cc<c+SZW&&cc<ncols;cc++) {val[ng] = DVon(cc); ng++;}
if(ng<2) continue;
TabSort(ng,val);
double dvcut = val[int(ng*frackeep)];
//
ng = 0;
double x = 0.;
for(int cc=c;cc<c+SZW&&cc<ncols;cc++) if(DVon(cc)<dvcut) {val[ng] = acs(r,cc); x+=cc; ng++;}
if(ng!=0) x /= (double)ng;
////cout<<"rawgain2: r="<<r<<" c="<<c<<" dvcut="<<dvcut<<" ng="<<ng<<" x="<<x<<" v="<<val[ng/2]<<endl;
if(ng<2) continue;
TabSort(ng,val);
vsplx.push_back(x);
vsply.push_back(val[ng/2]);
}
TVector<r_8> vvx(vsplx), vvy(vsply);
CSpline csp(vvx.Size(),vvx.Data(),vvy.Data());
csp.ComputeCSpline();
for(int c=0;c<ncols;c++) gain(r,c) = csp.CSplineInt(c);
}
delete [] val;
return gain;
}
/*--Fonction--*/
Vector normalizeGain2(Matrix& gain)
{
char buff[16];
Vector gnorm(gain.NRows());
for(sa_size_t r=0;r<gain.NRows();r++) {
double gn=gain.Row(r).Sum();
double fcg=1000./gn; // we decide arbitray to renormalize the gain to have the sum[gains(i)]=1000
gain.Row(r) *= fcg;
sprintf(buff,"GNORM%d",(int)r);
gain.Info()[buff]=gn;
gnorm(r)=gn/1000.;
size_t badcnt=0;
#define LOWGTHR 1.e-9
for(sa_size_t c=0; c<gain.NCols(); c++)
if (gain(r,c)<LOWGTHR) {
badcnt++; gain(r,c)=LOWGTHR;
}
cout << " normalizeGain2()/Info GNorm["<<r<<"]= "<<gn<<" (badcount="<<badcnt<<")"<<endl;
}
return gnorm;
}
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