Changeset 3593 in Sophya for trunk


Ignore:
Timestamp:
Apr 7, 2009, 3:02:12 PM (16 years ago)
Author:
ansari
Message:

Ajout traitement fichiers de donnees FFT ds mfits2psec.cc , Reza 07/04/2009

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/AddOn/TAcq/mfits2spec.cc

    r3592 r3593  
    55/* ------------------------------------------------------------------
    66   Programme de calcul de spectre moyenne a partir des fichiers fits 
    7    d'acquisition de BAORadio - donnees brutes (non FFT)
     7   d'acquisition de BAORadio - donnees brutes ou fft
    88   R. Ansari, C. Magneville
    99   V1 : Juillet 2008  , V2 : Avril 2009
    … …  
    3737#include "brpaqu.h"
    3838
     39
     40// Definition de type pour les tableaux de float
     41#define DTF  r_4
     42
    3943//---- Declaration des fonctions de calcul ----
    4044// Fonction d'analyse 1ere version, pas d'entete ds le fichier, 1 voie
    … …  
    4448// Fonction d'analyse 2eme version , 2 voies / paquet
    4549int ana_data_2(vector<string>& infiles, string& oufile);
     50// Fonction d'analyse 2eme version , 2 voies / paquet
     51int ana_data_fft_2(vector<string>& infiles, string& oufile);
    4652
    4753//----------------------------------------------------
    … …  
    5056{
    5157  if (narg < 4) {
    52     cout << " ---Calcul spectres moyennes a partir de fits BAORadio " << endl;
     58    cout << " ---Mean spectrum computation from BAORadio FITS files" << endl;
    5359    cout << " Usage: mfits2spec ACT OutPPF file1 [file2 file3 ...] " << endl;
    54     cout << " ACT=-0,-1,-2 ==> 0: Nancay-Juil2008, - 1,2 : 1/2 voies / paquet " << endl;
    55     cout << " OutPPF : Output PPF file name, file1,file2 ... Input FITS files " << endl; 
     60    cout << "  Or :  mfits2spec ACT OutPPF -infits DirName Imin Imax " << endl;
     61    cout << " ACT=-0,-1,-2,-fft2 \n"
     62         << "   -0: Nancay-July2008 \n"
     63         << "   -1,-2 :  Raw data  1 ou 2 channels / packet(or frame) \n"
     64         << "   -fft2 :  FFT data 2 channels / packet " << endl;
     65    cout << " OutPPF : Output PPF file name " << endl; 
     66    cout << " file1,file2 ... : Input FITS files " << endl; 
     67    cout << " -infiles DirName Imin Imax : Input fits directory and sequence numbers \n"
     68         << "    FileNames=DirName/signalII.fits Imin<=II<=Imax " << endl; 
    5669    return 1;
    5770  }
    … …  
    6477    string outppf = arg[2];
    6578    vector<string> infiles;
    66     for(int i=3; i<narg; i++) infiles.push_back(arg[i]);
     79    if (strcmp(arg[3],"-infits")==0) {
     80      if (narg<7)  {
     81        cout << " mfits2spec/Error arguments - type mfits2spec for help" << endl;
     82        return 2;
     83      }
     84      char nbuff[1024];
     85      char* dirname = arg[4];
     86      int imin = atoi(arg[5]);
     87      int imax = atoi(arg[6]);
     88      for(int ii=imin; ii<=imax; ii++) {
     89        sprintf(nbuff,"%s/signal%d.fits",dirname,ii);
     90        infiles.push_back(nbuff);
     91      }
     92    }
     93    else // Liste de noms de fichiers fournis directement
     94      for(int i=3; i<narg; i++) infiles.push_back(arg[i]);
     95
    6796    cout << " ---------- mfits2spec.cc Start - ACT= " << act << " ------------- " << endl;
    6897    ResourceUsage resu;
    … …  
    7099    else if (act == "-1") ana_data_1(infiles, outppf);
    71100    else if (act == "-2") ana_data_2(infiles, outppf);
     101    else if (act == "-fft2") ana_data_fft_2(infiles, outppf);
    72102    else cout << " mfits2spec.cc / Bad argument ACT=" << act << " -> exit" << endl;
    73103    cout << resu ;
    … …  
    93123}
    94124
    95 inline r_4 Zmod2(complex<r_4> z)
     125inline DTF Zmod2(complex<DTF> z)
    96126{ return (z.real()*z.real()+z.imag()*z.imag()); }
    97127
    98128/*--Nouvelle-Fonction--*/
    99129int ana_data_0(vector<string>& infiles, string& outfile)
     130// Calcul spectre moyen 1 voie, donnees brutes, Donnees Juillet 2008
    100131{
    101132  TVector<float> spectre;
    … …  
    115146
    116147    Byte* data = new Byte[sx];
    117     TVector<r_4> vx(sx);
    118     TVector< complex<r_4> > cfour;
     148    TVector<DTF> vx(sx);
     149    TVector< complex<DTF> > cfour;
    119150    FFTPackServer ffts;
    120151     
    … …  
    122153      mff.ReadB(data, sx, j*sx);
    123154      // On convertit en float
    124       for(sa_size_t ix=0; ix<vx.Size(); ix++)  vx(ix) = (r_4)(data[ix]);
     155      for(sa_size_t ix=0; ix<vx.Size(); ix++)  vx(ix) = (DTF)(data[ix]);
    125156      // On fait le FFT
    126157      ffts.FFTForward(vx, cfour);
    … …  
    139170  }   
    140171  cout << "---- Ecriture spectre moyenne ds " << outfile << endl;     
    141   spectre /= (r_4)(nzm);
     172  spectre /= (DTF)(nzm);
    142173  POutPersist po(outfile);
    143174  po << spectre;
    … …  
    147178/*--Nouvelle-Fonction--*/
    148179int ana_data_1(vector<string>& infiles, string& outfile)
     180// Calcul spectre moyen 1 voie, donnees brutes 
    149181{
    150182  TVector<float> spectre;
    151183  float freq0 = 0; 
    152184  int paqsz = 0;
    153   int nfileok;
     185  int nfileok=0;
    154186  sa_size_t nzm = 0;   // Nb de spectres moyennes
    155187  Byte* data = NULL;
    156188  FFTPackServer ffts;
     189
     190  Timer tm("ana_data_1");
     191  size_t totnbytesrd = 0;
     192
    157193  for(int ifile=0; ifile<infiles.size(); ifile++) {
    158194    string ffname = infiles[ifile];
    … …  
    181217    size_t sy = mff.NAxis2();
    182218    BRPaquet paq(NULL, data, paqsz);
    183     TVector<r_4> vx(paq.DataSize());
    184     TVector< complex<r_4> > cfour;
     219    TVector<DTF> vx(paq.DataSize());
     220    TVector< complex<DTF> > cfour;
    185221     
    186222    for(int j=0; j<sy; j++) {
    187223      mff.ReadB(data, sx, j*sx);
    188224      // On convertit en float
    189       for(sa_size_t ix=0; ix<vx.Size(); ix++)  vx(ix) = (r_4)(paq.Data1()[ix])-127.5;
     225      for(sa_size_t ix=0; ix<vx.Size(); ix++)  vx(ix) = (DTF)(paq.Data1()[ix])-127.5;
    190226      // On fait le FFT
    191227      ffts.FFTForward(vx, cfour);
    … …  
    197233      nzm++;  freq0 += cfour(0).real();
    198234      }
    199     nfileok++;
    200     cout << "---- FIN lecture " << ffname << " NFileOK=" << nfileok << endl;     
     235    nfileok++;   
     236    size_t nbytesrd = sx*sy;
     237    totnbytesrd += nbytesrd;
     238    tm.Split();
     239    cout << "---- FIN lecture " << ffname << " NFileOK=" << nfileok << endl;
     240    cout << "  TotalDiskRead" << totnbytesrd/(1024*1024) << " MBytes Disk-Read rate= "
     241         << (double)(nbytesrd)/(1024.*1024.)/tm.PartialCPUTime() << " MB/s" << endl;   
    201242  }   
    202243  if (data) delete[] data;
    … …  
    206247  }
    207248  cout << "---- Ecriture spectre moyenne ds " << outfile  << endl;     
    208   spectre /= (r_4)(nzm);
     249  spectre /= (DTF)(nzm);
    209250  spectre.Info().Comment() = " SpectreMoyenne (Moyenne module^2) ";
    210251  spectre.Info()["NMOY"] = nzm; // Nombre de spectres moyennes
    … …  
    217258/*--Nouvelle-Fonction--*/
    218259int ana_data_2(vector<string>& infiles, string& outfile)
     260// Calcul spectres moyens 2 voies, donnees brutes 
    219261{
    220262  TVector<float> specV1, specV2;
    … …  
    223265  float freq0v2 = 0; 
    224266  int paqsz = 0;
    225   int nfileok;
     267  int nfileok=0;
    226268  sa_size_t nzm = 0;   // Nb de spectres moyennes
    227269  Byte* data = NULL;
    228270  FFTPackServer ffts;
     271
     272  Timer tm("ana_data_2");
     273  size_t totnbytesrd = 0;
     274
    229275  for(int ifile=0; ifile<infiles.size(); ifile++) {
    230276    string ffname = infiles[ifile];
    … …  
    254300    size_t sy = mff.NAxis2();
    255301    BRPaquet paq(NULL, data, paqsz);
    256     TVector<r_4> vx1(paq.DataSize()/2);
    257     TVector<r_4> vx2(paq.DataSize()/2);
    258     TVector< complex<r_4> > cfour1,cfour2;
     302    TVector<DTF> vx1(paq.DataSize()/2);
     303    TVector<DTF> vx2(paq.DataSize()/2);
     304    TVector< complex<DTF> > cfour1,cfour2;
    259305     
    260306    for(int j=0; j<sy; j++) {
    … …  
    262308     // if (j%50 == 0) cout << " *DBG* mff.ReadB() , j=" << j << endl;
    263309     // On convertit en float
    264       for(sa_size_t ix=0; ix<vx1.Size(); ix++)  vx1(ix) = (r_4)(paq.Data1()[ix])-127.5;
    265       for(sa_size_t ix=0; ix<vx2.Size(); ix++)  vx2(ix) = (r_4)(paq.Data2()[ix])-127.5;
     310      for(sa_size_t ix=0; ix<vx1.Size(); ix++)  vx1(ix) = (DTF)(paq.Data1()[ix])-127.5;
     311      for(sa_size_t ix=0; ix<vx2.Size(); ix++)  vx2(ix) = (DTF)(paq.Data2()[ix])-127.5;
    266312      // On fait le FFT
    267313      ffts.FFTForward(vx1, cfour1);
    … …  
    280326      }
    281327    nfileok++;
    282     cout << "---- FIN lecture " << ffname << " NFileOK=" << nfileok << endl;     
     328    size_t nbytesrd = sx*sy;
     329    totnbytesrd += nbytesrd;
     330    tm.Split();
     331    cout << "---- FIN lecture " << ffname << " NFileOK=" << nfileok << endl;
     332    cout << "  TotalDiskRead" << totnbytesrd/(1024*1024) << " MBytes Disk-Read rate= "
     333         << (double)(nbytesrd)/(1024.*1024.)/tm.PartialCPUTime() << " MB/s" << endl;   
    283334  } 
    284335  if (data) delete[] data;
    … …  
    288339  }
    289340  cout << "---- Ecriture spectre moyenne ds " << outfile  << endl;     
    290   specV1 /= (r_4)(nzm);
    291   specV2 /= (r_4)(nzm);
    292   cxspecV12 /= complex<r_4>((r_4)(nzm),0.);
    293   specV1.Info().Comment() = " SpectreMoyenne (Moyenne module^2) - Voie 1 (/2)";
    294   specV2.Info().Comment() = " SpectreMoyenne (Moyenne module^2) - Voie 2 (/2)";
     341  specV1 /= (DTF)(nzm);
     342  specV2 /= (DTF)(nzm);
     343  cxspecV12 /= complex<DTF>((DTF)(nzm),0.);
     344  specV1.Info().Comment() = " SpectreMoyenne (Moyenne module^2) - Raw-Voie 1 (/2)";
     345  specV2.Info().Comment() = " SpectreMoyenne (Moyenne module^2) - Raw-Voie 2 (/2)";
    295346  specV1.Info()["NMOY"] = specV2.Info()["NMOY"] = nzm; // Nombre de spectres moyennes
    296347  specV1.Info()["Freq0"] = freq0v1;
    … …  
    302353  return 0;
    303354}
     355
     356
     357inline DTF Zmod2TwoByte(TwoByteComplex z)
     358{ return (z.realD()*z.realD()+z.imagD()*z.imagD()); }
     359
     360/*--Nouvelle-Fonction--*/
     361int ana_data_fft_2(vector<string>& infiles, string& outfile)
     362// Calcul spectres moyens 2 voies, donnees brutes 
     363{
     364  TVector<float> specV1, specV2;
     365  TVector< complex<float> > cxspecV12;
     366  float freq0v1 = 0; 
     367  float freq0v2 = 0; 
     368  int paqsz = 0;
     369  int nfileok=0;
     370  sa_size_t nzm = 0;   // Nb de spectres moyennes
     371  Byte* data = NULL;
     372
     373  Timer tm("ana_data_fft_2");
     374  size_t totnbytesrd = 0;
     375
     376  for(int ifile=0; ifile<infiles.size(); ifile++) {
     377    string ffname = infiles[ifile];
     378// -------------- Lecture de bytes     
     379    cout << "ana_data_fft_2[" << ifile <<"]Ouverture/lecture fichier " << ffname << endl;
     380    MiniFITSFile mff(ffname, MF_Read);
     381    cout << "... Type=" << mff.DataTypeToString() << " NAxis1=" << mff.NAxis1()
     382         << " NAxis2=" << mff.NAxis2() << endl;
     383    if (mff.DataType() != MF_Byte) {
     384      cout << " PB : DataType!=MF_Byte --> skipping " << endl;
     385      continue;
     386    }
     387// Les fichier FITS contiennent l'entet (24 bytes), mais pas le trailer (16 bytes) ...
     388    if (paqsz == 0)  {  // premier passage, on fixe la taille de paquet et on alloue le buffer
     389      paqsz = mff.NAxis1()+16;
     390      cout << " ana_data_fft_2/ Allocating data , PaqSz=" << paqsz << endl;
     391      data = new Byte[paqsz];
     392      for(int ib=0; ib<paqsz; ib++) data[ib]=0;
     393    }
     394    else {
     395      if (paqsz != mff.NAxis1()+16) {
     396      cout << " PB : paqsz=" << paqsz << " != mff.NAxis1()+16 --> skipping " << endl;
     397      continue;
     398      }
     399    }
     400    size_t sx = mff.NAxis1();
     401    size_t sy = mff.NAxis2();
     402    BRPaquet paq(NULL, data, paqsz);
     403    TVector<DTF> vx1(paq.DataSize()/2);
     404    TVector<DTF> vx2(paq.DataSize()/2);
     405    TVector< complex<DTF> > cfour1,cfour2;
     406
     407    if (!specV1.IsAllocated())  specV1.SetSize(paq.DataSize()/4+1);
     408    if (!specV2.IsAllocated())  specV2.SetSize(paq.DataSize()/4+1);
     409    if (!cxspecV12.IsAllocated())  cxspecV12.SetSize(paq.DataSize()/4+1);
     410     
     411    for(int j=0; j<sy; j++) {
     412      mff.ReadB(data, sx, j*sx);
     413
     414      TwoByteComplex* zz1 = (TwoByteComplex*)paq.Data1();
     415      TwoByteComplex* zz2 = (TwoByteComplex*)paq.Data2();
     416      // Composante continue, partie reelle uniquement
     417      specV1(0) += zz1[0].realD()*zz1[0].realD(); 
     418      specV2(0) += zz2[0].realD()*zz2[0].realD(); 
     419      cxspecV12(0) += zz1[0].realD()*zz2[0].realD();
     420
     421      for(sa_size_t jf=1; jf<specV1.Size()-1; jf++)  {
     422        specV1(jf) += Zmod2TwoByte(zz1[jf]);
     423        specV2(jf) += Zmod2TwoByte(zz2[jf]);
     424        cxspecV12(jf) += (DTF)(zz1[jf].realD()*zz2[jf].realD()+zz1[jf].imagD()*zz2[jf].imagD());
     425      }
     426      // Freq. Nyquist a N/2 
     427      specV1(specV1.Size()-1) += zz1[0].imagD()*zz1[0].imagD();  // Freq. Nyquist a N/2
     428      specV2(specV2.Size()-1) += zz2[0].imagD()*zz2[0].imagD();  // Freq. Nyquist a N/2
     429      cxspecV12(cxspecV12.Size()-1) += zz1[0].imagD()*zz2[0].imagD();
     430
     431      nzm++;  freq0v1 += zz1[0].realD();  freq0v2 += zz2[0].realD();
     432      }
     433    nfileok++;
     434    size_t nbytesrd = sx*sy;
     435    totnbytesrd += nbytesrd;
     436    tm.Split();
     437    cout << "---- FIN lecture " << ffname << " NFileOK=" << nfileok << endl;
     438    cout << "  TotalDiskRead" << totnbytesrd/(1024*1024) << " MBytes Disk-Read rate= "
     439         << (double)(nbytesrd)/(1024.*1024.)/tm.PartialCPUTime() << " MB/s" << endl;   
     440  } 
     441  if (data) delete[] data;
     442  if (nzm <= 0) {
     443    cout << " ana_data_fft_2/ERROR : nzm=" << nzm << " spectres moyennes !" << endl;
     444    return nzm;
     445  }
     446  cout << "---- Ecriture spectre moyenne ds " << outfile  << endl;     
     447  specV1 /= (DTF)(nzm);
     448  specV2 /= (DTF)(nzm);
     449  cxspecV12 /= complex<DTF>((DTF)(nzm),0.);
     450  specV1.Info().Comment() = " SpectreMoyenne (Moyenne module^2) - FFT-Voie 1 (/2)";
     451  specV2.Info().Comment() = " SpectreMoyenne (Moyenne module^2) - FFT-Voie 2 (/2)";
     452  specV1.Info()["NMOY"] = specV2.Info()["NMOY"] = nzm; // Nombre de spectres moyennes
     453  specV1.Info()["Freq0"] = freq0v1;
     454  specV2.Info()["Freq0"] = freq0v2;
     455  POutPersist po(outfile);
     456  po << PPFNameTag("specV1") << specV1;
     457  po << PPFNameTag("specV2") << specV2;
     458  po << PPFNameTag("cxspecV12") << cxspecV12;
     459  return 0;
     460}
Note: See TracChangeset for help on using the changeset viewer.