diff --git a/StationSim/src/ArrayConfig.cc b/StationSim/src/ArrayConfig.cc
index b1369a3a8a5267218f22d6b4f440a4a9b7ba6739..afab72bd72c2c46d55e797500020c3a69f7f0389 100644
--- a/StationSim/src/ArrayConfig.cc
+++ b/StationSim/src/ArrayConfig.cc
@@ -25,6 +25,8 @@
ArrayConfig::ArrayConfig (string config_file)
{
ifstream s (config_file.c_str (), ifstream::in);
+
+ conf_file = config_file;
int n;
diff --git a/StationSim/src/ArrayConfig.h b/StationSim/src/ArrayConfig.h
index 5f93a5d2bdc88c8571b7659a1bf8c9cfe727c9cd..81521a0fc72945396deda86bcc390c031d065809 100644
--- a/StationSim/src/ArrayConfig.h
+++ b/StationSim/src/ArrayConfig.h
@@ -36,6 +36,7 @@ class ArrayConfig
public:
ArrayConfig (string config_file);
+ string conf_file;
// returns the coordinates
const LoVec_double& getPointX () const
diff --git a/StationSim/src/MDL.cc b/StationSim/src/MDL.cc
index 4245416e3583f65c263c959036c0d5681d95f876..ebf54d9ee8829bcb75717221a42eeee22fd50b27 100644
--- a/StationSim/src/MDL.cc
+++ b/StationSim/src/MDL.cc
@@ -19,9 +19,6 @@
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
-// Chris Broekema, november 2002.
-
-
/**************************************************************************/
/* This file implements the Minimun Description Length algorithm. This is */
@@ -31,84 +28,52 @@
#include <MDL.h>
-/* TODO: Check if the cast to real needs to be done earlier */
-using namespace blitz;
+//using namespace blitz;
-unsigned int mdl (LoMat_double d_mat, unsigned int nantl, unsigned int nsnsh)
+int mdl (const LoVec_double& input, int M, int N)
{
- double nom;
- double denom;
- double mdlmin;
- int nnsnsh=-nsnsh;
- unsigned int kmin = 0;
- LoVec_double d;
- LoVec_double MDL(d_mat.cols());
-
- nsnsh=1;
-
- if (d_mat.cols() == d_mat.rows()) {
- d.resize(d_mat.cols());
-
- d = LCSMath::diag( d_mat );
- d = LCSMath::sort( d );
- d.reverseSelf(firstDim);
- for ( int k = d.lbound(firstDim); k < d.ubound(firstDim); k++ ) {
-
- denom = (1 /(double)(d.ubound(firstDim)-k)) *
- sum(d( Range( k+1, d.ubound(firstDim) ) ) ) ;
- nom = product(d(Range(k+1,d.ubound(firstDim)))) ;
- MDL(k) = nnsnsh * log(nom/(pow(denom,(1/(nantl-(k+1)))))) +
- (k+1)/2 * (2*nantl-(k+1)) * log(nsnsh);
- if ( k == 0 ) {
- // init
- kmin = k+1;
- mdlmin = MDL(k);
- } else {
- if (MDL(k) < mdlmin ) {
- kmin = k+1;
- mdlmin = MDL(k);
+ double nom;
+ double denom;
+ double mdlmin;
+ int min = 0;
+
+ LoVec_double MDL (input.size ());
+ LoVec_double d = input.copy();
+ //d = LCSMath::sort (d);
+ d.reverseSelf (blitz::firstDim);
+
+ for (int Ns = 0; Ns < M; Ns++)
+ {
+ denom = 0;
+ for (int i = Ns; i < M; i++) {
+ denom += d (i);
+ }
+ denom *= (1 / (double) (M - Ns));
+ denom = pow (denom, (double) (M - Ns));
+
+ // cout << denom << endl;
+
+ nom = 1;
+ for (int i = Ns; i < M; i++) {
+ nom *= d (i);
+ }
+
+ // cout << nom << endl;
+
+ MDL (Ns) = -1 * N * log (nom / denom) + Ns / 2 * (2 * M - Ns) * log (N);
+
+ cout << MDL (Ns) << endl;
+
+ if (Ns == 0) {
+ mdlmin = MDL(Ns);
+ } else {
+ if (MDL (Ns) < mdlmin) {
+ min = Ns;
+ mdlmin = MDL (Ns);
+ }
+ }
}
- }
- }
- } else {
- cout << "MDL encountered non-square matrix" << endl;
- }
- return kmin;
-}
-
-unsigned int mdl (LoVec_double d_mat, unsigned int nantl, unsigned int nsnsh)
-{
- double nom;
- double denom;
- double mdlmin;
- unsigned int kmin = 0;
- LoVec_double d;
- LoVec_double MDL(d_mat.size());
-
- d.resize(d_mat.size());
- d = LCSMath::sort(d_mat);
- d.reverseSelf(firstDim);
- nsnsh=1;
-
- int nnsnsh=-nsnsh;
- for ( int k = d.lbound(firstDim); k < d.ubound(firstDim); k++ ) {
-
- denom = (1 /(double)(d.ubound(firstDim)-k)) *
- sum(d( Range( k+1, d.ubound(firstDim) ) ) ) ;
- nom = product(d(Range(k+1,d.ubound(firstDim)))) ;
- MDL(k) = nnsnsh * log(nom/(pow(denom,(1/(nantl-(k+1)))))) +
- (k+1)/2 * (2*nantl-(k+1)) * log(nsnsh);
- if ( k == 0 ) {
- // init
- kmin = k+1;
- mdlmin = MDL(k);
- } else {
- if (MDL(k) < mdlmin ) {
- kmin = k+1;
- mdlmin = MDL(k);
- }
- }
- }
- return kmin;
+ // cout << MDL << endl;
+ return min;
}
diff --git a/StationSim/src/MDL.h b/StationSim/src/MDL.h
index db519b4b86a41f75b65e38af50ad0ebd2bcef85f..df8d126c2738b5067759b2302861de31e915d195 100644
--- a/StationSim/src/MDL.h
+++ b/StationSim/src/MDL.h
@@ -25,7 +25,8 @@
#define STATIONSIM_MDL_H
#include <Math/LCSMath.h>
#include <Common/Lorrays.h>
+#include <cmath>
+
+int mdl (const LoVec_double& d, int M, int N) ;
-unsigned int mdl (LoMat_double d, unsigned int nantl, unsigned int nsnsh) ;
-unsigned int mdl (LoVec_double d, unsigned int nantl, unsigned int nsnsh) ;
#endif
diff --git a/StationSim/src/PhaseShift.cc b/StationSim/src/PhaseShift.cc
index c9aaf5825c13758a8f8b7e137ca0fa7ee917d5ba..08dd4d35ca68ac9089a6f143a5dbc66a34cc7bff 100644
--- a/StationSim/src/PhaseShift.cc
+++ b/StationSim/src/PhaseShift.cc
@@ -105,7 +105,7 @@ namespace PhaseShift
dprintf1 (1) ("foward fft done\n");
// compute vector of doas (per each antenna)
- LoVec_double doa = DOA (antennas.getPointX (), antennas.getPointY (), theta, phi);
+ LoVec_double doa = DOA (antennas.getPointX (), antennas.getPointY (), theta, phi, nfft);
// this is the output matrix: one column of signal per each antenna
LoMat_dcomplex phased_signal (antennas.size (), siglen);
@@ -140,13 +140,13 @@ namespace PhaseShift
return fs;
}
- LoVec_double DOA (const LoVec_double& px, const LoVec_double& py, double theta, double phi)
+ LoVec_double DOA (const LoVec_double& px, const LoVec_double& py, double theta, double phi, int nfft)
{
FailWhen1 (px.size () != py.size (), "vector size mismatch");
LoVec_double res (px.size ());
res = -2 * M_PI * (px * sin (theta) * cos (phi) + py * sin (theta) * sin (phi));
- return res;
+ return res / nfft;
}
}; // namespace PhaseShift
diff --git a/StationSim/src/PhaseShift.h b/StationSim/src/PhaseShift.h
index 86fccc7f002874c28a4a307c16cea9594b03a906..9728cc9d7b97816b36411e1e5e1072d82a10108c 100644
--- a/StationSim/src/PhaseShift.h
+++ b/StationSim/src/PhaseShift.h
@@ -74,7 +74,8 @@ namespace PhaseShift
LoVec_double DOA (const LoVec_double& px,
const LoVec_double& py,
double theta,
- double phi);
+ double phi,
+ int nfft);
}; // namespace PhaseShift
diff --git a/StationSim/src/StationSim.cc b/StationSim/src/StationSim.cc
index 5b1dcfdde3bae0d7d573caf6693063132276c454..330525da214bd1f4bb79f6c9a6518af8933224ad 100644
--- a/StationSim/src/StationSim.cc
+++ b/StationSim/src/StationSim.cc
@@ -97,6 +97,8 @@ void StationSim::define (const ParamBlock& params)
const int nbeam = params.getInt ("nbeam", 1);
const int maxNtarget = params.getInt ("maxntarget", 1);
const int maxNrfi = params.getInt ("maxnrfi", 1);
+ const int STArate = params.getInt ("starate", 100);
+ const int WDrate = params.getInt ("wdrate", 10000);
const int delayMod = modulationWindowSize - 1;
const int delayPhase = nfft_phaseshift - 1;
const int delaySubFilt = nrcu * (nsubband - 1);
@@ -191,7 +193,7 @@ void StationSim::define (const ParamBlock& params)
}
simul.addStep (add_signals);
- // Create the NoEMI data readers
+// // Create the NoEMI data readers
// for (int i = 0; i < nrcu; ++i) {
// sprintf (suffix, "%d", i);
@@ -224,7 +226,8 @@ void StationSim::define (const ParamBlock& params)
sprintf (suffix, "%d", i);
// The beamformer object
- Step beam (WH_BeamFormer(suffix, nrcu + 1, 1, nrcu, nbeam, maxNtarget, maxNrfi),
+ Step beam (WH_BeamFormer(suffix, nrcu + 1, 1, nrcu, nbeam, maxNtarget, maxNrfi,
+ STArate),
string("beam_former_") + suffix, false);
for (int j = 0; j < nrcu; ++j) {
beam.getInData (j).setReadDelay (delaySubFilt);
@@ -234,6 +237,11 @@ void StationSim::define (const ParamBlock& params)
}
beam.getInData (nrcu).setReadDelay (delaySubFilt + delayBeamForm);
beam.setRate(nsubband);
+ if (STArate <= WDrate) {
+ beam.setInRate (STArate, nrcu);
+ } else {
+ beam.setInRate (WDrate, nrcu);
+ }
simul.addStep (beam);
@@ -248,30 +256,60 @@ void StationSim::define (const ParamBlock& params)
sta.getOutData (j).setWriteDelay (delaySubFilt);
}
sta.setRate(nsubband);
+ sta.setOutRate(nsubband + STArate);
simul.addStep (sta);
// the Weight Determination Object
- Step weight_det (WH_WeightDetermination(suffix, 0, 1, nrcu, bfDipoleFile),
+ Step weight_det (WH_WeightDetermination(suffix, 0, 1, nrcu, bfDipoleFile, 0.3, 0.3),
string("weight_det_") + suffix, false);
weight_det.getOutData (0).setWriteDelay (delaySubFilt);
- weight_det.setRate(nsubband);
+ weight_det.setRate(nsubband + WDrate);
simul.addStep(weight_det);
-
+
+// // Deterministic nulls should have the same rate as STA !!!!!!!
+// // add a deterministic null
+// Step det_null (WH_WeightDetermination(suffix, 0, 1, nrcu, bfDipoleFile, 2.583, 0.8238),
+// string("det_null_") + suffix, false);
+
+// det_null.getOutData (0).setWriteDelay (delaySubFilt);
+// det_null.setRate(nsubband + STArate);
+// simul.addStep(det_null);
+
+
+// // add a another deterministic null
+// Step det_null_2 (WH_WeightDetermination(suffix, 0, 1, nrcu, bfDipoleFile, -2.9709, 0.7128),
+// string("det_null2_") + suffix, false);
+
+// det_null_2.getOutData (0).setWriteDelay (delaySubFilt);
+// det_null_2.setRate(nsubband + STArate);
+// simul.addStep(det_null_2);
+
// the projection object
- Step projection (WH_Projection(suffix, 3, 1, nrcu, maxNrfi),
+ int noutProj = 3;
+ Step projection (WH_Projection(suffix, noutProj, 1, nrcu, maxNrfi),
string("projection_") + suffix, false);
- for (int j = 0; j < 3; ++j) { // One input + mdl + eigvectors
+ for (int j = 0; j < noutProj; ++j) { // One input + two det null + mdl + eigvectors
projection.getInData (j).setReadDelay (delaySubFilt);
}
for (int j = 0; j < 1; ++j) {
projection.getOutData (j).setWriteDelay (delaySubFilt);
}
-
projection.setRate(nsubband);
+ for (int i = 0; i < noutProj - 2; ++i) {
+ projection.setInRate (nsubband + WDrate);
+ }
+ for (int i = noutProj - 2; i < noutProj; ++i) {
+ projection.setInRate(nsubband + STArate, i);
+ }
+ if (STArate <= WDrate) {
+ projection.setOutRate (STArate);
+ } else {
+ projection.setOutRate (WDrate);
+ }
simul.addStep(projection);
}
@@ -343,11 +381,20 @@ void StationSim::define (const ParamBlock& params)
string ("projection_") + suffix2 + string (".in_mdl"));
// connect the Weight determinator to Projection
- simul.connect (string ("weight_det_") + suffix2 + string (".out"),
+ simul.connect (string ("weight_det_") + suffix2 + string (".out_0"),
string ("projection_") + suffix2 + string (".in_0"));
+
+// // connect the deterministic null to the projection
+// simul.connect (string ("det_null_") + suffix2 + string (".out_0"),
+// string ("projection_") + suffix2 + string(".in_1"));
+
+// // connect the deterministic null to the projection
+// simul.connect (string ("det_null2_") + suffix2 + string (".out_0"),
+// string ("projection_") + suffix2 + string(".in_2"));
+
}
- // Connect the data_reader to the subband filterbank
+// // Connect the data_reader to the subband filterbank
// for (int i = 0; i < nrcu; ++i) {
// sprintf (suffix, "%d", i);
diff --git a/StationSim/src/WH_AddSignals.cc b/StationSim/src/WH_AddSignals.cc
index 48e6b82c9923e5491424e2c277a7b833375068e4..8b79e5f7c4e9036d5ca8d18b1afc08b36267cbd4 100644
--- a/StationSim/src/WH_AddSignals.cc
+++ b/StationSim/src/WH_AddSignals.cc
@@ -59,8 +59,8 @@ WH_AddSignals::WH_AddSignals (const string& name,
}
// DEBUG
itsCount = 0;
- itsFileOutReal.open ("/home/alex/gerdes/add_real.txt");
- itsFileOutComplex.open ("/home/alex/gerdes/add_complex.txt");
+ itsFileOutReal.open ("/home/alex/data/add_real.txt");
+ itsFileOutComplex.open ("/home/alex/data/add_complex.txt");
itsFileOutReal.precision(20);
itsFileOutComplex.precision(20);
}
@@ -112,6 +112,8 @@ void WH_AddSignals::process ()
}
itsFileOutReal << endl;
itsFileOutComplex << endl;
+
+ //cout << itsCount++ << endl;
}
}
diff --git a/StationSim/src/WH_BandSep.cc b/StationSim/src/WH_BandSep.cc
index 7a7d1e2ab3a3af418a42dca7af6b2780cb7b7803..a6d006a4525e0e71d91d4d1b040f1c1104d34e9a 100644
--- a/StationSim/src/WH_BandSep.cc
+++ b/StationSim/src/WH_BandSep.cc
@@ -56,9 +56,9 @@ WH_BandSep::WH_BandSep (const string& name, unsigned int nsubband,
}
// DEBUG
- handle = gnuplot_init();
- itsFileOutReal.open ((string ("/home/alex/gerdes/subband_real_") + name + string (".txt")).c_str ());
- itsFileOutComplex.open ((string ("/home/alex/gerdes/subband_complex_") + name + string (".txt")).c_str ());
+// handle = gnuplot_init();
+// itsFileOutReal.open ((string ("/home/alex/gerdes/subband_real_") + name + string (".txt")).c_str ());
+// itsFileOutComplex.open ((string ("/home/alex/gerdes/subband_complex_") + name + string (".txt")).c_str ());
}
WH_BandSep::~WH_BandSep()
@@ -70,9 +70,9 @@ WH_BandSep::~WH_BandSep()
delete [] itsOutHolders;
// DEBUG
- gnuplot_close (handle);
- itsFileOutReal.close ();
- itsFileOutComplex.close ();
+// gnuplot_close (handle);
+// itsFileOutReal.close ();
+// itsFileOutComplex.close ();
}
WorkHolder* WH_BandSep::construct(const string& name, int ninput, int noutput, const ParamBlock& params)
@@ -120,12 +120,12 @@ void WH_BandSep::process()
}
// DEBUG
- for (int i = 0; i < itsNsubband; i++) {
- itsFileOutReal << real (itsOutHolders[i]->getBuffer ()[0]) << " ";
- itsFileOutComplex << imag (itsOutHolders[i]->getBuffer ()[0]) << " ";
- }
- itsFileOutReal << endl;
- itsFileOutComplex << endl;
+// for (int i = 0; i < itsNsubband; i++) {
+// itsFileOutReal << real (itsOutHolders[i]->getBuffer ()[0]) << " ";
+// itsFileOutComplex << imag (itsOutHolders[i]->getBuffer ()[0]) << " ";
+// }
+// itsFileOutReal << endl;
+// itsFileOutComplex << endl;
}
}
}
diff --git a/StationSim/src/WH_BeamFormer.cc b/StationSim/src/WH_BeamFormer.cc
index 8d8efb90b9cdd41dc2b265aaaa812edfcbb00a1f..67038b534676818cb9496e29c51c0437f8381163 100644
--- a/StationSim/src/WH_BeamFormer.cc
+++ b/StationSim/src/WH_BeamFormer.cc
@@ -20,7 +20,6 @@
//#
//# $Id$
//#
-//# Chris Broekema, january 2003
//#
@@ -39,7 +38,9 @@ WH_BeamFormer::WH_BeamFormer (const string& name,
unsigned int nout,
unsigned int nrcu,
unsigned int nbeam,
- unsigned int maxNtarget, unsigned int maxNrfi)
+ unsigned int maxNtarget,
+ unsigned int maxNrfi,
+ int delay)
: WorkHolder (nin, nout, name,"WH_BeamFormer"), // Check number of inputs and outputs
itsInHolders (0),
itsOutHolders (0),
@@ -69,11 +70,13 @@ WH_BeamFormer::WH_BeamFormer (const string& name,
}
//DEBUG
-// itsFileOutReal.open ((string ("/home/alex/gerdes/BF_real_") + name + string (".txt")).c_str ());
-// itsFileOutComplex.open ((string ("/home/alex/gerdes/BF_complex_") + name + string (".txt")).c_str ());
-// itsFileInput.open ("/home/alex/gerdes/testvector.txt");
-// itsTestVector.resize (itsNrcu);
+ itsFileOutReal.open ((string ("/home/alex/gerdes/BF_real_") + name + string (".txt")).c_str ());
+ itsFileOutComplex.open ((string ("/home/alex/gerdes/BF_complex_") + name + string (".txt")).c_str ());
+// itsFileInput.open ("/home/alex/temp/test_vectorEssai.txt");
+ itsPos = 0;
+// itsTestVector.resize(itsNrcu, 25445);
// itsFileInput >> itsTestVector;
+ itsDelay = delay;
}
@@ -87,9 +90,10 @@ WH_BeamFormer::~WH_BeamFormer()
delete itsOutHolders[i];
}
delete [] itsOutHolders;
-// // DEBUG
-// itsFileOutReal.close ();
-// itsFileOutComplex.close ();
+
+ // DEBUG
+ itsFileOutReal.close ();
+ itsFileOutComplex.close ();
// itsFileInput.close ();
}
@@ -101,23 +105,28 @@ WorkHolder* WH_BeamFormer::construct (const string& name,
params.getInt ("nrcu", 10),
params.getInt ("nbeam", 10),
params.getInt ("maxntarget", 10),
- params.getInt ("maxnrfi", 10));
+ params.getInt ("maxnrfi", 10),
+ 10);
}
WH_BeamFormer* WH_BeamFormer::make (const string& name) const
{
return new WH_BeamFormer (name, getInputs(), getOutputs(),
- itsNrcu, itsNbeam, itsMaxNtarget, itsMaxNrfi);
+ itsNrcu, itsNbeam, itsMaxNtarget, itsMaxNrfi, itsDelay);
}
void WH_BeamFormer::process()
{
if (getOutputs () > 0) {
+ dcomplex temp;
for (int i = 0; i < itsNrcu; i++) {
sample(i) = (dcomplex)itsInHolders[i]->getBuffer()[0];
+
+// //DEBUG
+// sample(i) = itsTestVector(i, itsPos);
}
- LoVec_dcomplex weight(itsWeight.getBuffer(), itsNrcu, duplicateData);
+ LoVec_dcomplex weight(itsWeight.getBuffer(), itsNrcu, duplicateData);
dcomplex output = 0;
for (int i = 0; i < itsNrcu; i++) {
@@ -129,11 +138,16 @@ void WH_BeamFormer::process()
itsOutHolders[j]->getBuffer ()[0] = output;
}
- // // DEBUG
- // cout << sample << " " << real(sample(0)) << " " << imag(sample(0)) << endl;
- // itsFileOutReal << real(sample) << endl;
- // itsFileOutComplex << imag (sample) << " " << endl;
- // cout << sample << endl;
+ // DEBUG
+// if (itsDelay == 0) {
+ itsPos++;
+// itsPos %= itsTestVector.cols();
+ itsFileOutReal << real(output) << endl;
+ itsFileOutComplex << imag (output) << endl;
+// } else {
+// itsDelay--;
+// }
+ cout << itsPos << endl;
}
}
diff --git a/StationSim/src/WH_BeamFormer.h b/StationSim/src/WH_BeamFormer.h
index fef2bc793dbe58dfe01fa7602dd53c4b82ef46ff..577b44b7e36b67673f319e402822d89e3bf386f0 100644
--- a/StationSim/src/WH_BeamFormer.h
+++ b/StationSim/src/WH_BeamFormer.h
@@ -46,7 +46,8 @@ public:
/// The first WorkHolder should have nin=0.
WH_BeamFormer (const string& name,
unsigned int nin, unsigned int nout, unsigned int nrcu,
- unsigned int nbeam, unsigned int maxNtarget, unsigned int maxNrfi);
+ unsigned int nbeam, unsigned int maxNtarget, unsigned int maxNrfi,
+ int delay);
virtual ~WH_BeamFormer();
@@ -91,10 +92,12 @@ private:
LoVec_dcomplex sample; // current sample in Blitz format
//DEBUG
+ int itsPos;
ofstream itsFileOutReal;
ofstream itsFileOutComplex;
ifstream itsFileInput;
- LoVec_double itsTestVector;
+ LoMat_dcomplex itsTestVector;
+ int itsDelay;
};
#endif
diff --git a/StationSim/src/WH_Projection.cc b/StationSim/src/WH_Projection.cc
index 9e4397364bfd112cea2e6319604f06551e4daefc..2464a4abba1b3ad493ff1a473d126a68861340c1 100644
--- a/StationSim/src/WH_Projection.cc
+++ b/StationSim/src/WH_Projection.cc
@@ -18,7 +18,6 @@
//# along with this program; if not, write to the Free Software
//# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//#
-//# Chris Broekema, january 2003
//#
//# $Id$
//#
@@ -44,7 +43,11 @@ WH_Projection::WH_Projection (const string& name, unsigned int nin, unsigned int
itsRFISources ("in_rfi", nant, maxnrfi),
itsNrcu (nant),
itsMaxRFI (maxnrfi),
- itsWeight (itsNrcu)
+ itsDetectedRFIs (nin - 3),
+ itsWeight (itsNrcu),
+ itsV (0,0),
+ itsA (itsNrcu)
+
{
char str[8];
// nin is the number of steering vectors that will be used as an input (at least one
@@ -83,44 +86,47 @@ WH_Projection* WH_Projection::make (const string& name) const
return new WH_Projection (name, getInputs(), getOutputs(), itsNrcu, itsMaxRFI);
}
+
void WH_Projection::process()
{
if (getOutputs() > 0) {
- // Get the number of detected RFI sources from the STA comp. (internally calc. by MDL)
- int detectedRFIs = (int)(itsNumberOfRFIs.getBuffer()[0]);
- detectedRFIs = 2;
-
- // Read in the eigenvectors from the STA component. They are the detected rfi sources and
- // should be nulled.
- // LoMat_dcomplex V (itsInHolders[0]->getBuffer(), shape(itsNrcu, itsMaxRFI), duplicateData);
- LoMat_dcomplex V (itsRFISources.getBuffer(), shape(itsNrcu, detectedRFIs), duplicateData);
-
- // Get the steering vector form the weight determination comp. assume that only one will be put in.
- LoVec_dcomplex a (itsInHolders[0]->getBuffer(), itsNrcu, duplicateData);
-
- // Calculate the weights using the eigenvectors and the steering vector
- itsWeight = getWeights(V, a);
-
- // DEBUG matrix inverse
-// int N = 3;
-// LoMat_dcomplex Orig(N,N);
-// LoMat_dcomplex A(N, N);
-// LoVec_dcomplex Vec(N);
-
-// Vec = 1, 1, 1;
-// Orig = 1,2,3,
-// 4,5,6,
-// 7,8,9;
-// Orig(0,0)=dcomplex(0,1);
-
-// cout << Orig << endl;
-// A = LCSMath::invert(Orig);
-// cout << A << endl;
-
-// cout << LCSMath::matMult(A, Orig) << endl;
-
-// cout << Orig << endl << Vec << endl << vm_mult(Vec,Orig) << endl;
- // END DEBUG matrix inverse
+ if (itsRFISources.doHandle ()) {
+ // Get the number of detected RFI sources from the STA comp. (internally calc. by MDL)
+ int NumberOfEigenVectors = (int)(itsNumberOfRFIs.getBuffer()[0]);
+ itsDetectedRFIs = NumberOfEigenVectors + getInputs() - 3;
+
+ // Read in the eigenvectors from the STA component. They are the detected rfi sources and
+ // should be nulled. The deterministic nulls should be put in this matrix.
+ LoMat_dcomplex V (itsRFISources.getBuffer(), shape(NumberOfEigenVectors, itsNrcu), duplicateData);
+ itsV.resize(V.shape());
+ itsV = V;
+ itsV.transposeSelf(secondDim, firstDim);
+
+ // Merge the deterministic nulls in the eigenvectors matrix
+ itsV.resizeAndPreserve(itsNrcu, itsDetectedRFIs); // Add number of deterministic nulls
+ for (int i = 1; i < getInputs() - 2; i++) {
+ LoVec_dcomplex det_rfi (itsInHolders[i]->getBuffer(), itsNrcu, duplicateData);
+ itsV(blitz::Range::all(), itsDetectedRFIs - i) = det_rfi;
+ }
+ }
+
+ if (itsInHolders[0]->doHandle ()) {
+ // Get the steering vector form the weight determination comp. assume that only one will be put in.
+ LoVec_dcomplex steerv (itsInHolders[0]->getBuffer(), itsNrcu, duplicateData);
+ itsA = steerv;
+ }
+
+ if (itsRFISources.doHandle () || itsInHolders[0]->doHandle ()) {
+ // Calculate the weights using the eigenvectors and the steering vector
+ if (itsDetectedRFIs == 1) {
+ LoVec_dcomplex Vvec = itsV(blitz::Range::all(), itsV.lbound(secondDim));
+ itsWeight = getWeights(Vvec, itsA);
+ } else if (itsDetectedRFIs > 1) {
+ itsWeight = getWeights(itsV, itsA);
+ } else {
+ itsWeight = itsA;
+ }
+ }
// Copy output to the next step
for (int i = 0; i < getOutputs(); i++) {
@@ -178,22 +184,15 @@ LoVec_dcomplex WH_Projection::getWeights (LoMat_dcomplex V, LoVec_dcomplex a)
LoVec_dcomplex w (V.rows());
w = 1;
- LoMat_dcomplex Pv (V.shape());
LoMat_dcomplex I = LCSMath::diag(w); // Create Identity matrix
- //cout << w << endl << Pv << endl << I << endl << V << endl << a << endl;
-
LoMat_dcomplex temp = LCSMath::matMult(LCSMath::hermitianTranspose(V), V);
- // cout << temp << endl;
- Pv = LCSMath::invert(temp); // (VHV)^-1
- // cout << Pv << endl;
- Pv = LCSMath::matMult(V, LCSMath::matMult(Pv, LCSMath::hermitianTranspose(V))); // VPvVH
- // cout << Pv << endl;
+ temp = LCSMath::invert(temp); // (VHV)^-1
+ LoMat_dcomplex temp2 = LCSMath::matMult(temp, LCSMath::hermitianTranspose(V));
+ LoMat_dcomplex Pv = LCSMath::matMult(V, temp2); // VPvVH
Pv = I - Pv;
- // cout << Pv << endl;
w = mv_mult(Pv, a); // w = Pv a;
- // cout << w << endl;
return w;
}
diff --git a/StationSim/src/WH_Projection.h b/StationSim/src/WH_Projection.h
index d8b43cc55e146d64743082a40e3213706df7b77a..e24d5bb629deb8979d2416c06b88d7f2a02f5b7c 100644
--- a/StationSim/src/WH_Projection.h
+++ b/StationSim/src/WH_Projection.h
@@ -58,7 +58,7 @@ public:
/// Make a fresh copy of the WH object.
virtual WH_Projection* make (const string& name) const;
-
+
/// Do a process step.
virtual void process();
@@ -89,7 +89,10 @@ private:
/// Length of buffers.
unsigned int itsNrcu;
unsigned int itsMaxRFI;
+ int itsDetectedRFIs;
LoVec_dcomplex itsWeight;
+ LoMat_dcomplex itsV;
+ LoVec_dcomplex itsA;
LoVec_dcomplex WH_Projection::getWeights (LoVec_dcomplex B, LoVec_dcomplex d) ;
LoVec_dcomplex WH_Projection::getWeights (LoMat_dcomplex V, LoVec_dcomplex a) ;
diff --git a/StationSim/src/WH_ReadSignal.cc b/StationSim/src/WH_ReadSignal.cc
index cb2443e32184e43baf21b25b54c677c2038db431..db7eb97439263896820be9c048e1e2ef0929b972 100644
--- a/StationSim/src/WH_ReadSignal.cc
+++ b/StationSim/src/WH_ReadSignal.cc
@@ -84,7 +84,14 @@ void WH_ReadSignal::process ()
if (!itsFile.eof ()) {
itsFile >> sample;
} else {
- sample = 0;
+ char s[256];
+ itsFile.close();
+ itsFile.open(itsFileName.c_str());
+
+ while (strncmp (s, "Data :", 6) != 0)
+ itsFile.getline (s, 256);
+
+ itsFile >> sample;
}
for (int i = 0; i < getOutputs (); i++) {
diff --git a/StationSim/src/WH_STA.cc b/StationSim/src/WH_STA.cc
index 33af7946353bedb3724f174e3ce1316388f2c239..db8a8905fb4997d94b2d51e3d77d45e22586d914 100644
--- a/StationSim/src/WH_STA.cc
+++ b/StationSim/src/WH_STA.cc
@@ -18,7 +18,6 @@
//# along with this program; if not, write to the Free Software
//# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//#
-//# Chris Broekema, january 2003.
//#
#include <stdio.h> // for sprintf
@@ -33,18 +32,20 @@
#include <blitz/blitz.h>
-WH_STA::WH_STA (const string& name, unsigned int nin, unsigned int nout,
- unsigned int nant, unsigned int maxnrfi, unsigned int buflength)
+WH_STA::WH_STA (const string& name, int nin, int nout,
+ int nant, int maxnrfi, int buflength)
: WorkHolder (nin, nout, name, "WH_STA"),
itsInHolders (0),
itsOutHolders (0),
itsNumberOfRFIs ("out_mdl", 1, 1),
itsNrcu (nant),
itsMaxRFI (maxnrfi),
- itsBufLength (4),
- itsBuffer (itsNrcu, itsBufLength),
- itsSnapshot (itsNrcu),
- itsCurPos (0)
+ itsBufLength (buflength),
+ itsBuffer (nant, buflength),
+ itsPos (0),
+ itsEvectors (nant,nant),
+ itsEvalues (nant),
+ itsAcm (nant, nant)
{
char str[8];
if (nant > 0) {
@@ -63,6 +64,12 @@ WH_STA::WH_STA (const string& name, unsigned int nin, unsigned int nout,
// RFI sources. This should be implemented more elegantly.
itsOutHolders[i] = new DH_SampleC (string("out_") + str, itsNrcu, itsMaxRFI);
}
+ //DEBUG
+// delay = 3;
+ itsTestVector.resize(itsNrcu, 300);
+ itsFileInput.open ("/home/alex/gerdes/test_vectorSTA_sin.txt");
+ itsFileInput >> itsTestVector;
+ itsCount = 0;
}
WH_STA::~WH_STA()
@@ -93,75 +100,73 @@ void WH_STA::process()
{
using namespace blitz;
+ // Place the next incoming sample vector in the snapshot buffer
+ // Keep a cylic buffer for the input snapshots
for (int i = 0; i < itsNrcu; i++) {
- itsSnapshot(i) = *itsInHolders[i]->getBuffer();
+ //itsBuffer(i, itsPos) = itsInHolders[i]->getBuffer()[0];
+ itsBuffer(i, itsPos) = itsTestVector(i, itsCount);
}
-
- unsigned int ub = itsBuffer.ubound(secondDim);
- unsigned int lb = itsBuffer.lbound(secondDim);
-
- itsBuffer(Range::all(), Range(lb, ub-1)) =
- itsBuffer(Range::all(), Range(lb+1,ub));
- itsBuffer(Range::all(), ub) = itsSnapshot;
-
-
- // This is only a single snapshot.. Make a buffer to really to something
- // useful
-// DH_SampleC::BufferType* bufin = itsInHolders->getBuffer();
-
-// itsBuffer = *bufin;
+ itsPos = (itsPos + 1) % itsBuffer.cols();
+ itsCount = (itsCount + 1) % itsTestVector.cols();
// // Select the appropriate algorithm
// // See if the PASTd algorithm need updating
// // Use either EVD or SVD for updating
- //DEBUG AG: put in a testvector into the ACM calc.
- LoMat_dcomplex testVector(itsNrcu, itsNrcu);
- testVector = 0;
- testVector(0,Range::all()) = 1;
-
- // EVD - first calculate the ACM
- LoMat_dcomplex itsAcm (itsNrcu, itsNrcu) ;
-// itsAcm = LCSMath::acm(itsBuffer);
- itsAcm = LCSMath::acm(testVector);
- // cout << itsAcm << endl;
-
-// // EVD - using the ACM, calculate eigen vectors and values.
-// LoMat_dcomplex itsEvectors;
-// LoVec_double itsEvalues;
-// LCSMath::eig(itsAcm, itsEvectors, itsEvalues);
-
-// // EVD - Determine the number of sources in the signal
-// // TODO: fit MDL into Common library.
-// unsigned int RFI = mdl(itsEvalues, itsNrcu, itsBufLength);
-
-// // EVD - Find the appropriate Eigen vectors
-// // Assume the most powerfull sources are in front. This may not be true!!
-// // TODO check if this is a good assumption -- possible overflow
-// LoMat_dcomplex B(itsNrcu, RFI);
-// // B = itsEvectors(Range::all(), Range(itsEvectors.lbound(firstDim),
-// // itsEvectors.lbound(firstDim) + RFI - 1));
-
-// LoVec_dcomplex w (itsNrcu) ; // the weight vector
-// dcomplex *w_ptr = w.data() ;
-
-// LoVec_dcomplex d(itsNrcu);
- // w = getWeights(B, d);
-
- // Now assign the calculated weight vector to the output
-
-
-// memcpy(itsOutHolders[0]->getBuffer(), d.data(), itsNrcu * sizeof(DH_SampleC::BufferType));
-// memcpy(itsOutHolders[1]->getBuffer(), mdl, sizeof(DH_SampleC::BufferType));
-
-
-// cout << itsBuffer << endl;
-// cout << itsAcm << endl;
+ if (getOutHolder(0)->doHandle ()) {
+ // Create contigeous buffer
+ itsBuffer = CreateContigeousBuffer (itsBuffer, itsPos);
+ itsPos = 0; // Set cyclic buffer pointer to the first entry because it is ordered now
+
+ // calculate the ACM
+ itsAcm = LCSMath::acm (itsBuffer);
+
+ // using the ACM, calculate eigen vectors and values.
+ LCSMath::eig (itsAcm, itsEvectors, itsEvalues);
+
+ // Determine the number of sources in the signal
+
+ // put in a vector with two large eigen values and see what MDL finds
+ itsEvalues = -1;
+// itsEvalues(itsNrcu-1) = 200000000.0;
+// itsEvalues(itsNrcu-2) = 100000000.0;
+// itsEvalues(itsNrcu-3) = 10000.0;
+ // itsEvalues(itsNrcu-4) = 1000.0;
+
+ double RFI = (double) mdl (itsEvalues, itsNrcu, itsBufLength);
+ // cout << RFI << endl;
+ if (RFI > itsMaxRFI) RFI = itsMaxRFI;
+
+ // DEBUG
+// RFI = 1;
+
+// // DEBUG
+// for (int i = 0; i < itsNrcu; i++) {
+// cout << itsEvalues(i) << endl;
+// if (itsEvalues(i) > 10) {
+// for (int j = 0; j < itsNrcu; j++) {
+// cout << itsEvectors(i, j) << endl;
+// }
+// }
+// }
+
+ // Find the appropriate Eigen vectors
+ LoMat_dcomplex B = itsEvectors(Range(itsEvectors.ubound(secondDim) - RFI + 1,
+ itsEvectors.ubound(secondDim)), Range::all());
+
+ // cout << B << endl;
+ // cout << itsEvalues << endl;
+
+ // Now assign the Eigen vectors to the output
+ for (int i = 0; i < getOutputs() - 1; i++) {
+ memcpy(itsOutHolders[i]->getBuffer(), B.data(), itsNrcu * (int)RFI * sizeof(DH_SampleC::BufferType));
+ }
+ memcpy(itsNumberOfRFIs.getBuffer(), &RFI, sizeof(DH_SampleR::BufferType));
+ }
// SVD
- // PASTd
-
+ // PASTd
}
void WH_STA::dump() const
@@ -183,3 +188,19 @@ DataHolder* WH_STA::getOutHolder (int channel)
return &itsNumberOfRFIs;
}
+LoMat_dcomplex WH_STA::CreateContigeousBuffer (const LoMat_dcomplex& aBuffer, int pos)
+{
+ using namespace blitz;
+ AssertStr (pos <= aBuffer.ubound(secondDim), "Can't create contigeous buffer, pos too high!");
+ Range all = Range::all();
+ LoMat_dcomplex output (aBuffer.shape());
+
+ int p = 0;
+ for (int i = pos - 1; i >= 0; i--) {
+ output (all, p++) = aBuffer(all, i);
+ }
+ for (int i = aBuffer.ubound(secondDim); i >= pos; i--) {
+ output (all, p++) = aBuffer(all, i);
+ }
+ return output;
+}
diff --git a/StationSim/src/WH_STA.h b/StationSim/src/WH_STA.h
index e9ddbcf429b64ac110c2cae21826d0baa30b3e6c..ef92408f0bff8ad4fef030a19283d20b3c47907c 100644
--- a/StationSim/src/WH_STA.h
+++ b/StationSim/src/WH_STA.h
@@ -42,11 +42,11 @@ public:
/// are created and how many elements there are in the buffer.
/// The first WorkHolder should have nin=0.
WH_STA (const string& name,
- unsigned int nin,
- unsigned int nout,
- unsigned int nant,
- unsigned int maxnrfi,
- unsigned int buflength
+ int nin,
+ int nout,
+ int nant,
+ int maxnrfi,
+ int buflength
);
virtual ~WH_STA();
@@ -80,7 +80,7 @@ private:
/// Forbid assignment.
WH_STA& operator= (const WH_STA&);
- /// Calculate a steer vector
+ LoMat_dcomplex CreateContigeousBuffer (const LoMat_dcomplex& aBuffer, int pos);
/// In- and OutHolders
DH_SampleC** itsInHolders;
@@ -88,12 +88,18 @@ private:
DH_SampleR itsNumberOfRFIs;
/// Length of buffers.
- unsigned int itsNrcu;
- unsigned int itsMaxRFI;
- unsigned int itsBufLength;
-
+ int itsNrcu;
+ int itsMaxRFI;
+ int itsBufLength;
+
LoMat_dcomplex itsBuffer;
- LoVec_dcomplex itsSnapshot;
- unsigned int itsCurPos;
+ int itsPos;
+ LoMat_dcomplex itsEvectors;
+ LoVec_double itsEvalues;
+ LoMat_dcomplex itsAcm;
+ int delay;
+ ifstream itsFileInput;
+ LoMat_dcomplex itsTestVector;
+ int itsCount;
};
#endif
diff --git a/StationSim/src/WH_WeightDetermination.cc b/StationSim/src/WH_WeightDetermination.cc
index dc040a069668f7935ba4b3b4f762e40e7e98294d..c743067aa39988de78012ffa19f6c3cf123f5a56 100644
--- a/StationSim/src/WH_WeightDetermination.cc
+++ b/StationSim/src/WH_WeightDetermination.cc
@@ -18,7 +18,6 @@
//# along with this program; if not, write to the Free Software
//# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//#
-//# Chris Broekema, januari 2003.
//#
@@ -31,18 +30,22 @@
#include <Math/LCSMath.h>
#include <blitz/blitz.h>
-WH_WeightDetermination::WH_WeightDetermination(const string& name, unsigned int nin, unsigned int nout, unsigned int nant, string s)
- : WorkHolder (nin, nout, name, "WH_WeightDetermination"),
- itsOutHolder (0),
- itsNrcu (0),
- itsArray (0),
- itsArrayFile (s)
+WH_WeightDetermination::WH_WeightDetermination(const string& name, unsigned int nin, unsigned int nout,
+ unsigned int nant, string s, double phi, double theta)
+ : WorkHolder (nin, nout, name, "WH_WeightDetermination"),
+ itsOutHolders (0),
+ itsNrcu (nant),
+ itsArray (s),
+ itsPhi (phi),
+ itsTheta (theta)
{
- itsNrcu = nant;
- itsArray = new ArrayConfig (s);
-
+ char str[8];
if (nout > 0) {
- itsOutHolder = new DH_SampleC("out", itsNrcu, 1);
+ itsOutHolders = new DH_SampleC* [nout];
+ }
+ for (unsigned int i = 0; i < nout; i++) {
+ sprintf (str, "%d", i);
+ itsOutHolders[i] = new DH_SampleC (string("out_") + str, itsNrcu, 1);
}
}
@@ -53,7 +56,7 @@ WH_WeightDetermination::~WH_WeightDetermination()
WH_WeightDetermination* WH_WeightDetermination::make (const string& name) const
{
- return new WH_WeightDetermination (name, getInputs(), getOutputs(), itsNrcu, itsArrayFile);
+ return new WH_WeightDetermination (name, getInputs(), getOutputs(), itsNrcu, itsArray.conf_file, itsPhi, itsTheta);
}
void WH_WeightDetermination::preprocess()
@@ -63,24 +66,23 @@ void WH_WeightDetermination::preprocess()
void WH_WeightDetermination::process()
{
- double phi = 0.1;
- double theta = -0.2;
-
-
- LoVec_dcomplex d(itsNrcu);
- d = steerv(phi, theta, itsArray->getPointX(), itsArray->getPointY());
-
- memcpy(itsOutHolder->getBuffer(), d.data(), itsNrcu * sizeof(DH_SampleC::BufferType));
+ if (getOutputs() > 0) {
+ LoVec_dcomplex d = steerv(itsPhi, itsTheta, itsArray.getPointX(), itsArray.getPointY());
+
+ for (int i = 0; i < getOutputs(); i++) {
+ memcpy(itsOutHolders[i]->getBuffer(), d.data(), itsNrcu * sizeof(DH_SampleC::BufferType));
+ }
+ }
}
void WH_WeightDetermination::dump() const
{
using namespace blitz;
- LoVec_dcomplex weight(itsOutHolder->getBuffer(), itsNrcu, duplicateData);
+ LoVec_dcomplex weight(itsOutHolders[0]->getBuffer(), itsNrcu, duplicateData);
- cout << "Weight vector Buffer: " << endl;
- cout << weight<< endl;
+ // cout << "Weight vector Buffer: " << endl;
+ // cout << weight<< endl;
}
@@ -92,9 +94,8 @@ DataHolder* WH_WeightDetermination::getInHolder (int channel)
DH_SampleC* WH_WeightDetermination::getOutHolder (int channel)
{
- AssertStr (channel < getOutputs(),
- "output channel too high");
- return itsOutHolder;
+ AssertStr (channel < getOutputs(), "output channel too high");
+ return itsOutHolders[channel];
}
LoVec_dcomplex WH_WeightDetermination::steerv (double phi, double theta, LoVec_double px, LoVec_double py) {
diff --git a/StationSim/src/WH_WeightDetermination.h b/StationSim/src/WH_WeightDetermination.h
index d80ede654f87ad67c40ae93cee402ef7f0f1e679..3058e8521939646f36a2cefa568c1c86d4b5c1d4 100644
--- a/StationSim/src/WH_WeightDetermination.h
+++ b/StationSim/src/WH_WeightDetermination.h
@@ -44,10 +44,12 @@ public:
/// are created and how many elements there are in the buffer.
/// The first WorkHolder should have nin=0.
WH_WeightDetermination (const string& name,
- unsigned int nin,
- unsigned int nout,
- unsigned int nant,
- string s);
+ unsigned int nin,
+ unsigned int nout,
+ unsigned int nant,
+ string s,
+ double phi,
+ double theta);
virtual ~WH_WeightDetermination();
@@ -72,8 +74,6 @@ public:
/// Get a pointer to the i-th output DataHolder.
virtual DH_SampleC* getOutHolder (int channel);
- LoVec_dcomplex steerv (double u, double v, LoVec_double px, LoVec_double py) ;
-
private:
/// Forbid copy constructor.
WH_WeightDetermination (const WH_WeightDetermination&);
@@ -81,18 +81,14 @@ private:
/// Forbid assignment.
WH_WeightDetermination& operator= (const WH_WeightDetermination&);
- /// Calculate a steer vector
+ LoVec_dcomplex steerv (double u, double v, LoVec_double px, LoVec_double py) ;
- /// In- and OutHolders
- DH_SampleC* itsOutHolder;
+ // OutHolders
+ DH_SampleC** itsOutHolders;
- /// Length of buffers.
unsigned int itsNrcu;
-
- string itsConfigFile;
-
- vector<float> itsDipoleLoc;
- ArrayConfig *itsArray;
- string itsArrayFile;
+ ArrayConfig itsArray;
+ double itsPhi;
+ double itsTheta;
};
#endif