Changeset 3592 in Sophya for trunk/AddOn/TAcq


Ignore:
Timestamp:
Apr 6, 2009, 11:57:39 PM (16 years ago)
Author:
ansari
Message:

Amelioration programme mfits2pec.cc - Reza 06/04/2009

Location:
trunk/AddOn/TAcq
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/AddOn/TAcq/brpaqu.cc

    r3538 r3592  
    1111}
    1212
    13 BRPaquet::BRPaquet(Byte* src, Byte* dst, int paqsz)
     13/* --Methode__ */
     14BRPaquet::BRPaquet(Byte* src, Byte* dst, int paqsz, BRDataFmtConv fgswap)
     15  // swapall = true -> on swap tout le paquet, sinon swap entete seulement
    1416{
    1517  dst_ = dst;
    1618  sz_ = paqsz;
    17   if ((src==NULL)||(dst==NULL))  return;
    18   // On copie l'entete sans swapper
    19   for(int k=0; k<BRHDRSIZE; k++)  dst_[k] = src[k];
    20   // On copie la zone donnees en faisant un byte-swap correspondant a 4 octets
    21    for(int ka=BRHDRSIZE; ka<sz_-BRTRLSIZE; ka+=4) {
    22     for(int kb=0; kb<4; kb++)
    23       dst_[ka+kb] = src[ka+3-kb];
     19  if ((src == NULL) || (dst == NULL)) return;
     20  // Il faut mettre une protection (throw) si dst==NULL ou sz==0
     21
     22  UInt32* src32 = (UInt32*)src;
     23  UInt32* dst32 = (UInt32*)dst;
     24
     25  switch ( fgswap ) {
     26  case BR_DoNothing :  // rien a faire
     27    break;
     28  case BR_Copy : // copie directe
     29    memcpy(dst_, src, sz_);
     30    break;
     31  case BR_Swap32 : // On swappe toutes les donnees du paquet
     32  //  les bytes sont dans l'ordre par paquet de 4 octets (Int32) , les deux Int32 de
     33  // On copie la zone donnees en faisant un byte-swap correspondant a 8 octets (4->8 Reza/firmware SGDMA)
     34    for(int ka=0; ka<sz_/4; ka+=2) {
     35      dst32[ka] = src32[ka+1];
     36      dst32[ka+1] = src32[ka];
     37    }
     38   
     39    break;
     40  case BR_SwapAll:
     41  for(int ka=0; ka<sz_; ka+=8) {
     42      for(int kb=0; kb<4; kb++) {
     43        dst_[ka+kb] = src[ka+3-kb+4];   
     44        dst_[ka+kb+4] = src[ka+3-kb];   
     45      }
     46    }
     47    for(int ka=HeaderSize()+DataSize(); ka<sz_; ka+=8) {
     48      for(int kb=0; kb<4; kb++) {
     49        dst_[ka+kb] = src[ka+3-kb+4];   
     50        dst_[ka+kb+4] = src[ka+3-kb];   
     51      }
     52    }
     53    break;
     54  case BR_SwapHDR :
     55  case BR_FFTOneChan :
     56  case BR_FFTTwoChan :
     57    // ByteSwap 8 (4->8 Reza/firmware SGDMA) de l'enete
     58    for(int ka=0; ka<BRHDRSIZE; ka+=8) {
     59      for(int kb=0; kb<4; kb++) {
     60        dst_[ka+kb] = src[ka+3-kb+4];   
     61        dst_[ka+kb+4] = src[ka+3-kb];   
     62      }
     63    }
     64
     65    if (fgswap == BR_FFTOneChan) ReorderFFTData(src+HeaderSize(), dst_+HeaderSize(), DataSize());
     66    else if (fgswap == BR_FFTTwoChan) {
     67      ReorderFFTData(src+HeaderSize(), dst_+HeaderSize(), DataSize()/2);
     68      ReorderFFTData(src+HeaderSize()+DataSize()/2, dst_+HeaderSize()+DataSize()/2, DataSize()/2);
     69    }
     70    for(int ka=HeaderSize()+DataSize(); ka<sz_; ka+=8) {
     71      for(int kb=0; kb<4; kb++) {
     72        dst_[ka+kb] = src[ka+3-kb+4];   
     73        dst_[ka+kb+4] = src[ka+3-kb];   
     74      }
     75    }
     76    break;
     77 
     78  case BR_FFTOneChan32 :
     79  case BR_FFTTwoChan32 :
     80    // on swappe tout  en attendant le swap general sinon il faut encore créer une fonctiondifferente de ReorderFFT
     81    // swapp du header uniquement
     82    for(int ka=0; ka<BRHDRSIZE/4; ka+=2) {
     83      dst32[ka] = src32[ka+1];
     84      dst32[ka+1] = src32[ka];
     85    }
     86 
     87    // on reoordonne et on swappe en mem temps
     88    if (fgswap == BR_FFTOneChan32) ReorderFFTData32(src+HeaderSize(), dst_+HeaderSize(), DataSize());
     89    else if (fgswap == BR_FFTTwoChan32) {
     90      ReorderFFTData32(src+HeaderSize(), dst_+HeaderSize(), DataSize()/2);
     91      ReorderFFTData32(src+HeaderSize()+DataSize()/2, dst_+HeaderSize()+DataSize()/2, DataSize()/2);
     92    }
     93    // swap du trailler uniquement
     94    for(int ka=(HeaderSize()+DataSize())/4;ka < sz_/4; ka+=2) {
     95      dst32[ka] = src32[ka+1];
     96      dst32[ka+1] = src32[ka];
     97    }
     98   
     99  case BR_FFTOneChanNoSwap :
     100  case BR_FFTTwoChanNoSwap :
     101    // on a plus de swapdonc il faut copier dans dst
     102    memcpy(dst_, src, sz_);
     103    // on reoordonne et on swappe en mem temps
     104    if (fgswap == BR_FFTOneChanNoSwap) ReorderFFTDataNoSwap(src+HeaderSize(), dst_+HeaderSize(), DataSize());
     105    else if (fgswap == BR_FFTTwoChanNoSwap) {
     106      ReorderFFTDataNoSwap(src+HeaderSize(), dst_+HeaderSize(), DataSize()/2);
     107      ReorderFFTDataNoSwap(src+HeaderSize()+DataSize()/2, dst_+HeaderSize()+DataSize()/2, DataSize()/2);
     108    }
     109   
     110    break;
     111  }  // Fin switch
     112
     113}
     114
     115/* --Methode__ */
     116BRPaquet::BRPaquet(Byte* srcdst, int paqsz)
     117{
     118  dst_ = srcdst;
     119  sz_ = paqsz;
     120  // Il faut mettre une protection (throw) si srcdst==NULL ou sz==0
     121}
     122
     123
     124/* --Methode__ */
     125UInt16 BRPaquet::ChannelID()
     126{
     127
     128  UInt16 ChnId=ChanId();
     129  UInt16 ChpId=ChipId();
     130 
     131     
     132  if (ChpId == 2)
     133    {
     134      if (ChnId == 1) ChnId = Ch3;
     135      if (ChnId == 2) ChnId = Ch4;
     136      if (ChnId == 3) ChnId = Ch3_4;
     137    }
     138  return(ChnId);
     139}
     140
     141/* --Methode__ */
     142UInt16 BRPaquet::ModeAcquisition()
     143{
     144  UInt16 ModAq;
     145  printf("Mod Acq %x \n",ModeAcq());
     146  ModAq = ((ModeAcq() & 0x30)>> 4);
     147  return(ModAq);
     148 
     149}
     150
     151/* --Methode__ */
     152void BRPaquet::Print(ostream & os, int nelt, bool prht)
     153{
     154  os << endl << "BRPaquet::Print() PaqSz=" << PaquetSize() << " DataSize=" << DataSize()
     155     << " dst_pointer=(hex)" << hex << (unsigned long)dst_ << dec << endl;
     156  if (dst_ == NULL) {
     157    os << " ...NULL paquet " << endl;
     158    return;
    24159  }
    25   // On copie le trailer sans swapper
    26   for(int k=sz_-1; k>=sz_-BRTRLSIZE; k--)  dst_[k] = src[k];
    27 
    28 }
    29 
    30 
    31 void BRPaquet::Print(ostream & os, int nelt, bool prht)
    32 {
    33   os << "BRPaquet::Print() PaqSz=" << PaquetSize() << " DataSize=" << DataSize()
    34      << " dst_pointer=(hex)" << hex << (unsigned long)dst_ << dec << endl;
     160
     161  os << endl << " BR AcqMode: " << ModeAcquisition() << " Channel: " << ChannelID()  << endl;
     162
    35163  if (TrailerSize() > 0)
    36164    os << " ...HDRMarker(hex)=" << hex <<  HDRMarker() << " TRLMarker=" << TRLMarker() << dec << endl;
    … …  
    38166    os << " ...HDRMarker(hex)=" << hex <<  HDRMarker() << " NO TRLMarker=" << dec << endl;
    39167  UInt32 tt1, tt2;
    40 //  tt2 = (TimeTag()&0xFFFFFFFF00000000) >> 32;
    41   tt2=0;
    42   tt1 = TimeTag()&0x00000000FFFFFFFF;
    43   os << " ...TimeTag (hex)=" << hex << TimeTag() << " TT1= " << tt1 << " TT2=" << tt2 << dec << endl;
     168  tt2 = TimeTag1();
     169  tt1 = TimeTag2();
     170  os << " ...TimeTag (hex)=" << hex << "TimeTag()" << " TT1= " << tt1 << " TT2=" << tt2 << dec << endl;
     171  // os << " ...Position Chariot (hex)= " << hex << PositionChariot() << endl;
    44172  if (nelt > DataSize()/2) nelt = DataSize()/2;
    45173  os << " ...Data[1.." << nelt << "]= ";
    46   for(int k=0; k<nelt; k++) os << (int)(*(Data()+k)) << " , ";
     174
     175  for(int k=0; k<nelt; k++) os << (int)(*(Data1()+k)) << " , ";
    47176  os << endl;
    48177  os << " ...Data[" << DataSize()-nelt << ".." << DataSize()-1 << "]= ";
    49   for(int k=DataSize()-nelt; k<DataSize(); k++) os << (int)(*(Data()+k)) << " , ";
     178  for(int k=DataSize()-nelt; k<DataSize(); k++) os << (int)(*(Data1()+k)) << " , ";
    50179  os << endl;
    51180  if (prht) {   // Impression header / trailer
    … …  
    64193  }
    65194}
     195
     196
     197// ---------------------------------------------------------
     198//  **** REMARQUE   N/2+1 complexes -> N/2 complexes *****
     199// N = Nb d'echantillon en temps -> N/2 paires (real, imag)
     200// Il y a le continu, et N/2 frequences ---> N/2+1 nombres complexes,
     201// mais avec la contrainte Z(0).imag = 0  Z(N/2).imag = 0
     202// f(i) i=0...N-1   ===> Z(k) ( Z complexe , k=0...N/2 )
     203// mais avec la contrainte Z(0).imag = 0  Z(N/2).imag = 0
     204// On peut donc tout mettre ds N/2 complexes en choisissant
     205// de mettre ds Z(0).imag  Z(N/2).real
     206// ----------------------------------------------------------
     207
     208// Fonction magique qui donne le pointeur permettant de tenir compte du byte-swp sur 8 octets
     209static inline int IndexByteSwap8(int idx)
     210{
     211  return ( (idx-(idx%8))+(7-idx%8) ) ;
     212}
     213
     214/* --Methode__ */
     215void BRPaquet::ReorderFFTData(Byte* src, Byte* dst, int N)
     216{
     217  // Code recopie depuis /Dev/DisplayData/HistoWindow.cc
     218  // fonction TraceWind::DisplayBaoDatasFFT() et adapte aux structures BRPaquet et Cie
     219  // Modif par rapport au code de Bruno : N/2 elements complexes au lieu de N/2+1 - Remarque ci-dessus
     220 
     221  int nCoef = N / 2; // to change
     222  int debutIndex = N / 4 + 1;
     223  int fifoSize =  N / 4 - 1;
     224  int i;
     225
     226  TwoByteComplex* dstcmplx = (TwoByteComplex*)dst;
     227 
     228  //   cout << " Display BAO Datas FFT (" << N << ")" << " : from 0 to "<< nCoef << endl;
     229  //   cout << " Variables : debutIndex, fifoSize " << debutIndex << ", " << fifoSize << endl;
     230
     231 
     232  // Sortie 1
     233  for (i = 0; i < fifoSize ; i++)
     234    {
     235      dstcmplx[debutIndex + i].realB() = src[IndexByteSwap8(2*i)];
     236      dstcmplx[debutIndex + i].imagB() = src[IndexByteSwap8(2*i + 1)];
     237    }
     238
     239  // element au milieu
     240  dstcmplx[N / 4].realB() = src[IndexByteSwap8(2*fifoSize)];
     241  dstcmplx[N / 4].imagB() = src[IndexByteSwap8(2*fifoSize + 1)];
     242
     243  // Sortie 2
     244  for (i = 0; i < fifoSize ; i++)
     245    {
     246      dstcmplx[fifoSize - i].realB() = src[IndexByteSwap8(nCoef + 2*i)];
     247      dstcmplx[fifoSize - i].imagB() = src[IndexByteSwap8(nCoef + 2*i + 1)];
     248    }
     249
     250  // k = 0 et k = N/2 
     251  dstcmplx[0].realB() = src[IndexByteSwap8(N - 2)];
     252  // Voir Remarque ci-dessus Z(N/2).real -> Z(0).image
     253  dstcmplx[0].imagB() = src[IndexByteSwap8(N - 1)];  // Attention, on met ici la real(fmax)
     254
     255  return ;
     256}
     257
     258static inline int IndexByteSwap8_32(int idx)
     259{
     260  return ( (idx-(idx%8))+((4+idx%8)%8) ) ;
     261}
     262
     263void BRPaquet::ReorderFFTData32(Byte* src, Byte* dst, int N)
     264{
     265  // Code recopie depuis /Dev/DisplayData/HistoWindow.cc
     266  // fonction TraceWind::DisplayBaoDatasFFT() et adapte aux structures BRPaquet et Cie
     267  // Modif par rapport au code de Bruno : N/2 elements complexes au lieu de N/2+1 - Remarque ci-dessus
     268
     269  int nCoef = N / 2; // to change
     270  int debutIndex = N / 4 + 1;
     271  int fifoSize =  N / 4 - 1;
     272  int i;
     273
     274  TwoByteComplex* dstcmplx = (TwoByteComplex*)dst;
     275 
     276  //   cout << " Display BAO Datas FFT (" << N << ")" << " : from 0 to "<< nCoef << endl;
     277  //   cout << " Variables : debutIndex, fifoSize " << debutIndex << ", " << fifoSize << endl;
     278
     279 
     280  // Sortie 1
     281  for (i = 0; i < fifoSize ; i++)
     282    {
     283      dstcmplx[debutIndex + i].realB() = src[IndexByteSwap8_32(2*i)];
     284      dstcmplx[debutIndex + i].imagB() = src[IndexByteSwap8_32(2*i + 1)];
     285    }
     286
     287  // element au milieu
     288  dstcmplx[N / 4].realB() = src[IndexByteSwap8_32(2*fifoSize)];
     289  dstcmplx[N / 4].imagB() = src[IndexByteSwap8_32(2*fifoSize + 1)];
     290
     291  // Sortie 2
     292  for (i = 0; i < fifoSize ; i++)
     293    {
     294      dstcmplx[fifoSize - i].realB() = src[IndexByteSwap8_32(nCoef + 2*i)];
     295      dstcmplx[fifoSize - i].imagB() = src[IndexByteSwap8_32(nCoef + 2*i + 1)];
     296    }
     297
     298  // k = 0 et k = N/2 
     299  dstcmplx[0].realB() = src[IndexByteSwap8_32(N - 2)];
     300  // Voir Remarque ci-dessus Z(N/2).real -> Z(0).image
     301  dstcmplx[0].imagB() = src[IndexByteSwap8_32(N - 1)];  // Attention, on met ici la real(fmax)
     302
     303  return ;
     304}
     305void BRPaquet::ReorderFFTDataNoSwap(Byte* src, Byte* dst, int N)
     306{
     307  // Code recopie depuis /Dev/DisplayData/HistoWindow.cc
     308  // fonction TraceWind::DisplayBaoDatasFFT() et adapte aux structures BRPaquet et Cie
     309  // Modif par rapport au code de Bruno : N/2 elements complexes au lieu de N/2+1 - Remarque ci-dessus
     310 
     311  int nCoef = N / 2; // to change
     312  int debutIndex = N / 4 + 1;
     313  int fifoSize =  N / 4 - 1;
     314  int i;
     315
     316  TwoByteComplex* dstcmplx = (TwoByteComplex*)dst;
     317 
     318  //   cout << " Display BAO Datas FFT (" << N << ")" << " : from 0 to "<< nCoef << endl;
     319  //   cout << " Variables : debutIndex, fifoSize " << debutIndex << ", " << fifoSize << endl;
     320
     321 
     322  // Sortie 1
     323  for (i = 0; i < fifoSize ; i++)
     324    {
     325      dstcmplx[debutIndex + i].realB() = src[(2*i)];
     326      dstcmplx[debutIndex + i].imagB() = src[(2*i + 1)];
     327    }
     328
     329  // element au milieu
     330  dstcmplx[N / 4].realB() = src[(2*fifoSize)];
     331  dstcmplx[N / 4].imagB() = src[(2*fifoSize + 1)];
     332
     333  // Sortie 2
     334  for (i = 0; i < fifoSize ; i++)
     335    {
     336      dstcmplx[fifoSize - i].realB() = src[(nCoef + 2*i)];
     337      dstcmplx[fifoSize - i].imagB() = src[(nCoef + 2*i + 1)];
     338    }
     339
     340  // k = 0 et k = N/2 
     341  dstcmplx[0].realB() = src[(N - 2)];
     342  // Voir Remarque ci-dessus Z(N/2).real -> Z(0).image
     343  dstcmplx[0].imagB() = src[(N - 1)];  // Attention, on met ici la real(fmax)
     344
     345  return ;
     346}
     347
     348/* --Methode__ */
     349const char* BRPaquet::FmtConvToString(BRDataFmtConv fgswap)
     350{
     351  const char * rs="";
     352  switch ( fgswap ) {
     353  case BR_DoNothing : 
     354    rs = "BR_DoNothing";
     355    break;
     356  case BR_Copy :
     357    rs = "BR_Copy";
     358    break;
     359  case BR_SwapAll :
     360    rs = "BR_SwapAll";
     361    break;
     362  case BR_SwapHDR :
     363    rs = "BR_SwapHDR";
     364    break;
     365  case BR_FFTOneChan :
     366    rs = "BR_FFTOneChan";
     367    break;
     368  case BR_FFTTwoChan :
     369    rs = "BR_FFTTwoChan";
     370    break;
     371  case BR_Swap32 :
     372    rs = "BR_Swap32";
     373    break;
     374  case BR_FFTOneChan32 :
     375    rs = "BR_FFTOneChan32";
     376    break;
     377  case BR_FFTTwoChan32 :
     378    rs = "BR_FFTTwoChan32";
     379    break;
     380  case BR_FFTOneChanNoSwap :
     381    rs = "BR_FFTOneChanNoSwap";
     382    break;
     383  case BR_FFTTwoChanNoSwap :
     384    rs = "BR_FFTTwoChanNoSwap";
     385    break;
     386  default:
     387    rs = "?????";
     388    break;
     389  }  // Fin switch
     390  return rs;
     391}
     392
     393// --------------------------------------------------------------------------
     394// Classe pour effectuer des verifications d'integrite sur les paquets/frames
     395// --------------------------------------------------------------------------
     396
     397BRPaqChecker::BRPaqChecker()
     398{
     399  totnframes = 0;
     400  nframeok = 0;
     401  lostframes = 0;
     402  frclst = 0;
     403}
     404
     405BRPaqChecker::~BRPaqChecker()
     406{
     407}
     408
     409bool BRPaqChecker::Check(BRPaquet& paq)
     410{
     411  totnframes++;
     412  if (paq.HDRMarker() != 0x76543210) return false;
     413  unsigned int curfc = paq.FrameCounter();
     414  unsigned int delfc = 0;
     415  if (nframeok > 0) {
     416    if (curfc>frclst)  delfc = (curfc-frclst);
     417    else delfc = (65535-frclst+curfc);
     418    lostframes += (unsigned long long)delfc - 1;
     419  }
     420  nframeok++; frclst = curfc;
     421  return true;
     422}
     423
     424ostream& BRPaqChecker::Print(ostream& os)
     425{
     426  os << "BRPaqChecker:  Tot.Nb.Frames.Proc=" << totnframes << " NbFrameOK=" << nframeok
     427     << " LostFrames=" << lostframes
     428     << " Loss=" << (double)lostframes*100./(double)totnframes << " %" << endl;
     429  return os;
     430}

  • trunk/AddOn/TAcq/brpaqu.h

    r3538 r3592  
    1212#include "brtypes.h"
    1313
     14
    1415using namespace std;
    1516
     17// On definit une classe TwoByteComplex() pour manipuler une paire de byte
     18// representant partie relle et imaginaire d'un nombre complexe
     19// Remarque Byte = unsigned char
     20// Remarque SByte = char = signed char
     21class TwoByteComplex {
     22 public:
     23  explicit TwoByteComplex() { val_[0] = val_[1] = 0; }
     24  //  explicit TwoByteComplex(TwoByteComplex const & a) { val_[0] = a.val_[0]; val_[1] = 0; }
     25  explicit TwoByteComplex(Byte re, Byte im) { val_[0] = re; val_[1] = im; }
     26  explicit TwoByteComplex(SByte re, SByte im) { val_[0] = re; val_[1] = im; }
     27  explicit TwoByteComplex(int re, int im) { val_[0] = re; val_[1] = im; }
     28  explicit TwoByteComplex(float re, float im) { val_[0] = re; val_[1] = im; }
     29  explicit TwoByteComplex(double re, double im) { val_[0] = re; val_[1] = im; }
     30 
     31  inline Byte& realB() { return  val_[0]; }
     32  inline Byte& imagB() { return  val_[1]; }
     33
     34  inline SByte realSB() { return  ((SByte*)val_)[0]; }
     35  inline SByte imagSB() { return  ((SByte*)val_)[1]; }
     36
     37  inline int realI() { return  ((SByte*)val_)[0]; }
     38  inline int imagI() { return  ((SByte*)val_)[1]; }
     39
     40  inline double realD() { return  ((SByte*)val_)[0]; }
     41  inline double imagD() { return  ((SByte*)val_)[1]; }
     42
     43  Byte val_[2];
     44};
     45
     46// Actuellement apres test Nancay un mot de 4octects est mis avant le header
     47#define OFFSET 0
    1648// Taille entete/trailer en octets
    1749#define BRHDRSIZE 24
    1850// Actuellement le trailer est vire ...
    19 /* #define BRTRLSIZE 16  */
    20 #define BRTRLSIZE 0
     51// Reza, 28 Jan 2009  le trailer fait 16 octets
     52#define BRTRLSIZE 16 
     53/* REZA passe a 16 octets  #define BRTRLSIZE 40 */
     54//Offset Mode Acq
     55#define BRMODACQOFF  20 /*RezaMOD #define BRMODACQOFF  24 */
     56// Offset ChanId
     57#define BRCHANIDOFF 20 /*RezaMOD #define BRCHANIDOFF 24 */
    2158// Offset du time-tag
    22 #define BRTMTAGOFF 8
     59#define BRTMTAGOFF 12 /*RezaMOD #define BRTMTAGOFF 16 */
     60// offset FrameCounter
     61#define BRFRCPTOFF 16 /*RezaMOD#define BRFRCPTOFF 20 */
     62// offset Paquet Id
     63#define BRPKTIDOFF 20  /*RezaMOD #define BRPKTIDOFF 24 */
     64// offset paquet lenght = FrameDataLength
     65#define BRPKTLENOFF 8 /*RezaMOD #define BRPKTLENOFF 12 */
     66// offset position chariot ???
     67//#define BRPCOFF 16
     68
     69enum ChannelID
     70  {
     71    None=0,
     72    Ch1,
     73    Ch2,
     74    Ch1_2,
     75    Ch3,
     76    Ch4,
     77    Ch3_4
     78  };
     79enum ModeAcq
     80  {
     81    RawData=1,
     82    FFT
     83  };
     84
     85// Definition des actions sur la conversion de format (swap...) des donnees arrivant
     86enum BRDataFmtConv {
     87  BR_DoNothing,
     88  BR_Copy,
     89  BR_SwapAll,
     90  BR_SwapHDR,
     91  BR_FFTOneChan,
     92  BR_FFTTwoChan,
     93  BR_Swap32,
     94  BR_FFTOneChan32,
     95  BR_FFTTwoChan32,
     96  BR_FFTOneChanNoSwap,
     97  BR_FFTTwoChanNoSwap
     98};
    2399
    24100// Structure correspondant a HEADER-DATA-TRAILER
    25101class BRPaquet {
    26102 public:
    27   BRPaquet(Byte* src, Byte* dst, int paqsz=4096);
     103  BRPaquet(Byte* src, Byte* dst, int paqsz=4096, BRDataFmtConv fgswap=BR_SwapAll);
     104  BRPaquet(Byte* srcdst, int paqsz=4096);
    28105  //  ~BRPaquet();
    29106
    … …  
    31108  inline int PaquetSize() { return sz_; }
    32109  inline int DataSize() { return sz_-(BRHDRSIZE+BRTRLSIZE); }
    33 
    34110  inline int HeaderSize() { return BRHDRSIZE; }
    35111  inline int TrailerSize() { return BRTRLSIZE; }
    36112
    37113  // Acces differentes zone memoire
    38   inline Byte* Data() { return dst_+BRHDRSIZE; }
    39   inline Byte* Header() { return dst_; }
    40   inline Byte* Trailer() { return (dst_+sz_-BRTRLSIZE); }
     114  inline Byte* Data() { return dst_+BRHDRSIZE+OFFSET; }
     115  inline Byte* Data1() { return dst_+BRHDRSIZE+OFFSET; }
     116  inline Byte* Data2() { return dst_+BRHDRSIZE+(DataSize()/2)+OFFSET; }
     117  inline Byte* Header() { return dst_+OFFSET; }
     118  inline Byte* Trailer() { return (dst_+sz_-BRTRLSIZE+OFFSET); }
    41119 
    42120  // Valeurs differentes zones HDR/TRL
    43   inline UInt64 HDRMarker() {return *((UInt64*)dst_);}
    44   inline UInt64 TRLMarker() {return *((UInt64*)(dst_+sz_-BRTRLSIZE));}
    45   inline UInt64 TimeTag() {return *((UInt64*)(dst_+BRTMTAGOFF));}
     121  inline UInt32 HDRMarker() {return *((UInt32*)(dst_+OFFSET));}
     122  inline UInt32 TRLMarker() {return *((UInt32*)(dst_+(sz_-BRTRLSIZE+OFFSET)));}
     123  inline UInt16 ModeAcq()   {return *((UInt16*)(dst_+(BRMODACQOFF+OFFSET)));}
     124  inline UInt16 ChanId()    {return (( *((UInt16*)(dst_+(BRCHANIDOFF+OFFSET))) & 0x1800)>> 12) ;}
     125  inline UInt16 ChipId()    {return (( *((UInt16*)(dst_+(BRCHANIDOFF+OFFSET))) & 0xC00)>> 10) ;}
     126
     127  inline UInt32 TimeTag1()   {return *((UInt32*)(dst_+(BRTMTAGOFF+OFFSET)));}
     128  inline UInt32 FrameCounter()   {return ((*((UInt32*)(dst_+(BRFRCPTOFF+OFFSET))) &0xFFFF0000) >> 16);}
     129  inline UInt32 TimeTag2() {return (*((UInt32*)(dst_+(BRFRCPTOFF+OFFSET))) &0xFFFF);}
     130
     131  inline UInt32 PaqId()     {return *((UInt32*)(dst_+(BRPKTIDOFF+OFFSET)));}
     132  inline UInt16 PaqLen()    {return *((UInt16*)(dst_+(BRPKTLENOFF+OFFSET)));}
     133
     134
     135  // inline unsigned short PositionChariot() {return *((unsigned short*)(dst_+BRPCOFF+OFFSET));}
    46136
    47137  // pour faire un print de la structure
    48138  void Print(ostream & os, int nelt=8, bool prht=true);
     139  UInt16  ChannelID();
     140  UInt16  ModeAcquisition();
    49141
    50   // protected:
     142  // fonction appelee par le constructeur pour reordonner les donnees FFT
     143  static void ReorderFFTData(Byte* src, Byte* dst, int sz);
     144  static void ReorderFFTData32(Byte* src, Byte* dst, int sz);
     145  static void ReorderFFTDataNoSwap(Byte* src, Byte* dst, int sz);
     146  static const char* FmtConvToString(BRDataFmtConv fgswap);
     147// protected:
    51148  // donnees membres
    52149  int sz_; //taille du paquet
    … …  
    54151};
    55152
     153// --------------------------------------------------------------------------
     154// Classe pour effectuer des verifications d'integrite sur les paquets/frames
     155// --------------------------------------------------------------------------
     156
     157class BRPaqChecker {
     158public:
     159  BRPaqChecker();
     160  ~BRPaqChecker();
     161  // Verifie le paquet, renvoie true si OK
     162  bool Check(BRPaquet& paq);
     163  // Imprime le compte de paquets ...
     164  ostream & Print(ostream& os);
     165
     166  unsigned long long totnframes;    // Nombre totale de frames/paquets traites
     167  unsigned long long  nframeok;     // Nombre totale de frames/paquets avec HDR/TRL OK
     168  unsigned long long lostframes;    // Nombre totale de frames/paquets perdus
     169  unsigned int frclst;         // derniere valeur du frame-counter
     170};
     171
    56172#endif
     173 
     174 

  • trunk/AddOn/TAcq/makefile

    r3591 r3592  
    11include $(SOPHYABASE)/include/sophyamake.inc
    22
    3 all : traidio  tmtfft tstminifits tpciew tbrpaq tmtacq tstrdfits mfits2spec tsok sigana
     3all : traidio  tmtfft tstminifits tpciew tbrpaq tmtacq tstrdfits mfits2spec tsok
    44
    55clean :
    66        rm *.o traidio tmtfft tstminifits tpciew tbrpaq tmtacq tstrdfits
    7         rm mfits2spec tsok sigana *.ppf *.fits
     7        rm mfits2spec tsok *.ppf *.fits
    88
    99########################################################
    … …  
    4646
    4747##   --------------
    48 sigana : sigana.o sigana.o
    49         $(CXXLINK) -o sigana sigana.o minifits.o  $(SOPHYAEXTSLBLIST)
    5048
    51 sigana.o : sigana.cc minifits.h
    52         $(CXXCOMPILE) -o sigana.o sigana.cc
     49mfits2spec : mfits2spec.o minifits.o brpaqu.o
     50        $(CXXLINK) -o mfits2spec mfits2spec.o minifits.o brpaqu.o $(SOPHYAEXTSLBLIST)
    5351
    54 mfits2spec : mfits2spec.o minifits.o
    55         $(CXXLINK) -o mfits2spec mfits2spec.o minifits.o  $(SOPHYAEXTSLBLIST)
    56 
    57 mfits2spec.o : mfits2spec.cc minifits.h
     52mfits2spec.o : mfits2spec.cc minifits.h brpaqu.h
    5853        $(CXXCOMPILE) -o mfits2spec.o mfits2spec.cc
    5954

  • trunk/AddOn/TAcq/mfits2spec.cc

    r3538 r3592  
    22#include "sopnamsp.h"
    33#include "machdefs.h"
     4
     5/* ------------------------------------------------------------------
     6   Programme de calcul de spectre moyenne a partir des fichiers fits 
     7   d'acquisition de BAORadio - donnees brutes (non FFT)
     8   R. Ansari, C. Magneville
     9   V1 : Juillet 2008  , V2 : Avril 2009
     10   ------------------------------------------------------------------ */
    411
    512// include standard c/c++
    … …  
    2633
    2734
    28 // include mini-fits lib
     35// include mini-fits lib , et structure paquet BAORadio
    2936#include "minifits.h"
    30 
    31 
    32 inline r_4 Zmod2(complex<r_4> z)
    33 { return (z.real()*z.real()+z.imag()*z.imag()); }
    34 
     37#include "brpaqu.h"
     38
     39//---- Declaration des fonctions de calcul ----
     40// Fonction d'analyse 1ere version, pas d'entete ds le fichier, 1 voie
     41int ana_data_0(vector<string>& infiles, string& outfile);
     42// Fonction d'analyse 2eme version , 1 voie / paquet
     43int ana_data_1(vector<string>& infiles, string& oufile);
     44// Fonction d'analyse 2eme version , 2 voies / paquet
     45int ana_data_2(vector<string>& infiles, string& oufile);
     46
     47//----------------------------------------------------
     48//----------------------------------------------------
    3549int main(int narg, char* arg[])
    3650{
    37   if (narg < 2) {
    38     cout << " Usage: mfits2spec file1.fits [file2.fits ...] " << endl;
     51  if (narg < 4) {
     52    cout << " ---Calcul spectres moyennes a partir de fits BAORadio " << endl;
     53    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; 
    3956    return 1;
    4057  }
    41   cout << " ---------- mfits2spec.cc Start ------------- " << endl;
    4258
    4359  TArrayInitiator  _inia;
    … …  
    4561  int rc = 0;
    4662  try {
    47         ResourceUsage resu;
    48         TVector<float> spectre;
    49         sa_size_t nzm = 0;   // Nb de spectres moyennes
    50     for(int ifile=1; ifile<narg; ifile++) {
    51       string ffname = arg[ifile];
    52 // -------------- Lecture de bytes     
    53       cout << ifile <<"-Ouverture/lecture fichier " << ffname << endl;
    54       MiniFITSFile mff(ffname, MF_Read);
    55       cout << "... Type=" << mff.DataTypeToString() << " NAxis1=" << mff.NAxis1()
    56            << " NAxis2=" << mff.NAxis2() << endl;
    57       if (mff.DataType() != MF_Byte) {
    58         cout << " PB : DataType!=MF_Byte --> skipping " << endl;
    59       }
    60       size_t sx = mff.NAxis1();
    61       size_t sy = mff.NAxis2();
    62 
    63       Byte* data = new Byte[sx];
    64       TVector<r_4> vx(sx);
    65       TVector< complex<r_4> > cfour;
    66       FFTPackServer ffts;
    67      
    68       for(int j=0; j<sy; j++) {
    69             mff.ReadB(data, sx, j*sx);
    70         // On convertit en float
    71         for(sa_size_t ix=0; ix<vx.Size(); ix++)  vx(ix) = (r_4)(data[ix]);
    72         // On fait le FFT
    73         ffts.FFTForward(vx, cfour);
    74         if (!spectre.IsAllocated())  spectre.SetSize(cfour.Size());
    75        
    76         // On cumule les modules carres
    77         for(sa_size_t jf=1; jf<spectre.Size(); jf++)
    78           spectre(jf) += Zmod2(cfour(jf));
    79         nzm++;
    80 
    81 //      cout << " j=" << j << " data[0...3]=" << (int)data[0] << " , "
    82 //           << (int)data[1] << " , " <<  (int)data[2] << endl;
    83       }
    84       cout << "---- FIN lecture " << ffname << endl;     
    85     }   
    86     cout << "---- Ecriture spectre moyenne ds mspectre.ppf " << endl;     
    87     spectre /= (r_4)(nzm);
    88     POutPersist po("mspectre.ppf");
    89     po << spectre;
     63    string act = arg[1];
     64    string outppf = arg[2];
     65    vector<string> infiles;
     66    for(int i=3; i<narg; i++) infiles.push_back(arg[i]);
     67    cout << " ---------- mfits2spec.cc Start - ACT= " << act << " ------------- " << endl;
     68    ResourceUsage resu;
     69    if (act == "-0") ana_data_0(infiles, outppf);
     70    else if (act == "-1") ana_data_1(infiles, outppf);
     71    else if (act == "-2") ana_data_2(infiles, outppf);
     72    else cout << " mfits2spec.cc / Bad argument ACT=" << act << " -> exit" << endl;
    9073    cout << resu ;
    9174  }
    … …  
    10992
    11093}
     94
     95inline r_4 Zmod2(complex<r_4> z)
     96{ return (z.real()*z.real()+z.imag()*z.imag()); }
     97
     98/*--Nouvelle-Fonction--*/
     99int ana_data_0(vector<string>& infiles, string& outfile)
     100{
     101  TVector<float> spectre;
     102  sa_size_t nzm = 0;   // Nb de spectres moyennes
     103  for(int ifile=0; ifile<infiles.size(); ifile++) {
     104    string ffname = infiles[ifile];
     105// -------------- Lecture de bytes     
     106    cout << ifile <<"-Ouverture/lecture fichier " << ffname << endl;
     107    MiniFITSFile mff(ffname, MF_Read);
     108    cout << "... Type=" << mff.DataTypeToString() << " NAxis1=" << mff.NAxis1()
     109         << " NAxis2=" << mff.NAxis2() << endl;
     110    if (mff.DataType() != MF_Byte) {
     111      cout << " PB : DataType!=MF_Byte --> skipping " << endl;
     112    }
     113    size_t sx = mff.NAxis1();
     114    size_t sy = mff.NAxis2();
     115
     116    Byte* data = new Byte[sx];
     117    TVector<r_4> vx(sx);
     118    TVector< complex<r_4> > cfour;
     119    FFTPackServer ffts;
     120     
     121    for(int j=0; j<sy; j++) {
     122      mff.ReadB(data, sx, j*sx);
     123      // On convertit en float
     124      for(sa_size_t ix=0; ix<vx.Size(); ix++)  vx(ix) = (r_4)(data[ix]);
     125      // On fait le FFT
     126      ffts.FFTForward(vx, cfour);
     127      if (!spectre.IsAllocated())  spectre.SetSize(cfour.Size());
     128       
     129      // On cumule les modules carres
     130      for(sa_size_t jf=1; jf<spectre.Size(); jf++)
     131        spectre(jf) += Zmod2(cfour(jf));
     132      nzm++;
     133
     134//      cout << " j=" << j << " data[0...3]=" << (int)data[0] << " , "
     135//           << (int)data[1] << " , " <<  (int)data[2] << endl;
     136      }
     137    cout << "---- FIN lecture " << ffname << endl;     
     138    delete[] data;
     139  }   
     140  cout << "---- Ecriture spectre moyenne ds " << outfile << endl;     
     141  spectre /= (r_4)(nzm);
     142  POutPersist po(outfile);
     143  po << spectre;
     144  return 0;
     145}
     146
     147/*--Nouvelle-Fonction--*/
     148int ana_data_1(vector<string>& infiles, string& outfile)
     149{
     150  TVector<float> spectre;
     151  float freq0 = 0; 
     152  int paqsz = 0;
     153  int nfileok;
     154  sa_size_t nzm = 0;   // Nb de spectres moyennes
     155  Byte* data = NULL;
     156  FFTPackServer ffts;
     157  for(int ifile=0; ifile<infiles.size(); ifile++) {
     158    string ffname = infiles[ifile];
     159// -------------- Lecture de bytes     
     160    cout << "ana_data_1[" << ifile <<"]Ouverture/lecture fichier " << ffname << endl;
     161    MiniFITSFile mff(ffname, MF_Read);
     162    cout << "... Type=" << mff.DataTypeToString() << " NAxis1=" << mff.NAxis1()
     163         << " NAxis2=" << mff.NAxis2() << endl;
     164    if (mff.DataType() != MF_Byte) {
     165      cout << " PB : DataType!=MF_Byte --> skipping " << endl;
     166      continue;
     167    }
     168// Les fichier FITS contiennent l'entet (24 bytes), mais pas le trailer (16 bytes) ...
     169    if (paqsz == 0)  {  // premier passage, on fixe la taille de paquet et on alloue le buffer
     170      paqsz = mff.NAxis1()+16;
     171      data = new Byte[paqsz];
     172      for(int ib=0; ib<paqsz; ib++) data[ib]=0;
     173    }
     174    else {
     175      if (paqsz != mff.NAxis1()+16) {
     176      cout << " PB : paqsz=" << paqsz << " != mff.NAxis1()+16 --> skipping " << endl;
     177      continue;
     178      }
     179    }
     180    size_t sx = mff.NAxis1();
     181    size_t sy = mff.NAxis2();
     182    BRPaquet paq(NULL, data, paqsz);
     183    TVector<r_4> vx(paq.DataSize());
     184    TVector< complex<r_4> > cfour;
     185     
     186    for(int j=0; j<sy; j++) {
     187      mff.ReadB(data, sx, j*sx);
     188      // On convertit en float
     189      for(sa_size_t ix=0; ix<vx.Size(); ix++)  vx(ix) = (r_4)(paq.Data1()[ix])-127.5;
     190      // On fait le FFT
     191      ffts.FFTForward(vx, cfour);
     192      if (!spectre.IsAllocated())  spectre.SetSize(cfour.Size());
     193       
     194      // On cumule les modules carres  - en sautant la frequence zero
     195      for(sa_size_t jf=1; jf<spectre.Size(); jf++)
     196        spectre(jf) += Zmod2(cfour(jf));
     197      nzm++;  freq0 += cfour(0).real();
     198      }
     199    nfileok++;
     200    cout << "---- FIN lecture " << ffname << " NFileOK=" << nfileok << endl;     
     201  }   
     202  if (data) delete[] data;
     203  if (nzm <= 0) {
     204    cout << " ana_data_1/ERROR : nzm=" << nzm << " spectres moyennes !" << endl;
     205    return nzm;
     206  }
     207  cout << "---- Ecriture spectre moyenne ds " << outfile  << endl;     
     208  spectre /= (r_4)(nzm);
     209  spectre.Info().Comment() = " SpectreMoyenne (Moyenne module^2) ";
     210  spectre.Info()["NMOY"] = nzm; // Nombre de spectres moyennes
     211  spectre.Info()["Freq0"] = freq0;
     212  POutPersist po(outfile);
     213  po << PPFNameTag("specV1") << spectre;
     214  return nfileok;
     215}
     216
     217/*--Nouvelle-Fonction--*/
     218int ana_data_2(vector<string>& infiles, string& outfile)
     219{
     220  TVector<float> specV1, specV2;
     221  TVector< complex<float> > cxspecV12;
     222  float freq0v1 = 0; 
     223  float freq0v2 = 0; 
     224  int paqsz = 0;
     225  int nfileok;
     226  sa_size_t nzm = 0;   // Nb de spectres moyennes
     227  Byte* data = NULL;
     228  FFTPackServer ffts;
     229  for(int ifile=0; ifile<infiles.size(); ifile++) {
     230    string ffname = infiles[ifile];
     231// -------------- Lecture de bytes     
     232    cout << "ana_data_2[" << ifile <<"]Ouverture/lecture fichier " << ffname << endl;
     233    MiniFITSFile mff(ffname, MF_Read);
     234    cout << "... Type=" << mff.DataTypeToString() << " NAxis1=" << mff.NAxis1()
     235         << " NAxis2=" << mff.NAxis2() << endl;
     236    if (mff.DataType() != MF_Byte) {
     237      cout << " PB : DataType!=MF_Byte --> skipping " << endl;
     238      continue;
     239    }
     240// Les fichier FITS contiennent l'entet (24 bytes), mais pas le trailer (16 bytes) ...
     241    if (paqsz == 0)  {  // premier passage, on fixe la taille de paquet et on alloue le buffer
     242      paqsz = mff.NAxis1()+16;
     243      cout << " ana_data_2/ Allocating data , PaqSz=" << paqsz << endl;
     244      data = new Byte[paqsz];
     245      for(int ib=0; ib<paqsz; ib++) data[ib]=0;
     246    }
     247    else {
     248      if (paqsz != mff.NAxis1()+16) {
     249      cout << " PB : paqsz=" << paqsz << " != mff.NAxis1()+16 --> skipping " << endl;
     250      continue;
     251      }
     252    }
     253    size_t sx = mff.NAxis1();
     254    size_t sy = mff.NAxis2();
     255    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;
     259     
     260    for(int j=0; j<sy; j++) {
     261      mff.ReadB(data, sx, j*sx);
     262     // if (j%50 == 0) cout << " *DBG* mff.ReadB() , j=" << j << endl;
     263     // 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;
     266      // On fait le FFT
     267      ffts.FFTForward(vx1, cfour1);
     268      ffts.FFTForward(vx2, cfour2);
     269      if (!specV1.IsAllocated())  specV1.SetSize(cfour1.Size());
     270      if (!specV2.IsAllocated())  specV2.SetSize(cfour2.Size());
     271      if (!cxspecV12.IsAllocated())  cxspecV12.SetSize(cfour1.Size());
     272       
     273      // On cumule les modules carres  - en sautant la frequence zero
     274      for(sa_size_t jf=1; jf<specV1.Size(); jf++) {
     275        specV1(jf) += Zmod2(cfour1(jf));
     276        specV2(jf) += Zmod2(cfour2(jf));
     277        cxspecV12(jf) += cfour1(jf)*cfour2(jf);
     278      }
     279      nzm++;  freq0v1 += cfour1(0).real();  freq0v2 += cfour2(0).real();
     280      }
     281    nfileok++;
     282    cout << "---- FIN lecture " << ffname << " NFileOK=" << nfileok << endl;     
     283  } 
     284  if (data) delete[] data;
     285  if (nzm <= 0) {
     286    cout << " ana_data_2/ERROR : nzm=" << nzm << " spectres moyennes !" << endl;
     287    return nzm;
     288  }
     289  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)";
     295  specV1.Info()["NMOY"] = specV2.Info()["NMOY"] = nzm; // Nombre de spectres moyennes
     296  specV1.Info()["Freq0"] = freq0v1;
     297  specV2.Info()["Freq0"] = freq0v2;
     298  POutPersist po(outfile);
     299  po << PPFNameTag("specV1") << specV1;
     300  po << PPFNameTag("specV2") << specV2;
     301  po << PPFNameTag("cxspecV12") << cxspecV12;
     302  return 0;
     303}

  • trunk/AddOn/TAcq/minifits.h

    r3538 r3592  
    1818  explicit MiniFITSException(const char * m) { msg = m; }
    1919  explicit MiniFITSException(const string& m) { msg = m; }
     20  virtual ~MiniFITSException() { }
    2021  virtual string const& Msg() const  {return msg;}
    2122 private:
Note: See TracChangeset for help on using the changeset viewer.