26 FilterID(1), Frequency(60.), Window(0.0), Stride(1.), nFirst(1), nLast(0),
27 nStep(1), nScan(20), nBand(5), nSubs(1), fBand(0.45), nLPF(-1), nWave(16),
28 clean(false), badData(false), noScan(false), nRIF(6),
SNR(2.), reFine(true),
29 dumpStart(0), FilterState(0), SeedFrequency(60.), CurrentTime(Time(0)),
30 StartTime(Time(0)),
Sample(0.0)
42 FilterID(1), Frequency(60.), Window(0.0), Stride(1.), nFirst(1), nLast(0),
43 nStep(1), nScan(20), nBand(5), nSubs(1), fBand(0.45), nLPF(-1), nWave(16),
44 clean(false), badData(false), noScan(false), nRIF(6),
SNR(2.), reFine(true),
45 dumpStart(0), FilterState(0), SeedFrequency(60.), CurrentTime(Time(0)),
46 StartTime(Time(0)),
Sample(0.0)
53 if(f < 0.)
clean =
true;
54 nSubs = (nT > 0) ? nT : 1;
60 : FilterID(x.FilterID), Frequency(x.Frequency), Window(x.Window),
61 Stride(x.Stride), nFirst(x.nFirst), nLast(x.nLast), nStep(x.nStep),
62 nScan(x.nScan), nBand(x.nBand), nSubs(x.nSubs), fBand(x.fBand),
63 nLPF(x.nLPF), nWave(x.nWave), clean(x.clean), badData(false),
64 noScan(x.noScan), nRIF(x.nRIF),
SNR(x.
SNR), reFine(x.reFine),
65 dumpStart(0), FilterState(0), SeedFrequency(60.), CurrentTime(Time(0)),
66 StartTime(Time(0)),
Sample(0.)
77 if (!ts.refDVect()->F_data() &&
78 !ts.refDVect()->S_data() &&
79 !ts.refDVect()->S_data() )
80 throw invalid_argument(
"Only float or short data accepted");
84 throw invalid_argument(
"Wrong frequency");
94 unsigned int imax = (
nLPF>0) ?
int(L/2)+1 :
int(L/4)+1;
95 if(
nFirst > imax) std::cout<<
"LineFilter: Invalid harmonic number.\n";
98 if(imax > (
unsigned)L/2) imax = L/2;
104 if(!(&ts))
return false;
108 catch (exception&
e) {
109 cerr <<
"Data error in LineFilter: " << e.what() << endl;
145 nBand = (nB>=2) ? nB : 2;
194 reFine = (sN>0) ?
true :
false;
197 noScan = (f < 0.) ?
true :
false;
215 int n =
int(m/(nb*L))*
L;
227 cout <<
" LineFilter::getPSD error: time series is too short to contain\n"
228 <<
" one cycle of fundamental harmonic " <<
Frequency <<
"\n";
232 c = (nb>1) ? 1./(nb-1) : 1.;
238 for(k = 0; k <
nSubs; k++){
241 for(j = 0; j < nb; j++){
243 psd.
data[0] -= tds.
Stack(TD, n, n*j+m*k);
254 td2.
cpf(tds); td2.
cpf(tds,L,0,L);
260 for (i = 2; i < L-1; i+=2){
262 psd.
data[i/2] += c*(a*a+b*b);
283 cout <<
" LineFilter::MakeFilter() error: badData flag is on\n";
293 cout <<
" LineFilter::MakeFilter() error: data length too short to contain\n"
294 <<
" at least one cycle of target frequency = " <<
Frequency <<
" Hz\n";
299 unsigned int imax =
maxLine(L);
347 cout <<
" getLine() error: invalid interference frequency"
355 unsigned int imax =
maxLine(L);
357 if (n/L == 0 || L < 4) {
358 cout <<
" getLine() error: input data length too short to contain\n"
359 <<
" one cycle of target frequency = " <<
Frequency <<
" Hz\n";
381 td2.
cpf(tds); td2.
cpf(tds,L,0,L);
394 for (
unsigned int i = 0;
i < (unsigned)L-1;
i+=2) {
400 amp.
data[
i] += (a*a + b*b)/nSubs;
409 a -= (a>0.) ? long(a+0.5) : long(a-0.5);
415 if((L&1) == 1) tds.
data[L-1] = 0.;
427 tds.getStatistics(a,b);
430 iF = (k == nSubs-1) ? TD.
size() : n*k+
n;
432 for (
int i = 0; i <
L; i++)
433 for(
int j = n*k+i;
j < iF;
j +=
L)
434 TD.
data[
j] = tds.data[i];
444 for (
unsigned int i = nFirst;
i < imax;
i += abs(
nStep)) {
474 register double C,
S,
c,
s, ch,sh,
x;
477 register double* p = TD.
data;
486 cout <<
" getLine() error: invalid interference frequency"
520 for (mode = nFirst; mode < imax; mode += abs(
nStep)) {
527 ch = cos(2*
PI/
double(n));
528 sh = sin(2*
PI/
double(n));
530 for (k = 0; k <
nSubs; k++) {
537 for (i = 0; i <
n; i++) {
546 else if(
FilterID==-2 || mode>nFirst){
549 for (i = 0; i <
n; i++) {
563 for (i = 0; i <
n; i++) {
570 amp.
data[k+
m] = 2.*sqrt(a*a + b*b)/double(n);
571 phi.
data[k+
m] = atan2(b,a);
587 for (mode = nFirst; mode < imax; mode += abs(
nStep)) {
594 for (k = 0; k <
nSubs; k++) {
613 for (i = 0; i < n/2; i++){
620 for (k = 0; k < nSubs-1; k++) {
622 x = phi.
data[k+
m] +omega*(k*n+n/2);
625 x = phi.
data[k+m+1]+omega*(k*n+n/2);
628 for (i = 0; i <
n; i++) {
630 b = amp.
data[k+1+
m]* u;
637 *(p++) -= (a*(n-i) + b*
i)/n;
641 a = amp.
data[nSubs-1+
m];
642 x = phi.
data[nSubs-1+
m]+omega*(TD.
size()-n/2);
645 for (i = TD.
size()-n/2; i<
N; i++){
685 cout <<
" getOmega() error: invalid interference frequency"
703 unsigned int imax =
maxLine(L);
705 if (n/L == 0 || L < 4) {
706 cout <<
" getOmega() error: input data length too short to contain\n"
707 <<
" one cycle of target frequency = " <<
Frequency <<
" Hz\n";
718 double step = n/TDR.
rate();
721 double FSNR =
SNR/(1.+
SNR);
723 for (
int k = 0;
k < nsub;
k++) {
729 td2.
cpf(tds); td2.
cpf(tds,L,0,L);
740 for (
unsigned int i = 2;
i < (unsigned)L-1;
i+=2) {
746 amp.
data[
i] += a*a+b*b;
748 phase +=
axb(wt/2.,
double(i/2));
749 phase -=
intw(phase);
752 F = phase - phi.
data[i+1];
754 phi.
data[
i] += (long(wt*(i/2)+0.5) +
F)/step/(i/2);
773 if(a < 0.0001) a = 0.0001;
775 omega += a*phi.
data[2*
i];
784 omega = (weight>1.) ? omega/weight/(nsub-1) : -
Frequency;
798 #define np 6 // number of points in interpolation scheme in resample()
813 double am = 1., ac = 0.;
836 if (TD.
rate() <= 0.) {
837 cout <<
" fScan() error: invalid sampling rate = "
838 << TD.
rate() <<
" Aborting calculation.\n";
844 cout <<
" fScan() error: invalid interference frequency = "
845 << ff <<
" Aborting calculation.\n";
860 cout <<
" Scanning frequency from "<< ff-mScan*fw/2. <<
" Hz to "
861 << ff+mScan*fw/2. <<
" Hz\n";
876 if (ip > 0 && ip < (mScan-1) && !
badData) {
878 fp += (d > 0.) ? 0.5*fw*(sw.
data[ip + 1] - sw.
data[ip - 1])/d : 0.;
894 for (
int m = 0;
m < 3;
m++) {
913 d=2.*ss[1]-(ss[2]+ss[0]);
916 ac = 0.5 * (ss[2] - ss[0]);
917 fwhm = sqrt(ac*ac+2.*ss[1]*d)/d;
922 ac = (ss[2] > ss[0])? am : -am;
929 if(
fabs(ac) < am) mode = 1;
930 if(
fabs(ac)<am/4. && fw/dfft>0.1) mode = 2;
934 delta = (fc-fp)/fw + ac;
942 if(
fabs(delta) < e1)
break;
943 if(
fabs(delta*fwhm) < e1 && fw/dfft<0.1)
break;
948 ss[0] = ss[1]; ss[1] = ss[2]; ks[2] = 1;
951 ss[2] = ss[1]; ss[1] = ss[0]; ks[0] = 1;
953 fc += (ac>0.) ? fw : -fw;
959 ss[0] = ss[1]; ks[1] = 1;
962 ss[2] = ss[1]; ks[1] = 1;
965 fc += (ac>0.) ? fw : -fw;
969 ks[0] = 1; ks[2] = 1;
971 if(fw/dfft<0.01) fw = 0.01*dfft;
1009 if ( TD.
rate() <= 0. || omega <= 0. ) {
1010 cout <<
" Interference() error: invalid interference frequency = "
1011 << omega <<
"\n Aborting calculation.\n";
1069 int n = ts.getNSample();
1087 TSeries
r(ts.getStartTime(),
Sample, ts.getNSample());
1089 float *vFloat = (
float*)r.refData();
1090 for(
int j = 0;
j<
n;
j++) vFloat[
j] =
float(wts.
data[
j]);
1108 if (!ts.
size())
return;
1109 if (!ts.
rate())
return;
1130 int nTS = ts.
size();
1134 NN = (nTS >> LPF) << LPF;
1153 cout <<
" LineFilter::apply() error: invalid time window "<<
Window<<
" sec.\n";
1162 while(i <= n-nn && nn > 0) {
1165 delta *= double(n - i)/nn;
1172 if(LPF) tx->
cpf(tw,nn,i);
1173 else tx->
cpf(ts,nn,i);
1180 if(omega < 0.) omega =
fScan(*tx);
1188 if(LPF) tw.
sub(*tx,nn,0,i);
1189 else ts.
sub(*tx,nn,0,i);
1207 if(nTS != (
int)ts.
size())
1208 cout <<
"LineFilter::apply(): is "<<ts.
size()<<
", should be: "<<nTS<<
"\n";
1227 list<lineData>::iterator it;
1238 out.
rate((step>0.) ? 1./step : 1.);
1240 if(l_size < 2) {
return out; }
1248 if(m > 0 && m <= v_size)
1254 for(
int i=0;
i<l_size;
i++) {
1258 F *= (m == 0) ? 1. : (v->
first+m-1);
1262 if(n) averF /= double(n);
1266 for(
int i=0;
i<l_size;
i++) {
1276 if(m == 0 || m > v_size)
1284 if(m == 0 || m > v_size){
1286 out.
data[
i] = v_size;
1295 a -= (a>0.) ? long(a+0.5) : long(a-0.5);
1300 if(c ==
'f') out.
data[
i] =
F;
1302 if(c ==
'f' && step<1.1*
Stride){
1304 a -= (a>0.) ? long(a+0.5) : long(a-0.5);
1305 out.
data[
i] = (long(
F*step+0.5) +
a)/step;
1311 out.
data[
i] = (m == 0) ? a : a*(v->
first+m-1);
1315 if(m == 0 || m > v_size)
1316 out.
data[
i] = v_size;
1323 if(m > v_size || v_size < 1)
break;
1325 for(
int j=0;
j<v_size;
j++){
1338 if(m > v_size || v_size < 1)
break;
1340 for(
int j=0;
j<v_size;
j++)
1348 if(m > v_size || v_size < 1)
break;
1350 for(
int j=0;
j<v_size;
j++)
1359 if(m > v_size || v_size < 1)
break;
1361 for(
int j=0;
j<v_size;
j++)
1370 if(m > v_size || v_size < 1)
break;
1372 for(
int j=0;
j<v_size;
j++){
1382 out.
data[
i] /= (a>0.) ? a : 1.;
1387 if(m > v_size || v_size < 1)
break;
1411 list<lineData>::iterator it =
lineList.begin();
1419 size_t max_size = 0;
1420 for(
size_t i=0;
i<l_size;
i++) {
1423 if(v_size > max_size) max_size = v_size;
1426 int m = (5*max_size+4);
1427 size_t n = m*(l_size+1);
1428 if(n < 4)
return false;
1433 out->
data[0] = max_size;
1434 out->
data[1] = l_size;
1445 double gpsTime = it->T_current.totalS();
1447 int gps =
int(gpsTime)/1000;
1448 out->
data[4] = float(gps);
1449 out->
data[5] = gpsTime - 1000.*double(gps);
1452 for(
size_t i=1;
i<=l_size;
i++) {
1461 for(
unsigned int j=0;
j<max_size;
j++){
1470 out->
data[
i*m+4+
j*5] = 0.;
1471 out->
data[
i*m+5+
j*5] = 0.;
1472 out->
data[
i*m+6+
j*5] = 0.;
1473 out->
data[
i*m+7+
j*5] = 0.;
1474 out->
data[
i*m+8+
j*5] = 0.;
1497 double gpsTime = 0.;
1498 unsigned int count = 0;
1501 if(in.
size() < 6)
return false;
1503 while(count < in.
size()){
1505 int m =
int(*(p+2) + 0.5);
1506 if(m <= 1)
return false;
1508 int l_size =
int(*(p+1) + 0.5);
1509 if(l_size < 1)
return false;
1511 int v_size =
int(*p + 0.5);
1512 count +=
int(*(p+3)+0.5);
1514 gpsTime = double(
int(*(p+4)+0.5))*1000. + *(p+5);
1521 v.
line.resize(v_size);
1522 v.
noise.resize(v_size);
1531 for(
int i=1;
i<=l_size;
i++) {
1533 if(*p != 0.) step = *p - last;
1536 v.
T_current = Time(
size_t(time+gpsTime));
1541 for(
int j=0;
j<v_size;
j++){
virtual size_t size() const
void reset()
Clear/release the internal History vector and reset the current time.
bool LoadTrend(const char *)
printf("total live time: non-zero lags = %10.1f \n", liveTot)
virtual void rate(double r)
The linefilter class containes methods to track and remove quasi- monochromatic lines.
double axb(double, double)
wavearray< double > a(hp.size())
LineFilter * clone(void) const
Clone a LineFilter.
double Interference(WaveData &, double)
std::vector< f_complex > amplitude
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\t layers : "<< nLAYERS<< "\t dF(hz) : "<< dF<< "\t dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1)*itime+ifreq;double time=itime *dT;double freq=(ifreq >0)?ifreq *dF:dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
wavearray< double > psd(33)
virtual void ReadBinary(const char *, int=0)
bool DumpTrend(const char *, int=0)
~LineFilter(void)
Destroy the LineFilter object and release the function storage.
wavearray< wavereal > WaveData
virtual void start(double s)
cout<< "Selected Pixels : "<< nPix<< endl;wc.cluster(1, 1);SSeries< double > ss
void setFilter(int nF=1, int nL=0, int nS=1, int nD=-1, int nB=5, int nR=6, int nW=8)
std::list< lineData > lineList
double getOmega(const WaveData &, int=2)
int Sample(Complex *H, double *t, int N, double totalMass, double Rate, wavearray< double > &hps, wavearray< double > &hxs)
void sub(const wavearray< DataType_t > &, int=0, int=0, int=0)
LineFilter(void)
Build an empty LineFilter.
lineData getHeteroLine(WaveData &)
unsigned int maxLine(int)
std::vector< float > filter
int getLayer(wavearray< DataType_t > &w, double n)
param: n - layer number
void dataCheck(const TSeries &ts) const
Check the data for validity.
lineData getLine(WaveData &)
TSeries operator()(const TSeries &ts)
WaveData getPSD(const WaveData &, int=1)
void setFScan(double f=0., double sn=2., double fS=0.45, int nS=20)
wavearray< float > getTrend(int, char)
bool isDataValid(const TSeries &ts) const
Check the data for validity.
virtual void DumpBinary(const char *, int=0)
double fScan(const WaveData &)
double makeFilter(const WaveData &, int=0)
wavearray< double > ts(N)
virtual void resample(double, int=6)
double fabs(const Complex &x)
void resample(const wavearray< DataType_t > &, double, int=6)
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
Meyer< double > S(1024, 2)
The LineFilter class containes methods to track and remove quasi- monochromatic lines.
double Stack(const wavearray< DataType_t > &, int)
std::complex< float > f_complex
std::vector< float > noise
virtual void resize(unsigned int)
void Inverse(int n=-1)
param: n - number of steps (-1 means full reconstruction)
std::vector< float > line
void putLayer(wavearray< DataType_t > &, double n)
param: n - layer number
TSeries apply(const TSeries &ts)
The argument time series is filtered to remove lines, and the argument TSeries ts is left unchanged...
void cpf(const wavearray< DataType_t > &, int=0, int=0, int=0)