source: Sophya/trunk/AddOn/TAcq/corrcrtadc.cc@ 4036

Last change on this file since 4036 was 4009, checked in by ansari, 14 years ago

code combinaison I-Q (codage fonction CombineIQPairs) , Reza 23/06/2011
File size: 13.0 KB
Line 
1// Utilisation de SOPHYA ...
2#include "sopnamsp.h"
3#include "machdefs.h"
4
5/* ---------------------------------------------------------------
6 Projet BAORadio-21cm intensity mapping - (C) LAL/IRFU 2008-2011
7
8 CRT (Pittsburgh-CMU) DUMP ADC read/processing program to compute
9 spectra and correlations du CRT
10 R. Ansari, C. Magneville - LAL/Irfu - Juin 2011
11 --------------------------------------------------------------- */
12
13// include standard c/c++
14#include <math.h>
15#include <stdio.h>
16#include <stdlib.h>
17#include <string.h>
18
19#include <iostream>
20#include <string>
21
22// SOPHYA include files
23#include "pexceptions.h"
24#include "tvector.h"
25#include "fioarr.h"
26// #include "tarrinit.h"
27#include "ntuple.h"
28#include "datatable.h"
29#include "histinit.h"
30#include "matharr.h"
31#include "timestamp.h"
32#include "fftwserver.h"
33#include "fftpserver.h"
34#include "utilarr.h"
35#include "histats.h"
36
37// include sophya mesure ressource CPU/memoire ...
38#include "resusage.h"
39#include "ctimer.h"
40#include "timing.h"
41
42
43//--------------------------- Functions in this file -------------------
44int Usage(bool fgshort=true);
45int DecodeProc(int narg, char* arg[]);
46int ProcessADCFiles(string& inoutpath, vector<int>& adcids, string& outname, bool fgcombiq,
47 int imin, int imax, int istep, int jf1, int jf2, int nfreq);
48int ComputeCorrel(TMatrix< complex<r_4> >& cfour, bool fgcombiq,
49 TMatrix< complex<r_8> >& correl, TMatrix< complex<r_8> >& autocorrel);
50int ADCFiles2Histo(string& inoutpath, vector<int>& adcids, string& outname, int imin, int imax, int istep);
51int CombineIQPairs(TMatrix< complex<r_4> >& cfour);
52//------------------------------------------------------------------------------------------------------------
53
54/* --Fonction-- */
55int Usage(bool fgshort)
56{
57 cout << " --- corrcrtadc.cc : Read CRT-ADC dump files and process to produce correlation matrices \n"
58 << " Usage: corrcrtadc -corr/-corciq/-hval InOutPath OutFile [AdcIds] [Imin,Imax,step]\n" << endl;
59 if (fgshort) {
60 cout << " corrcrtadc -h for detailed instructions" << endl;
61 return 1;
62 }
63 cout << " -corr/-corciq/-hval : Compute correlations or Fill histograms for ADC sample values \n"
64 << " -corciq : combine IQ pairs to separate side-bands ... \n"
65 << " InOutPath : Input/Output directory name \n"
66 << " ADCIds : ADC Id's to be processed (Example: 1,2,3,4, default=0,1) \n"
67 << " OutFile : Output PPF file name \n"
68 << " Imin,Imax,IStep: Input ADC dump files sequence number \n"
69 << " FileNames=InOutPath/adcdumpNJFII.fits Imin<=II<=Imax II+=IStep , J=ADCIds \n" << endl;
70 // << " NumFreq1,NumFreq2,NBinFreq: Freq Zone and number of frequency bins for ntuple" << endl;
71 return 1;
72}
73
74//----------------------------------------------------
75//----------------------------------------------------
76int main(int narg, char* arg[])
77{
78 if ((narg>1)&&(strcmp(arg[1],"-h")==0)) return Usage(false);
79 if (narg<4) return Usage(true);
80 HiStatsInitiator _inia;
81 // TArrayInitiator _inia;
82
83 int rc = 0;
84 try {
85 ResourceUsage resu;
86 DecodeProc(narg-1, arg+1);
87 resu.Update();
88 cout << resu;
89 }
90 catch (PException& exc) {
91 cerr << " corrcrtadc.cc catched PException " << exc.Msg() << endl;
92 rc = 77;
93 }
94 catch (std::exception& sex) {
95 cerr << "\n corrcrtadc.cc std::exception :"
96 << (string)typeid(sex).name() << "\n msg= "
97 << sex.what() << endl;
98 rc = 78;
99 }
100 catch (...) {
101 cerr << " corrcrtadc.cc catched unknown (...) exception " << endl;
102 rc = 79;
103 }
104
105 cout << ">>>> corrcrtadc.cc ------- END ----------- RC=" << rc << endl;
106 return rc;
107
108}
109
110
111//--------------------------------------------------------------------
112// Traitement fichiers produits par vismfib (V Nov09)
113/* --Fonction-- */
114int DecodeProc(int narg, char* arg[])
115{
116 // Decodage des arguments et traitement
117 string popt = arg[0];
118 bool fillhis=false;
119 bool fgcomiq=false;
120 if (popt=="-hval") fillhis=true;
121 else if (popt=="-corciq") fgcomiq=true;
122 string inoutpath = arg[1];
123 string outname = arg[2];
124 int aids[8]={0,1,-1,-1,-1,-1,-1,-1};
125 if (narg>3) sscanf(arg[3],"%d,%d,%d,%d,%d,%d,%d,%d",aids,aids+1,aids+2,aids+3,aids+4,aids+5,aids+6,aids+7);
126 vector<int> adcids;
127 for(int ii=0;ii<8;ii++)
128 if((aids[ii]>=0)&&(aids[ii]<16)) adcids.push_back(aids[ii]);
129 int imin=0;
130 int imax=0;
131 int istep=1;
132 if (narg>4) sscanf(arg[4],"%d,%d,%d",&imin,&imax,&istep);
133 int jf1=0;
134 int jf2=0;
135 int nfreq=0;
136 if (narg>5)
137 sscanf(arg[5],"%d,%d,%d",&jf1,&jf2,&nfreq);
138
139 cout << " ----- corrcrtadc/DecodeProc - Start " << ((fillhis)? " Fill ADC-val histogram" : " Compute Correlation") << endl;
140 cout << " InOutPath= " << inoutpath << " ADCIds: " ;
141 for(size_t ii=0; ii<adcids.size(); ii++) cout << adcids[ii] << " , ";
142 cout << " (NbADC=" << adcids.size() << ") OutFileName= " << outname
143 << " IMin,Max,Step=" << imin << "," << imax << "," << istep << endl;
144 // << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << " ------- " << endl;
145 int rc=0;
146 if (fillhis) rc = ADCFiles2Histo(inoutpath, adcids, outname, imin, imax, istep);
147 else rc=ProcessADCFiles(inoutpath, adcids, outname, fgcomiq, imin, imax, istep, jf1, jf2, nfreq);
148
149 return rc;
150}
151
152// Taille des fichiers 4 voies = 1024*1024*128 ===> 32 x 1024 x 1024 samples / voie
153// On lit par paquet de 4096 = 4 x 1024 / voie --> ADCRDBLKSZ=16x1024=16384
154#define ADCFLEN 4096
155#define ADCRDBLKSZ 16384
156#define ADCNREAD 8192
157
158static int PRTMODULO = 1000;
159
160// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
161/* --Fonction-- */
162int ProcessADCFiles(string& inoutpath, vector<int>& adcids, string& outname, bool fgcombiq, int imin, int imax, int istep, int jf1, int jf2, int nfreq)
163{
164 Timer tm("ProcessADCFiles");
165
166#define NFREQUSED 2040
167#define NFREQUSESTART 5
168#define NFREBIN 4
169
170 sa_size_t NFREQREBIN = NFREQUSED/NFREBIN;
171 sa_size_t NBADC = adcids.size();
172
173 sa_size_t NCHAN = 4*NBADC;
174 TMatrix< complex<r_4> > cfour(NCHAN, NFREQUSED);
175 sa_size_t NCORREL = NCHAN*(NCHAN+1)/2;
176 TMatrix< complex<r_8> > correl(NCORREL,NFREQUSED);
177 TMatrix< complex<r_8> > autocorrel(NCHAN,NFREQUSED);
178 double nsumcor=0.;
179
180 vector<FILE*> fip;
181 vector<int> adcchans;
182 for(size_t ka=0; ka<NBADC; ka++) {
183 fip.push_back(NULL);
184 for(size_t ica=0; ica<4; ica++)
185 adcchans.push_back(adcids[ka]*4+ica);
186 }
187
188 FFTPackServer ffts(false);
189
190 TVector<r_4> sigadc(ADCFLEN);
191 TVector< complex<r_4> > foursig;
192
193 char fname[512];
194
195 signed char rdbuff[ADCRDBLKSZ];
196 for(int ifile=imin; ifile<=imax; ifile+=istep) {
197 for(int ka=0; ka<NBADC; ka++) {
198 sprintf(fname, "%s/adcdumpN%dF%d.bd",inoutpath.c_str(),adcids[ka],ifile);
199 fip[ka]=fopen(fname,"rb");
200 if (fip[ka]==NULL) {
201 cout << "ProcessADCFiles/ERROR opening file " << fname << " AdcId=" << adcids[ka] << " I=" << ifile << endl;
202 throw IOExc(" ADC dump file open error!");
203 }
204 else
205 cout << "ProcessADCFiles/Info - file " << fname << " opened ... " << " AdcId=" << adcids[ka] << " I=" << ifile << endl;
206 }
207
208 for(int m=0; m<ADCNREAD; m++) {
209 sa_size_t irc=0;
210 for(int ka=0; ka<NBADC; ka++) {
211 fread(rdbuff,1,ADCRDBLKSZ,fip[ka]);
212 for(int jac=0; jac<4; jac++) {
213 for(sa_size_t ls=0; ls<ADCFLEN; ls++) sigadc(ls)=(float)rdbuff[ls*4+jac];
214 ffts.FFTForward(sigadc, foursig);
215 cfour.Row(irc)=foursig(Range(NFREQUSESTART,NFREQUSESTART+NFREQUSED-1)).Transpose(); irc++;
216 }
217 }
218 ComputeCorrel(cfour, fgcombiq, correl, autocorrel);
219 nsumcor++;
220 if (m%PRTMODULO==0) cout << " ProcessADCFiles/Info Done read+FFT+correl m=" << m << " / Max=" << ADCNREAD << endl;
221 }
222 for(int ka=0; ka<NBADC; ka++) fclose(fip[ka]);
223 }
224
225
226 correl /= nsumcor;
227 autocorrel /= nsumcor;
228
229 TMatrix< complex<r_8> > correbin(NCORREL,NFREQREBIN);
230 sa_size_t icc=0;
231 for(sa_size_t i=0; i<correbin.NCols(); i++) {
232 for(sa_size_t j=0; j<NFREBIN; j++) {
233 correbin.Column(i) += correl.Column(icc);
234 icc++;
235 }
236 }
237
238 sprintf(fname, "%s/%s",inoutpath.c_str(),outname.c_str());
239 cout << "ProcessADCFiles: Opening file " << fname << " for writing correlation matrix" << endl;
240 POutPersist po(fname);
241 po << PPFNameTag("correlmtx") << correl;
242 po << PPFNameTag("acorrmtx") << autocorrel;
243 po << PPFNameTag("rebincormtx") << correbin;
244
245 cout << "ProcessADCFiles: finished FFT + correlation computing " << endl;
246
247 cout << " --------------- Channel/Correlation number association --------------- " << endl;
248 cout << " CorreNumber- axb (MtxRow) : (Ca, Cb) -- ADCChan(c,d) -- I/Q(e,f)" << endl;
249 cout << " ----------------------------------------------------------------------- " << endl;
250 sa_size_t jcr=0;
251 char IQid[2];
252 char EWid[2];
253 int EWnum[2];
254 int Bid[2];
255 int Inpid[2];
256
257 for(sa_size_t j1=0; j1<cfour.NRows(); j1++) {
258 for(sa_size_t j2=j1; j2<cfour.NRows(); j2++) {
259 cout << jcr << " : (" << j1 << "," << j2 << ") -> ADCChans (" << adcchans[j1] << "," << adcchans[j2] << ") -- (";
260 Bid[0]=adcchans[j1]/4;
261 Bid[1]=adcchans[j2]/4;
262 Inpid[0]=adcchans[j1]%4;
263 Inpid[1]=adcchans[j2]%4;
264
265 if ((adcchans[j1]/4)%2==0) IQid[0]='I'; else IQid[0]='Q';
266 if ((adcchans[j2]/4)%2==0) IQid[1]='I'; else IQid[1]='Q';
267 if (adcchans[j1]%2==0) EWid[0]='E'; else EWid[0]='W';
268 if (adcchans[j2]%2==0) EWid[1]='E'; else EWid[1]='W';
269 cout << IQid[0] << '-' << EWid[0] << (Bid[0]/2)*2+Inpid[0]/2 << ","
270 << IQid[1] << '-' << EWid[1] << (Bid[1]/2)*2+Inpid[1]/2 << ")" << endl;
271 jcr++;
272 }
273 }
274 cout << " ----------------------------------------------------------------------- " << endl;
275
276 return 0;
277}
278
279/* --Fonction-- */
280int ComputeCorrel(TMatrix< complex<r_4> >& cfour, bool fgcombiq,
281 TMatrix< complex<r_8> >& correl, TMatrix< complex<r_8> >& autocorrel)
282{
283 if (fgcombiq) CombineIQPairs(cfour);
284 sa_size_t jcr=0;
285 for(sa_size_t j1=0; j1<cfour.NRows(); j1++) {
286 for(sa_size_t j2=j1; j2<cfour.NRows(); j2++) {
287 for(sa_size_t i=0; i<cfour.NCols(); i++)
288 correl(jcr,i) += complex<r_8> ( cfour(j1,i)*conj(cfour(j2,i)) );
289 if (j1==j2) autocorrel.Row(j1)=correl.Row(jcr);
290 jcr++;
291 }
292 }
293 return 0;
294}
295
296/* --Fonction-- */
297int CombineIQPairs(TMatrix< complex<r_4> >& cfour)
298{
299 // We assume the following channel numbering for the I-Q pairs
300 // Row pairs : (0,4) (1,5) (2,6) (3,7) (8,12) (9,13) (10,14) (11,15) ... are assumed to be
301 // IQ pairs in the matrix cfour
302 // We compute the I rows by the left side band, and the Q rows by the right side band
303
304 // Left = a I + b Q
305 // Right = c I + d Q
306 complex<r_4> za(1.,0.);
307 complex<r_4> zb(0.,1.);
308 complex<r_4> zc(1.,0.);
309 complex<r_4> zd(0.,-1.);
310
311 int nboards = cfour.NRows()/4;
312 for(int jb=0; jb<nboards-1; jb++) {
313 for(int c=0; c<4; c++) {
314 sa_size_t ii = jb*4+c; // sa_size_t are long integers
315 sa_size_t iq = ii+4;
316 // Combine the two I,Q fourier coefficient vectors to form the left and right sid bands
317 TVector< complex<r_4> > leftband = za*cfour.Row(ii)+zb*cfour.Row(iq);
318 TVector< complex<r_4> > rightband = zc*cfour.Row(ii)+zd*cfour.Row(iq);
319 // replace the two I,Q vectors by left and right computed side bands
320 cfour.Row(ii)=leftband;
321 cfour.Row(iq)=rightband;
322 }
323 }
324 return 0;
325}
326
327/* --Fonction-- */
328int ADCFiles2Histo(string& inoutpath, vector<int>& adcids, string& outname, int imin, int imax, int istep)
329{
330 Timer tm("ADCFiles2Histo");
331
332 sa_size_t NBADC = adcids.size();
333
334 sa_size_t NCHAN = 4*NBADC;
335 vector<Histo *> vhist;
336 for(int ka=0; ka<NCHAN; ka++)
337 vhist.push_back(new Histo(-128.5,128.5,257) );
338
339 vector<FILE*> fip;
340 for(size_t ka=0; ka<NBADC; ka++) fip.push_back(NULL);
341
342 char fname[512];
343
344 signed char rdbuff[ADCRDBLKSZ];
345 for(int ifile=imin; ifile<=imax; ifile+=istep) {
346 for(int ka=0; ka<NBADC; ka++) {
347 sprintf(fname, "%s/adcdumpN%dF%d.bd",inoutpath.c_str(),adcids[ka],ifile);
348 fip[ka]=fopen(fname,"rb");
349 if (fip[ka]==NULL) {
350 cout << "ADCFiles2Histo/ERROR opening file " << fname << " AdcId=" << adcids[ka] << " I=" << ifile << endl;
351 throw IOExc(" ADC dump file open error!");
352 }
353 else
354 cout << "ADCFiles2Histo/Info - file " << fname << " opened ... " << " AdcId=" << adcids[ka] << " I=" << ifile << endl;
355 }
356
357 for(int m=0; m<ADCNREAD; m++) {
358 sa_size_t irc=0;
359 for(int ka=0; ka<NBADC; ka++) {
360 fread(rdbuff,1,ADCRDBLKSZ,fip[ka]);
361 for(int jac=0; jac<4; jac++) {
362 for(sa_size_t ls=0; ls<ADCFLEN; ls++) vhist[irc]->Add((r_8)rdbuff[ls*4+jac]);
363 irc++;
364 }
365 }
366 if (m%PRTMODULO==0) cout << " ADCFiles2Histo/Info Done FillHist m=" << m << " / Max=" << ADCNREAD << endl;
367 }
368 for(int ka=0; ka<NBADC; ka++) fclose(fip[ka]);
369 }
370
371 sprintf(fname, "%s/%s",inoutpath.c_str(),outname.c_str());
372 cout << "ADCFiles2Histo: Opening file " << fname << " for writing ADC value histograms " << endl;
373 POutPersist po(fname);
374 char nametag[64];
375 for(size_t ka=0; ka<NCHAN; ka++) {
376 sprintf(nametag,"adcvalC%d",ka);
377 po << PPFNameTag(nametag) << *(vhist[ka]);
378 }
379 for(size_t ka=0; ka<NCHAN; ka++) delete vhist[ka];
380
381 cout << "ADCFiles2Histo: finished ADC value histogram fill " << endl;
382 return 0;
383}
Note: See TracBrowser for help on using the repository browser.