23 FilterID(1), Frequency(60.), Window(0.0), Stride(1.), nFirst(1), nLast(0),
24 nStep(1), nScan(20), nBand(5), nSubs(1), fBand(0.45), nLPF(-1), nWave(16),
25 clean(false), badData(false), noScan(false), nRIF(6),
SNR(2.), reFine(true),
26 dumpStart(0), FilterState(0), SeedFrequency(60.), CurrentTime(0.),
39 FilterID(1), Frequency(60.), Window(0.0), Stride(1.), nFirst(1), nLast(0),
40 nStep(1), nScan(20), nBand(5), nSubs(1), fBand(0.45), nLPF(-1), nWave(16),
41 clean(false), badData(false), noScan(false), nRIF(6),
SNR(2.), reFine(true),
42 dumpStart(0), FilterState(0), SeedFrequency(60.), CurrentTime(0.),
50 if(f < 0.)
clean =
true;
51 nSubs = (nT > 0) ? nT : 1;
57 : FilterID(x.FilterID), Frequency(x.Frequency), Window(x.Window),
58 Stride(x.Stride), nFirst(x.nFirst), nLast(x.nLast), nStep(x.nStep),
59 nScan(x.nScan), nBand(x.nBand), nSubs(x.nSubs), fBand(x.fBand),
60 nLPF(x.nLPF), nWave(x.nWave), clean(x.clean), badData(false),
61 noScan(x.noScan), nRIF(x.nRIF),
SNR(x.
SNR), reFine(x.reFine),
62 dumpStart(0), FilterState(0), SeedFrequency(60.), CurrentTime(0.),
73 unsigned int imax = (
nLPF>0) ?
int(L/2)+1 :
int(L/4)+1;
74 if(
nFirst > imax) std::cout<<
"linefilter: Invalid harmonic number.\n";
77 if(imax > (
unsigned)L/2) imax = L/2;
112 nBand = (nB>=2) ? nB : 2;
161 reFine = (sN>0) ?
true :
false;
164 noScan = (f < 0.) ?
true :
false;
182 int n =
int(m/(nb*L))*
L;
194 cout <<
" linefilter::getPSD error: time series is too short to contain\n"
195 <<
" one cycle of fundamental harmonic " <<
Frequency <<
"\n";
199 c = (nb>1) ? 1./(nb-1) : 1.;
205 for(k = 0; k <
nSubs; k++){
208 for(j = 0; j < nb; j++){
210 psd.
data[0] -= tds.
Stack(TD, n, n*j+m*k);
221 td2.
cpf(tds); td2.
cpf(tds,L,0,L);
227 for (i = 2; i < L-1; i+=2){
229 psd.
data[i/2] += c*(a*a+b*b);
250 cout <<
" linefilter::MakeFilter() error: badData flag is on\n";
260 cout <<
" linefilter::MakeFilter() error: data length too short to contain\n"
261 <<
" at least one cycle of target frequency = " <<
Frequency <<
" Hz\n";
266 unsigned int imax =
maxLine(L);
314 cout <<
" getLine() error: invalid interference frequency"
322 unsigned int imax =
maxLine(L);
324 if (n/L == 0 || L < 4) {
325 cout <<
" getLine() error: input data length too short to contain\n"
326 <<
" one cycle of target frequency = " <<
Frequency <<
" Hz\n";
348 td2.
cpf(tds); td2.
cpf(tds,L,0,L);
361 for (
unsigned int i = 0;
i < (unsigned)L-1;
i+=2) {
367 amp.
data[
i] += (a*a + b*b)/nSubs;
376 a -= (a>0.) ? long(a+0.5) : long(a-0.5);
382 if((L&1) == 1) tds.
data[L-1] = 0.;
394 tds.getStatistics(a,b);
397 iF = (k == nSubs-1) ? TD.
size() : n*k+
n;
399 for (
int i = 0; i <
L; i++)
400 for(
int j = n*k+i;
j < iF;
j +=
L)
401 TD.
data[
j] = tds.data[i];
411 for (
unsigned int i = nFirst;
i < imax;
i += abs(
nStep)) {
441 register double C,
S,
c,
s, ch,sh,
x;
444 register double* p = TD.
data;
453 cout <<
" getLine() error: invalid interference frequency"
487 for (mode = nFirst; mode < imax; mode += abs(
nStep)) {
494 ch = cos(2*
PI/
double(n));
495 sh = sin(2*
PI/
double(n));
497 for (k = 0; k <
nSubs; k++) {
504 for (i = 0; i <
n; i++) {
513 else if(
FilterID==-2 || mode>nFirst){
516 for (i = 0; i <
n; i++) {
530 for (i = 0; i <
n; i++) {
537 amp.
data[k+
m] = 2.*sqrt(a*a + b*b)/double(n);
538 phi.
data[k+
m] = atan2(b,a);
554 for (mode = nFirst; mode < imax; mode += abs(
nStep)) {
561 for (k = 0; k <
nSubs; k++) {
580 for (i = 0; i < n/2; i++){
587 for (k = 0; k < nSubs-1; k++) {
589 x = phi.
data[k+
m] +omega*(k*n+n/2);
592 x = phi.
data[k+m+1]+omega*(k*n+n/2);
595 for (i = 0; i <
n; i++) {
597 b = amp.
data[k+1+
m]* u;
604 *(p++) -= (a*(n-i) + b*
i)/n;
608 a = amp.
data[nSubs-1+
m];
609 x = phi.
data[nSubs-1+
m]+omega*(TD.
size()-n/2);
612 for (i = TD.
size()-n/2; i<
N; i++){
652 cout <<
" getOmega() error: invalid interference frequency"
670 unsigned int imax =
maxLine(L);
672 if (n/L == 0 || L < 4) {
673 cout <<
" getOmega() error: input data length too short to contain\n"
674 <<
" one cycle of target frequency = " <<
Frequency <<
" Hz\n";
685 double step = n/TDR.
rate();
688 double FSNR =
SNR/(1.+
SNR);
690 for (
int k = 0;
k < nsub;
k++) {
696 td2.
cpf(tds); td2.
cpf(tds,L,0,L);
707 for (
unsigned int i = 2;
i < (unsigned)L-1;
i+=2) {
713 amp.
data[
i] += a*a+b*b;
715 phase +=
axb(wt/2.,
double(i/2));
716 phase -=
intw(phase);
719 F = phase - phi.
data[i+1];
721 phi.
data[
i] += (long(wt*(i/2)+0.5) +
F)/step/(i/2);
740 if(a < 0.0001) a = 0.0001;
742 omega += a*phi.
data[2*
i];
751 omega = (weight>1.) ? omega/weight/(nsub-1) : -
Frequency;
765 #define np 6 // number of points in interpolation scheme in resample()
780 double am = 1., ac = 0.;
803 if (TD.
rate() <= 0.) {
804 cout <<
" fScan() error: invalid sampling rate = "
805 << TD.
rate() <<
" Aborting calculation.\n";
811 cout <<
" fScan() error: invalid interference frequency = "
812 << ff <<
" Aborting calculation.\n";
827 cout <<
" Scanning frequency from "<< ff-mScan*fw/2. <<
" Hz to "
828 << ff+mScan*fw/2. <<
" Hz\n";
843 if (ip > 0 && ip < (mScan-1) && !
badData) {
845 fp += (d > 0.) ? 0.5*fw*(sw.
data[ip + 1] - sw.
data[ip - 1])/d : 0.;
861 for (
int m = 0;
m < 3;
m++) {
880 d=2.*ss[1]-(ss[2]+ss[0]);
883 ac = 0.5 * (ss[2] - ss[0]);
884 fwhm = sqrt(ac*ac+2.*ss[1]*d)/d;
889 ac = (ss[2] > ss[0])? am : -am;
896 if(
fabs(ac) < am) mode = 1;
897 if(
fabs(ac)<am/4. && fw/dfft>0.1) mode = 2;
901 delta = (fc-fp)/fw + ac;
909 if(
fabs(delta) < e1)
break;
910 if(
fabs(delta*fwhm) < e1 && fw/dfft<0.1)
break;
915 ss[0] = ss[1]; ss[1] = ss[2]; ks[2] = 1;
918 ss[2] = ss[1]; ss[1] = ss[0]; ks[0] = 1;
920 fc += (ac>0.) ? fw : -fw;
926 ss[0] = ss[1]; ks[1] = 1;
929 ss[2] = ss[1]; ks[1] = 1;
932 fc += (ac>0.) ? fw : -fw;
936 ks[0] = 1; ks[2] = 1;
938 if(fw/dfft<0.01) fw = 0.01*dfft;
976 if ( TD.
rate() <= 0. || omega <= 0. ) {
977 cout <<
" Interference() error: invalid interference frequency = "
978 << omega <<
"\n Aborting calculation.\n";
1039 if (!ts.
size())
return;
1040 if (!ts.
rate())
return;
1061 int nTS = ts.
size();
1065 NN = (nTS >> LPF) << LPF;
1084 cout <<
" linefilter::apply() error: invalid time window "<<
Window<<
" sec.\n";
1093 while(i <= n-nn && nn > 0) {
1096 delta *= double(n - i)/nn;
1103 if(LPF) tx->
cpf(tw,nn,i);
1104 else tx->
cpf(ts,nn,i);
1111 if(omega < 0.) omega =
fScan(*tx);
1119 if(LPF) tw.
sub(*tx,nn,0,i);
1120 else ts.
cpf(*tx,nn,0,i);
1138 if(nTS != (
int)ts.
size())
1139 cout <<
"linefilter::apply(): is "<<ts.
size()<<
", should be: "<<nTS<<
"\n";
1158 list<linedata>::iterator it;
1169 out.
rate((step>0.) ? 1./step : 1.);
1171 if(l_size < 2) {
return out; }
1179 if(m > 0 && m <= v_size)
1185 for(
int i=0;
i<l_size;
i++) {
1189 F *= (m == 0) ? 1. : (v->
first+m-1);
1193 if(n) averF /= double(n);
1197 for(
int i=0;
i<l_size;
i++) {
1207 if(m == 0 || m > v_size)
1215 if(m == 0 || m > v_size){
1217 out.
data[
i] = v_size;
1226 a -= (a>0.) ? long(a+0.5) : long(a-0.5);
1231 if(c ==
'f') out.
data[
i] =
F;
1233 if(c ==
'f' && step<1.1*
Stride){
1235 a -= (a>0.) ? long(a+0.5) : long(a-0.5);
1236 out.
data[
i] = (long(
F*step+0.5) +
a)/step;
1242 out.
data[
i] = (m == 0) ? a : a*(v->
first+m-1);
1246 if(m == 0 || m > v_size)
1247 out.
data[
i] = v_size;
1254 if(m > v_size || v_size < 1)
break;
1256 for(
int j=0;
j<v_size;
j++){
1269 if(m > v_size || v_size < 1)
break;
1271 for(
int j=0;
j<v_size;
j++)
1279 if(m > v_size || v_size < 1)
break;
1281 for(
int j=0;
j<v_size;
j++)
1290 if(m > v_size || v_size < 1)
break;
1292 for(
int j=0;
j<v_size;
j++)
1301 if(m > v_size || v_size < 1)
break;
1303 for(
int j=0;
j<v_size;
j++){
1313 out.
data[
i] /= (a>0.) ? a : 1.;
1318 if(m > v_size || v_size < 1)
break;
1342 list<linedata>::iterator it =
lineList.begin();
1350 size_t max_size = 0;
1351 for(
size_t i=0;
i<l_size;
i++) {
1354 if(v_size > max_size) max_size = v_size;
1357 int m = (5*max_size+4);
1358 size_t n = m*(l_size+1);
1359 if(n < 4)
return false;
1364 out->
data[0] = max_size;
1365 out->
data[1] = l_size;
1376 double gpsTime = it->T_current;
1378 int gps =
int(gpsTime)/1000;
1379 out->
data[4] = float(gps);
1380 out->
data[5] = gpsTime - 1000.*double(gps);
1383 for(
size_t i=1;
i<=l_size;
i++) {
1392 for(
unsigned int j=0;
j<max_size;
j++){
1401 out->
data[
i*m+4+
j*5] = 0.;
1402 out->
data[
i*m+5+
j*5] = 0.;
1403 out->
data[
i*m+6+
j*5] = 0.;
1404 out->
data[
i*m+7+
j*5] = 0.;
1405 out->
data[
i*m+8+
j*5] = 0.;
1428 double gpsTime = 0.;
1429 unsigned int count = 0;
1432 if(in.
size() < 6)
return false;
1434 while(count < in.
size()){
1436 int m =
int(*(p+2) + 0.5);
1437 if(m <= 1)
return false;
1439 int l_size =
int(*(p+1) + 0.5);
1440 if(l_size < 1)
return false;
1442 int v_size =
int(*p + 0.5);
1443 count +=
int(*(p+3)+0.5);
1445 gpsTime = double(
int(*(p+4)+0.5))*1000. + *(p+5);
1452 v.
line.resize(v_size);
1453 v.
noise.resize(v_size);
1462 for(
int i=1;
i<=l_size;
i++) {
1464 if(*p != 0.) step = *p - last;
1472 for(
int j=0;
j<v_size;
j++){
std::list< linedata > lineList
virtual size_t size() const
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.
wavearray< double > a(hp.size())
void setFilter(int nF=1, int nL=0, int nS=1, int nD=-1, int nB=5, int nR=6, int nW=8)
void apply(WaveData &ts)
Operate on wavearray object.
std::vector< float > line
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)
wavearray< wavereal > WaveData
void reset()
Clear/release the internal History vector and reset the current time.
virtual void start(double s)
cout<< "Selected Pixels : "<< nPix<< endl;wc.cluster(1, 1);SSeries< double > ss
int Sample(Complex *H, double *t, int N, double totalMass, double Rate, wavearray< double > &hps, wavearray< double > &hxs)
linefilter * clone(void) const
Clone a linefilter.
std::vector< f_complex > amplitude
void sub(const wavearray< DataType_t > &, int=0, int=0, int=0)
bool LoadTrend(const char *)
int getLayer(wavearray< DataType_t > &w, double n)
param: n - layer number
linedata getLine(WaveData &)
wavearray< float > getTrend(int, char)
std::vector< float > filter
WaveData getPSD(const WaveData &, int=1)
void setFScan(double f=0., double sn=2., double fS=0.45, int nS=20)
double getOmega(const WaveData &, int=2)
virtual void DumpBinary(const char *, int=0)
double axb(double, double)
linefilter(void)
Build an empty linefilter.
bool DumpTrend(const char *, int=0)
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)
double fScan(const WaveData &)
linedata getHeteroLine(WaveData &)
double Stack(const wavearray< DataType_t > &, int)
std::complex< float > f_complex
~linefilter(void)
Destroy the linefilter object and release the function storage.
virtual void resize(unsigned int)
void Inverse(int n=-1)
param: n - number of steps (-1 means full reconstruction)
void putLayer(wavearray< DataType_t > &, double n)
param: n - layer number
double Interference(WaveData &, double)
std::vector< float > noise
void cpf(const wavearray< DataType_t > &, int=0, int=0, int=0)
unsigned int maxLine(int)