44 FilterID(1), Frequency(60.), Window(0.0), Stride(1.), nFirst(1), nLast(0),
45 nStep(1), nScan(20), nBand(5), nSubs(1), fBand(0.45), nLPF(-1), nWave(16),
46 clean(false), badData(false), noScan(false), nRIF(6),
SNR(2.), reFine(true),
47 dumpStart(0), FilterState(0), SeedFrequency(60.), CurrentTime(Time(0)),
48 StartTime(Time(0)),
Sample(0.0)
71 if(f < 0.)
clean =
true;
72 nSubs = (nT > 0) ? nT : 1;
95 if (!ts.refDVect()->F_data() &&
96 !ts.refDVect()->S_data() &&
97 !ts.refDVect()->S_data() )
98 throw invalid_argument(
"Only float or short data accepted");
102 throw invalid_argument(
"Wrong frequency");
112 unsigned int imax = (
nLPF>0) ?
int(L/2)+1 :
int(L/4)+1;
113 if(
nFirst > imax) std::cout<<
"LineFilter: Invalid harmonic number.\n";
116 if(imax > (
unsigned)L/2) imax = L/2;
122 if(!(&ts))
return false;
126 catch (exception&
e) {
127 cerr <<
"Data error in LineFilter: " << e.what() << endl;
163 nBand = (nB>=2) ? nB : 2;
212 reFine = (sN>0) ?
true :
false;
215 noScan = (f < 0.) ? true :
false;
233 int n =
int(m/(nb*L))*
L;
245 cout <<
" LineFilter::getPSD error: time series is too short to contain\n" 246 <<
" one cycle of fundamental harmonic " <<
Frequency <<
"\n";
250 c = (nb>1) ? 1./(nb-1) : 1.;
256 for(k = 0; k <
nSubs; k++){
259 for(j = 0; j < nb; j++){
261 psd.
data[0] -= tds.
Stack(TD, n, n*j+m*k);
272 td2.
cpf(tds); td2.
cpf(tds,L,0,L);
278 for (i = 2; i < L-1; i+=2){
280 psd.
data[i/2] += c*(a*a+b*b);
301 cout <<
" LineFilter::MakeFilter() error: badData flag is on\n";
311 cout <<
" LineFilter::MakeFilter() error: data length too short to contain\n" 312 <<
" at least one cycle of target frequency = " <<
Frequency <<
" Hz\n";
317 unsigned int imax =
maxLine(L);
365 cout <<
" getLine() error: invalid interference frequency" 373 unsigned int imax =
maxLine(L);
375 if (n/L == 0 || L < 4) {
376 cout <<
" getLine() error: input data length too short to contain\n" 377 <<
" one cycle of target frequency = " <<
Frequency <<
" Hz\n";
399 td2.
cpf(tds); td2.
cpf(tds,L,0,L);
412 for (
unsigned int i = 0;
i < (unsigned)L-1;
i+=2) {
418 amp.
data[
i] += (a*a + b*b)/nSubs;
427 a -= (a>0.) ? long(a+0.5) : long(a-0.5);
433 if((L&1) == 1) tds.
data[L-1] = 0.;
445 tds.getStatistics(a,b);
448 iF = (k == nSubs-1) ? TD.
size() : n*k+
n;
450 for (
int i = 0; i <
L; i++)
451 for(
int j = n*k+i;
j < iF;
j +=
L)
452 TD.
data[
j] = tds.data[i];
462 for (
unsigned int i = nFirst;
i < imax;
i += abs(
nStep)) {
492 register double C,
S,
c,
s, ch,sh,
x;
495 register double* p = TD.
data;
504 cout <<
" getLine() error: invalid interference frequency" 538 for (mode = nFirst; mode < imax; mode += abs(
nStep)) {
545 ch = cos(2*
PI/
double(n));
546 sh = sin(2*
PI/
double(n));
548 for (k = 0; k <
nSubs; k++) {
555 for (i = 0; i <
n; i++) {
564 else if(
FilterID==-2 || mode>nFirst){
567 for (i = 0; i <
n; i++) {
581 for (i = 0; i <
n; i++) {
588 amp.
data[k+
m] = 2.*sqrt(a*a + b*b)/double(n);
589 phi.
data[k+
m] = atan2(b,a);
605 for (mode = nFirst; mode < imax; mode += abs(
nStep)) {
612 for (k = 0; k <
nSubs; k++) {
631 for (i = 0; i < n/2; i++){
638 for (k = 0; k < nSubs-1; k++) {
640 x = phi.
data[k+
m] +omega*(k*n+n/2);
643 x = phi.
data[k+m+1]+omega*(k*n+n/2);
646 for (i = 0; i <
n; i++) {
648 b = amp.
data[k+1+
m]* u;
655 *(p++) -= (a*(n-i) + b*
i)/n;
659 a = amp.
data[nSubs-1+
m];
660 x = phi.
data[nSubs-1+
m]+omega*(TD.
size()-n/2);
663 for (i = TD.
size()-n/2; i<
N; i++){
703 cout <<
" getOmega() error: invalid interference frequency" 721 unsigned int imax =
maxLine(L);
723 if (n/L == 0 || L < 4) {
724 cout <<
" getOmega() error: input data length too short to contain\n" 725 <<
" one cycle of target frequency = " <<
Frequency <<
" Hz\n";
736 double step = n/TDR.
rate();
739 double FSNR =
SNR/(1.+
SNR);
741 for (
int k = 0;
k < nsub;
k++) {
747 td2.
cpf(tds); td2.
cpf(tds,L,0,L);
758 for (
unsigned int i = 2;
i < (unsigned)L-1;
i+=2) {
764 amp.
data[
i] += a*a+b*b;
766 phase +=
axb(wt/2.,
double(i/2));
767 phase -=
intw(phase);
770 F = phase - phi.
data[i+1];
772 phi.
data[
i] += (long(wt*(i/2)+0.5) +
F)/step/(i/2);
791 if(a < 0.0001) a = 0.0001;
793 omega += a*phi.
data[2*
i];
802 omega = (weight>1.) ? omega/weight/(nsub-1) : -
Frequency;
816 #define np 6 // number of points in interpolation scheme in resample() 831 double am = 1., ac = 0.;
854 if (TD.
rate() <= 0.) {
855 cout <<
" fScan() error: invalid sampling rate = " 856 << TD.
rate() <<
" Aborting calculation.\n";
862 cout <<
" fScan() error: invalid interference frequency = " 863 << ff <<
" Aborting calculation.\n";
878 cout <<
" Scanning frequency from "<< ff-mScan*fw/2. <<
" Hz to " 879 << ff+mScan*fw/2. <<
" Hz\n";
894 if (ip > 0 && ip < (mScan-1) && !
badData) {
896 fp += (d > 0.) ? 0.5*fw*(sw.
data[ip + 1] - sw.
data[ip - 1])/d : 0.;
912 for (
int m = 0;
m < 3;
m++) {
931 d=2.*ss[1]-(ss[2]+ss[0]);
934 ac = 0.5 * (ss[2] - ss[0]);
935 fwhm = sqrt(ac*ac+2.*ss[1]*d)/d;
940 ac = (ss[2] > ss[0])? am : -am;
947 if(
fabs(ac) < am) mode = 1;
948 if(
fabs(ac)<am/4. && fw/dfft>0.1) mode = 2;
952 delta = (fc-fp)/fw + ac;
960 if(
fabs(delta) < e1)
break;
961 if(
fabs(delta*fwhm) < e1 && fw/dfft<0.1)
break;
966 ss[0] = ss[1]; ss[1] = ss[2]; ks[2] = 1;
969 ss[2] = ss[1]; ss[1] = ss[0]; ks[0] = 1;
971 fc += (ac>0.) ? fw : -fw;
977 ss[0] = ss[1]; ks[1] = 1;
980 ss[2] = ss[1]; ks[1] = 1;
983 fc += (ac>0.) ? fw : -fw;
987 ks[0] = 1; ks[2] = 1;
989 if(fw/dfft<0.01) fw = 0.01*dfft;
1027 if ( TD.
rate() <= 0. || omega <= 0. ) {
1028 cout <<
" Interference() error: invalid interference frequency = " 1029 << omega <<
"\n Aborting calculation.\n";
1087 int n = ts.getNSample();
1105 TSeries
r(ts.getStartTime(),
Sample, ts.getNSample());
1107 float *vFloat = (
float*)r.refData();
1108 for(
int j = 0;
j<
n;
j++) vFloat[
j] =
float(wts.
data[
j]);
1126 if (!ts.
size())
return;
1127 if (!ts.
rate())
return;
1148 int nTS = ts.
size();
1152 NN = (nTS >> LPF) << LPF;
1171 cout <<
" LineFilter::apply() error: invalid time window "<<
Window<<
" sec.\n";
1180 while(i <= n-nn && nn > 0) {
1183 delta *= double(n - i)/nn;
1190 if(LPF) tx->
cpf(tw,nn,i);
1191 else tx->
cpf(ts,nn,i);
1198 if(omega < 0.) omega =
fScan(*tx);
1206 if(LPF) tw.
sub(*tx,nn,0,i);
1207 else ts.
sub(*tx,nn,0,i);
1225 if(nTS != (
int)ts.
size())
1226 cout <<
"LineFilter::apply(): is "<<ts.
size()<<
", should be: "<<nTS<<
"\n";
1245 list<lineData>::iterator it;
1256 out.
rate((step>0.) ? 1./step : 1.);
1258 if(l_size < 2) {
return out; }
1266 if(m > 0 && m <= v_size)
1272 for(
int i=0;
i<l_size;
i++) {
1276 F *= (m == 0) ? 1. : (v->
first+m-1);
1280 if(n) averF /= double(n);
1284 for(
int i=0;
i<l_size;
i++) {
1294 if(m == 0 || m > v_size)
1302 if(m == 0 || m > v_size){
1304 out.
data[
i] = v_size;
1313 a -= (a>0.) ? long(a+0.5) : long(a-0.5);
1318 if(c ==
'f') out.
data[
i] =
F;
1320 if(c ==
'f' && step<1.1*
Stride){
1322 a -= (a>0.) ? long(a+0.5) : long(a-0.5);
1323 out.
data[
i] = (long(
F*step+0.5) +
a)/step;
1329 out.
data[
i] = (m == 0) ? a : a*(v->
first+m-1);
1333 if(m == 0 || m > v_size)
1334 out.
data[
i] = v_size;
1341 if(m > v_size || v_size < 1)
break;
1343 for(
int j=0;
j<v_size;
j++){
1356 if(m > v_size || v_size < 1)
break;
1358 for(
int j=0;
j<v_size;
j++)
1366 if(m > v_size || v_size < 1)
break;
1368 for(
int j=0;
j<v_size;
j++)
1377 if(m > v_size || v_size < 1)
break;
1379 for(
int j=0;
j<v_size;
j++)
1388 if(m > v_size || v_size < 1)
break;
1390 for(
int j=0;
j<v_size;
j++){
1400 out.
data[
i] /= (a>0.) ?
a : 1.;
1405 if(m > v_size || v_size < 1)
break;
1429 list<lineData>::iterator it =
lineList.begin();
1437 size_t max_size = 0;
1438 for(
size_t i=0;
i<l_size;
i++) {
1441 if(v_size > max_size) max_size = v_size;
1444 int m = (5*max_size+4);
1445 size_t n = m*(l_size+1);
1446 if(n < 4)
return false;
1451 out->
data[0] = max_size;
1452 out->
data[1] = l_size;
1463 double gpsTime = it->T_current.totalS();
1465 int gps =
int(gpsTime)/1000;
1466 out->
data[4] = float(gps);
1467 out->
data[5] = gpsTime - 1000.*double(gps);
1470 for(
size_t i=1;
i<=l_size;
i++) {
1479 for(
unsigned int j=0;
j<max_size;
j++){
1488 out->
data[
i*m+4+
j*5] = 0.;
1489 out->
data[
i*m+5+
j*5] = 0.;
1490 out->
data[
i*m+6+
j*5] = 0.;
1491 out->
data[
i*m+7+
j*5] = 0.;
1492 out->
data[
i*m+8+
j*5] = 0.;
1515 double gpsTime = 0.;
1516 unsigned int count = 0;
1519 if(in.
size() < 6)
return false;
1521 while(count < in.
size()){
1522 p = &(in.
data[count]);
1523 int m =
int(*(p+2) + 0.5);
1524 if(m <= 1)
return false;
1526 int l_size =
int(*(p+1) + 0.5);
1527 if(l_size < 1)
return false;
1529 int v_size =
int(*p + 0.5);
1530 count +=
int(*(p+3)+0.5);
1532 gpsTime = double(
int(*(p+4)+0.5))*1000. + *(p+5);
1539 v.
line.resize(v_size);
1540 v.
noise.resize(v_size);
1549 for(
int i=1;
i<=l_size;
i++) {
1551 if(*p != 0.) step = *p - last;
1554 v.
T_current = Time(
size_t(time+gpsTime));
1559 for(
int j=0;
j<v_size;
j++){
void reset()
Clear/release the internal History vector and reset the current time.
bool LoadTrend(const char *)
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.
void dataCheck(const TSeries &ts) const
Check the data for validity.
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<< "\ layers : "<< nLAYERS<< "\ dF(hz) : "<< dF<< "\ 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<< "SNR "<< xsnr<< endl;wavearray< double > f
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)
virtual size_t size() const
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
printf("total live time: non-zero lags = %10.1f \, liveTot)
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)
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
bool isDataValid(const TSeries &ts) const
Check the data for validity.
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)