Logo coherent WaveBurst  
Library Reference Guide
Logo
WDM.cc
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 Sergey Klimenko, Valentin Necula
3 #
4 # This program is free software: you can redistribute it and/or modify
5 # it under the terms of the GNU General Public License as published by
6 # the Free Software Foundation, either version 3 of the License, or
7 # (at your option) any later version.
8 #
9 # This program is distributed in the hope that it will be useful,
10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 # GNU General Public License for more details.
13 #
14 # You should have received a copy of the GNU General Public License
15 # along with this program. If not, see <https://www.gnu.org/licenses/>.
16 */
17 
18 
19 // Wavelet Analysis Tool
20 //--------------------------------------------------------------------
21 // V. Necula & S. Klimenko, University of Florida
22 // Implementation of Fast Wilson-Daubechies-Meyer transform
23 // Reference: http://iopscience.iop.org/1742-6596/363/1/012032
24 //--------------------------------------------------------------------
25 
26 #include <strstream>
27 #include <stdexcept>
28 #include <xmmintrin.h>
29 #include "WDM.hh"
30 #include "sseries.hh"
31 #include "wavefft.hh"
32 #include <iostream>
33 #include <stdio.h>
34 #include "TFFTComplexReal.h"
35 #include "TFFTRealComplex.h"
36 #include "TMath.h"
37 #include <complex>
38 
39 using namespace std;
40 
41 
42 extern "C" void sse_dp4();
43 extern "C" void sse_dp5();
44 extern "C" void sse_dp6();
45 extern "C" void sse_dp7();
46 extern "C" void sse_dp8();
47 extern "C" void sse_dp9();
48 extern "C" void sse_dp10();
49 extern "C" void sse_dp11();
50 
51 extern float* watasm_data;
52 extern float* watasm_filter;
53 extern float watasm_xmm0[4];
54 
55 ClassImp(WDM<double>)
56 
57 static const double Pi = 3.14159265358979312;
58 
59 template<class DataType_t> double* WDM<DataType_t>::Cos[MAXBETA];
60 template<class DataType_t> double* WDM<DataType_t>::Cos2[MAXBETA];
61 template<class DataType_t> double* WDM<DataType_t>::SinCos[MAXBETA];
62 template<class DataType_t> double WDM<DataType_t>::CosSize[MAXBETA];
63 template<class DataType_t> double WDM<DataType_t>::Cos2Size[MAXBETA];
64 template<class DataType_t> double WDM<DataType_t>::SinCosSize[MAXBETA];
65 template<class DataType_t> int WDM<DataType_t>::objCounter = 0;
66 
67 // exp^{-...} i.e. negative sign
68 void FFT(double* a, double* b, int n)
69 {
70  wavefft(a, b, n, n, n, -1);
71 }
72 
73 void InvFFT(double* a, double* b, int n)
74 {
75  wavefft(a,b, n,n,n, 1);
76 }
77 
78 // precalculated transform scaling functions
79 template<class DataType_t>
81 {
82  #include "FourierCoefficients.icc"
83 }
84 
85 template<class DataType_t>
87 WaveDWT<DataType_t>(1, 1, 0, B_CYCLE)
88 {
89 // default constructor
90 
91  if(++objCounter==1) initFourier();
92  this->m_WaveType = WDMT;
93  this->m_Heterodine = 1;
94  this->BetaOrder = 4;
95  this->precision = 10;
96  this->m_Level = 0;
97  this->m_Layer = 0;
98  this->KWDM = 0;
99  this->LWDM = 0;
100  wdmFilter.resize(0);
101  this->m_L = 0;
102  this->m_H = 0;
103  TFMap00 = TFMap90 = 0;
104  td_buffer=0;
105 }
106 
107 template<class DataType_t>
108 WDM<DataType_t>::WDM(int M, int K, int iNu, int Precision) :
109 WaveDWT<DataType_t>(1, 1, 0, B_CYCLE)
110 {
111 // constructor
112 // M + 1: number of bands
113 // K : K=n*M, where n is integer. K defines the width of the 'edge'
114 // of the basis function in Fourier domain (see the paper)
115 // larger n, longer the WDM filter, lower the spectral leakage between the bands.
116 // iNu : defines the sharpness of the 'edge' of the basis function in Fourier domain (see paper)
117 // Precison : defines filter length by truncation error quantified by
118 // P = -log10(1 - norm_of_filter) (see the paper)
119 // after forward transformation the structure of the WDM sliced array is the following:
120 // f phase 0 phase 90
121 // M * * ... * * * ... *
122 // ... ... ... ... ... ... ...
123 // 1 * * .... * * * ... *
124 // 0 * * ... * * * ... *
125 // t 1 2 ... n n+1 n+2... 2n
126 // where t/f is the time/frequency index
127 // the global TF index in the linear array is i = t*(M+1)+f
128 
129  if(++objCounter==1) initFourier();
130 
131  this->m_WaveType = WDMT;
132  this->m_Heterodine = 1;
133  this->BetaOrder = iNu;
134  this->precision = Precision;
135  this->m_Level = 0;
136  this->m_Layer = M;
137  this->KWDM = K;
138  this->LWDM = 0;
139 
140  if(iNu<2){
141  printf("iNu too small, reset to 2\n");
142  iNu = 2;
143  }
144  if(iNu>7){
145  printf("iNu too large, reset to 7\n");
146  iNu = 7;
147  }
148 
149  int nMax = 3e5;
150  int M2 = M*2;
151  int N = 1;
152 
153  wavearray<double> filter = this->getFilter(nMax);
154 
155  double* tmp = filter.data;
156  double residual = 1 - tmp[0]*tmp[0];
157  double prec = pow(10., -double(Precision));
158 
159  do {
160  residual -= 2*tmp[N]*tmp[N];
161  N++;
162  //printf("%d %e\n", N, residual);
163  } while(residual>prec || (N-1)%M2 || N/M2<3);
164 
165  printf("Filter length = %d, norm = %.16f\n", N, 1.-residual);
166 
167  wdmFilter.resize(N);
168  wdmFilter.cpf(filter,N);
169  this->m_L = 0;
170  this->m_H = N;
171  TFMap00 = TFMap90 = 0;
172  td_buffer=0;
173 }
174 
175 template<class DataType_t>
177 WaveDWT<DataType_t>(1, 1, 0, B_CYCLE)
178 {
179 // constructor
180 // M + 1: number of bands
181 // K : K=n*M, where n is integer. K defines the width of the 'edge'
182 // of the basis function in Fourier domain (see the paper)
183 // larger n, longer the WDM filter, lower the spectral leakage between the bands.
184 // after forward transformation the structure of the WDM sliced array is the following:
185 // f phase 0 phase 90
186 // M * * ... * * * ... *
187 // ... ... ... ... ... ... ...
188 // 1 * * .... * * * ... *
189 // 0 * * ... * * * ... *
190 // t 1 2 ... n n+1 n+2... 2n
191 // where t/f is the time/frequency index
192 // the global TF index in the linear array is i = t*(M+1)+f
193 
194  if(++objCounter==1) initFourier();
195  int M = abs(m);
196 
197  this->m_WaveType = WDMT;
198  this->m_Heterodine = 1;
199  this->BetaOrder = 2;
200  this->m_Level = 0;
201  this->m_Layer = M;
202  this->KWDM = M;
203  this->LWDM = 0;
204 
205  int nMax = 3e5;
206  int M2 = M*6;
207  int K = m<0 ? M : 2*M-int(92*M/256.); // trancate basis function
208 
209  wavearray<double> filter = this->getFilter(M2+6);
210 
211  double* tmp = filter.data;
212  double residual = 1 - tmp[0]*tmp[0];
213  //cout<<filter.size()<<" "<<M2<<" "<<K<<endl;
214 
215  for(int i=1; i<=M2; i++) {
216  if(i<K) {
217  residual -= 2*tmp[i]*tmp[i];
218  //printf("%d %e %e\n", i, residual, tmp[i]);
219  } else {
220  tmp[i] = 0.;
221  }
222  }
223 
224  //printf("Filter length = %d, norm = %.16f\n", M2+1, 1.-residual);
225 
226  wdmFilter.resize(M2+1);
227  filter *= 1./sqrt(1.-residual);
228  this->precision = -log10(residual);
229  wdmFilter.cpf(filter,M2+1);
230  this->m_L = 0;
231  this->m_H = M2+1;
232  TFMap00 = TFMap90 = 0;
233  td_buffer=0;
234 }
235 
236 template<class DataType_t> WDM<DataType_t>::
238 {
239 // copy constructor
240 // w : other WDM object
241 
242  ++objCounter;
243  this->m_Heterodine = 1;
244  this->m_WaveType = WDMT;
245  this->BetaOrder = w.BetaOrder;
246  this->precision = w.precision;
247  this->m_Level = w.m_Level;
248  this->m_Layer = w.m_Layer;
249  this->m_L = w.m_L;
250  this->m_H = w.m_H;
251  this->KWDM = w.KWDM;
252  this->LWDM = w.LWDM;
253  this->nSTS = w.nSTS; // WaveDWT parent class
254  this->wdmFilter = w.wdmFilter;
255  this->T0 = w.T0;
256  this->Tx = w.Tx;
257  TFMap00 = TFMap90 = 0;
258  sinTD = w.sinTD;
259  cosTD = w.cosTD;
260  sinTDx = w.sinTDx;
261  td_buffer = 0;
262  if(T0.Last()) {
263  initSSEPointers();
264  for(int i=0; i<6; ++i) {
265  td_halo[i].Resize(T0[0].Last());
267  }
268  }
269 
270 }
271 
272 template<class DataType_t>
274 {
275 // Destructor
276 
277  if(--objCounter==0) {
278  for(int i=2; i<MAXBETA; ++i){
279  if(Cos[i]) delete [] Cos[i];
280  if(Cos2[i]) delete [] Cos2[i];
281  if(SinCos[i]) delete [] SinCos[i];
282  }
283  }
284 
285  if(TFMap00) delete [] TFMap00;
286  if(TFMap90) delete [] TFMap90;
287  if(td_buffer)free(td_buffer);
288 }
289 
290 //clone
291 template<class DataType_t>
293 {
294 // Clone the object
295 
296  return new WDM<DataType_t>(*this);
297 }
298 
299 
300 template<class DataType_t>
302 {
303 // light-weight clone without TD filters
304 
305  WDM<DataType_t>* pwdm = new WDM<DataType_t>();
306  pwdm->m_Heterodine = this->m_Heterodine;
307  pwdm->m_WaveType = this->m_WaveType;
308  pwdm->BetaOrder = this->BetaOrder;
309  pwdm->precision = this->precision;
310  pwdm->m_Level = this->m_Level;
311  pwdm->m_Layer = this->m_Layer;
312  pwdm->m_L = this->m_L;
313  pwdm->m_H = this->m_H;
314  pwdm->KWDM = this->KWDM;
315  pwdm->LWDM = this->LWDM;
316  pwdm->nSTS = this->nSTS;
317  pwdm->wdmFilter = this->wdmFilter;
318  pwdm->TFMap00 = pwdm->TFMap90 = 0;
319  return pwdm;
320 }
321 
322 template<class DataType_t>
324 {
325 // computes basis function (sampled)
326 // m : frequency index
327 // n : time index
328 // w : where to store the values
329 // return value: indicates the time translation from origin needed by w (n=0 corresponds to t=0)
330 
331  int N = wdmFilter.size() ;
332  int M = this->m_Layer;
333  w.Resize(N-1);
334  if(m==0){
335  if(n%2)for(int i=-N+1; i<N; ++i)w[i] = 0;
336  else{
337  for(int i=1; i<N; ++i)w[-i] = w[i] = wdmFilter[i];
338  w[0] = wdmFilter[0];
339  }
340  }
341  else if(m==M){
342  if( (n-M)%2) for(int i=-N+1; i<N; ++i)w[i] = 0;
343  else {
344  int s = (M%2) ? -1 : 1 ; //sign
345  w[0] = s*wdmFilter[0];
346  for(int i=1; i<N; ++i){
347  s = -s;
348  w[-i] = w[i] = s*wdmFilter[i];
349  }
350  }
351  }
352  else {
353  double ratio = m*Pi/M;
354  double sign = sqrt(2);
355  //if((m*n)%2)sign = -sqrt(2);
356 
357  if((m+n)%2){
358  w[0] = 0;
359  for(int i=1; i<N; ++i){
360  w[i] = - sign*sin(i*ratio)*wdmFilter[i];
361  w[-i] = - w[i];
362  }
363  }
364  else{
365  w[0] = sign*wdmFilter[0];
366  for(int i=1; i<N; ++i)
367  w[i] = w[-i] = sign*cos(i*ratio)*wdmFilter[i];
368  }
369  }
370  return n*M;
371 }
372 
373 template<class DataType_t>
375 {
376 // computes Quadrature basis function (sampled)
377 // m : frequency index
378 // n : time index
379 // w : where to store the values
380 // return value: indicates the time translation from origin needed by w (n=0 corresponds to t=0)
381 
382  int N = wdmFilter.size() ;
383  int M = this->m_Layer;
384  w.Resize(N-1);
385  if(m==0){
386  if(n%2){
387  for(int i=1; i<N; ++i)w[-i] = w[i] = wdmFilter[i];
388  w[0] = wdmFilter[0];
389  }
390  else for(int i=-N+1; i<N; ++i)w[i] = 0;
391  }
392  else if(m==M){
393  if( (n-M)%2){
394  int s = 1;
395  w[0] = wdmFilter[0];
396  for(int i=1; i<N; ++i){
397  s = -s;
398  w[-i] = w[i] = s*wdmFilter[i];
399  }
400  }
401  else for(int i=-N+1; i<N; ++i)w[i] = 0;
402  }
403  else {
404  double ratio = m*Pi/M;
405  double sign = sqrt(2);
406  //if((m*n)%2)sign = -sqrt(2);
407 
408  if((m+n)%2){
409  w[0] = sign*wdmFilter[0];
410  for(int i=1; i<N; ++i)
411  w[i] = w[-i] = sign*cos(i*ratio)*wdmFilter[i];
412  }
413  else{
414  w[0] = 0;
415  for(int i=1; i<N; ++i){
416  w[i] = sign*sin(i*ratio)*wdmFilter[i];
417  w[-i] = - w[i];
418  }
419  }
420  }
421  return n*M;
422 }
423 
424 
425 template<class DataType_t>
427 {
428 // computes basis function (sampled)
429 // j : TF index
430 // w : where to store the values
431 // Quad: true to return Quadrature basis function
432 // return value: indicates the time translation from origin needed by w (j in [0,M] correspond to t=0)
433 
434  int M1 = this->m_Layer+1;
435  int m = j%M1;
436  int n = j/M1;
438  int shift = Quad? getBaseWaveQ(m, n, w2) : getBaseWave(m, n, w2);
439  int nn = w2.Last();
440  w.resize(2*nn+1);
441  for(int i=-nn; i<=nn; ++i) w[nn+i] = w2[i];
442  return shift-nn;
443 }
444 
445 
446 template<class DataType_t>
448 {
449 // not used
450 
451  int M1 = this->m_Layer+1;
452  int N = this->nWWS/2/M1;
453  TFMap00 = new DataType_t*[N];
454  TFMap90 = new DataType_t*[N];
455  for(int i=0; i<N; ++i){
456  TFMap00[i] = this->pWWS + M1*i;
457  TFMap90[i] = TFMap00[i] + M1*N;
458  }
459 }
460 
461 
462 template<class DataType_t>
464 {
465 // access function for a frequency layer described by index
466 // return slice corresponding to the selected frequency band
467 
468  int M = this->m_Level;
469  int N = this->nWWS;
470  int layer = int(fabs(index)+0.001);
471 
472  if(layer>M){
473  printf("WDM::getSlice(): illegal argument %d. Should be no more than %d\n",layer,M);
474  exit(0);
475  }
476 
477  if(!this->allocate()){
478  std::invalid_argument("WDM::getSlice(): data is not allocated");
479  return std::slice(0,1,1);
480  }
481 
482  size_t n = N/(M+1); // number of samples
483  size_t k = N/this->nSTS>1 ? 1 : 0; // power flag when = 0
484  size_t s = M+1; // slice step
485  size_t i = index<0 ? layer+k*N/2 : layer; // first sample
486 
487  n /= k+1; // adjust size
488 
489  if(i+(n-1)*s+1 > this->nWWS){
490  std::invalid_argument("WaveDWT::getSlice(): invalide arguments");
491  return std::slice(0,1,1);
492  }
493 
494  return std::slice(i,n,s);
495 }
496 
497 template<class DataType_t>
499 {
500 // internal WDM function
501 // computes the WDM transformation filter
502 // n - number of the filter coefficients defined by WDM
503 // parameters and precision.
504 
506  double* Fourier = Cos[BetaOrder];
507  int nFourier = CosSize[BetaOrder];
508  int M = this->m_Layer;
509  int K = this->KWDM;
510  double B = Pi/K;
511  double A = (K-M)*Pi/2./K/M;
512  double K2 = K;
513  K2 *= K2;
514  double* filter = tmp.data;
515  double gNorm = sqrt(2*M)/Pi;
516  double fNorm = 0.;
517 
518  double* fourier = new double[nFourier];
519  for(int i=0; i<nFourier; ++i) fourier[i] = Fourier[i];
520 
521  filter[0] = (A + fourier[0]*B)*gNorm; // do n = 0 first
522  fNorm = filter[0]*filter[0];
523  fourier[0] /= sqrt(2); // this line requires to copy Fourier to fourier
524 
525  for(int i=1; i<n; ++i){ // other coefficients
526  double di = i;
527  double i2 = di*di;
528  double sumEven = 0, sumOdd = 0;
529  for(int j=0; j<nFourier; ++j){ // a_j in the Fourier expansion
530  if(i%K);
531  else if(j==i/K)continue;
532  if(j&1)sumOdd += di/(i2 - j*j*K2)*fourier[j];
533  else sumEven += di/(i2 - j*j*K2)*fourier[j];
534  }
535 
536  double intAB = 0; //integral on [A, A+B]
537  if(i%K == 0)
538  if(i/K < nFourier) intAB = fourier[i/K]*B/2*cos(i*A);
539 
540  intAB += 2*(sumEven*sin(i*B/2)*cos(i*Pi/(2*M)) - sumOdd*sin(i*Pi/(2*M))*cos(i*B/2));
541  filter[i] = gNorm* ( sin(i*A)/i + sqrt(2)*intAB ); //sqrt2 b/c of the Four. expansion
542  fNorm += 2*filter[i]*filter[i];
543  }
544 
545  delete [] fourier;
546  return tmp;
547 
548 }
549 
550 template<class DataType_t>
551 wavearray<double> WDM<DataType_t>::getTDFilter1(int n, int L) // ex: L=8 -> tau/8 step
552 {
553 // computes 1st integral needed for TD filters (see the reference)
554 // n : defines the number of TD filter coefficients
555 // L : upsample factor the same as for setTDFilter
556 
557  double* Fourier = Cos2[BetaOrder];
558  int nFourier = Cos2Size[BetaOrder];
559  int M = this->m_Layer;
560  int K = this->KWDM;
561  double B = Pi/K;
562  double A = (K-M)*Pi/2./K/M;
563  double K2 = K*K;
564  double gNorm = 2*M/Pi;
565  wavearray<double> tmp(n*L);
566  double* filter = tmp.data;
567 
568  double* fourier = new double[nFourier];
569  for(int i=0; i<nFourier; ++i) fourier[i] = Fourier[i];
570 
571  filter[0] = (A + fourier[0]*B)*gNorm; // do n = 0 first
572  fourier[0] /= sqrt(2); // this line requires to copy Fourier to fourier
573 
574  for(int i=1; i<n*L; ++i){ // other coefficients
575  double di = i*(1./L);
576  double i2 = di*di;
577  double sumEven = 0, sumOdd = 0;
578  for(int j=0; j<nFourier; ++j){ // a_j in the Fourier expansion
579  if(j*K*L==i)continue;
580  if(j&1)sumOdd += di/(i2 - j*j*K2)*fourier[j];
581  else sumEven += di/(i2 - j*j*K2)*fourier[j];
582  }
583 
584  double intAB = 0; //integral on [A, A+B]
585  if(i%(K*L) == 0)if(i/(K*L) < nFourier)
586  intAB = fourier[i/(K*L)]*B/2*cos(di*A);
587  intAB += 2*(sumEven*sin(di*B/2)*cos(di*Pi/(2*M)) - sumOdd*sin(di*Pi/(2*M))*cos(di*B/2));
588  filter[i] = gNorm* ( sin(di*A)/di + sqrt(2)*intAB ); //sqrt2 b/c of the Four. expansion
589  }
590 
591  delete [] fourier;
592  return tmp;
593 }
594 
595 
596 template<class DataType_t>
598 {
599 // computes 2nd integral needed for TD filters (see paper)
600 // n : defines the number of TD filter coefficients
601 // L : upsample factor the same as for setTDFilter
602 
603  wavearray<double> tmp(n*L);
604  double* res = tmp.data;
605  int M = this->m_Layer;
606  int K = this->KWDM;
607  double B = Pi/K;
608  double K2 = K*K;
609 
610  double* Fourier = SinCos[BetaOrder];
611  int nFourier = SinCosSize[BetaOrder];
612  double gNorm = M/Pi;
613 
614  double aux = Fourier[0];
615  res[0] = aux*B*gNorm;
616  Fourier[0] /= sqrt(2);
617  for(int i=1; i<n*L; ++i){ // psi_i (l in the writeup)
618  double di = i*(1./L);
619  double i2 = di*di;
620  double sum = 0;
621  for(int j=0; j<=nFourier; j+=2){ //j indexing the Fourier coeff
622  if(j*K*L==i)continue;
623  sum += di/(i2 - j*j*K2)*Fourier[j];
624  }
625  sum *= 2*sin(di*B/2);
626  if(i%(2*K*L) == 0)
627  if(i/(K*L) <= nFourier) { // j*K*L = i case (even j)
628  if( (i/(2*K*L)) & 1 ) sum -= Fourier[i/K/L]*B/2;
629  else sum += Fourier[i/K/L]*B/2;
630  }
631  res[i] = gNorm*sqrt(2)*sum; //sqrt2 b/c of the Four. expansion
632  }
633  Fourier[0] = aux;
634  return tmp;
635 }
636 
637 
638 template<class DataType_t>
639 void WDM<DataType_t>::setTDFilter(int nCoeffs, int L)
640 {
641 // initialization of the time delay filters
642 // nCoeffs : define the number of the filter coefficients
643 // L : upsample factor, defines the fundamental time delay step
644 // dt = tau/L , where tau is the sampling interval of the original
645 // time series
646 
647  int M = this->m_Layer;
648  this->LWDM = L;
649  T0.Resize(M*L);
650  Tx.Resize(M*L);
651  for(int i=0; i<6; ++i) {
652  td_halo[i].Resize(nCoeffs);
654  }
655 
656  wavearray<double> filt1 = getTDFilter1(M*(nCoeffs+1)+1, L);
657 
658  double* tmp = filt1.data ; //cos2
659  for(int i=M*L; i>=-M*L; --i){
660  T0[i].Resize(nCoeffs);
661  T0[i][0] = tmp[abs(i)];
662  for(int j=1; j<=nCoeffs; ++j){
663  T0[i][j] = tmp[i+ M*L*j];
664  T0[i][-j] = tmp[M*L*j - i];
665  }
666  T0[i].ZeroExtraElements();
667  }
668 
669  wavearray<double> filt2 = getTDFilter2(M*(nCoeffs+1)+1 , L); //sincos
670  tmp = filt2.data;
671  for(int i=M*L; i>=-M*L; --i){
672  Tx[i].Resize(nCoeffs);
673  Tx[i][0] = tmp[abs(i)]/2.; // divide by 2 b/c the filter is actually for sin(2x) = 2*sin*cos
674  for(int j=1; j<=nCoeffs; ++j){
675  Tx[i][j] = tmp[i+ M*L*j]/2.;
676  Tx[i][-j] = tmp[M*L*j - i]/2.;
677  }
678 
679  // change some signs to bring i^l factor to C_l form!!
680  for(int j=1; j<=nCoeffs; ++j)switch(j%4){
681  case 1: Tx[i][-j] *= -1; break;
682  case 2: Tx[i][j] *= -1; Tx[i][-j] *= -1; break;
683  case 3: Tx[i][j] *= -1;
684  }
685  Tx[i].ZeroExtraElements();
686  }
687 
688  initSSEPointers();
689 
690  sinTD.resize(2*L*M);
691  cosTD.resize(2*L*M);
692  for(int i=0; i<2*L*M; ++i){
693  sinTD[i] = sin(i*Pi/(L*M));
694  cosTD[i] = cos(i*Pi/(L*M));
695  }
696  sinTDx.resize(4*L*M);
697  for(int i=0; i<4*L*M; ++i)sinTDx[i] = sin(i*Pi/(2*L*M));
698 }
699 
700 template<class DataType_t>
702 {
703 // required before using SSE instructions to speed up time delay filter calculations
704 
705  if(td_buffer)free(td_buffer);
706  switch(T0[0].SSESize()){
707  case 64: SSE_TDF = sse_dp4; posix_memalign((void**)&td_buffer, 16, 64+16); break;
708  case 80: SSE_TDF = sse_dp5; posix_memalign((void**)&td_buffer, 16, 80+16); break;
709  case 96: SSE_TDF = sse_dp6; posix_memalign((void**)&td_buffer, 16, 96+16); break;
710  case 112: SSE_TDF = sse_dp7; posix_memalign((void**)&td_buffer, 16, 112+16); break;
711  case 128: SSE_TDF = sse_dp8; posix_memalign((void**)&td_buffer, 16, 128+16); break;
712  case 144: SSE_TDF = sse_dp9; posix_memalign((void**)&td_buffer, 16, 144+16); break;
713  case 160: SSE_TDF = sse_dp10; posix_memalign((void**)&td_buffer, 16, 160+16); break;
714  case 176: SSE_TDF = sse_dp11; posix_memalign((void**)&td_buffer, 16, 176+16); break;
715  default: printf("initSSEPointer error, size not found. Contact V.Necula\n"); SSE_TDF=0;
716  }
717  td_data = td_buffer + T0[0].SSESize()/8 + 4;
718 }
719 
720 /*
721 template<class DataType_t>
722 double WDM<DataType_t>::PrintTDFilters(int m, int dt, int nC)
723 { SymmArray<double>& sa0 = T0[dt];
724 SymmArray<double>& sax = Tx[dt];
725 int M = this->m_Layer;
726 double sin0 = sin((m*Pi*dt)/M);
727 double cos0 = cos((m*Pi*dt)/M);
728 double sinP = sin((2*m+1)*dt*Pi/(2.*M));
729 double sinM = sin((2*m-1)*dt*Pi/(2.*M));
730 if(nC<=0 || nC>sa0.Last())nC = sa0.Last();
731 
732 double res = 0, t0, tp, tm;
733 for(int i=-nC; i<=nC; ++i){
734 if(i&1)t0 = sa0[i]*sin0;
735 else t0 = sa0[i]*cos0;
736 tp = sax[i]*sinP;
737 tm = sax[i]*sinM;
738 //printf("i=%3d %+le %+le %+le\n", i, t0, tp, tm);
739 res += t0*t0 + tm*tm + tp*tp;
740 }
741 return res;
742 }
743 */
744 
745 // Quadrature: equivalent to multiplying C_l by -i (plus shifting the 0,M layers)
746 
747 template<class DataType_t>
748 double WDM<DataType_t>::getPixelAmplitude(int m, int n, int dT, bool Quad)
749 {
750 // computes one time dalay for one pixel without SSE instructions
751 // m : frequency index (0...M)
752 // n : time index (in pixels)
753 // dT : time delay (in units of the sampling period or fractions of it,
754 // depending on the TD Filters setting
755 // Quad: true to compute the delayed amplitude for the quadrature [m,n] pixel
756 
757 
758  int M = this->m_Layer;
759  int M1 = M+1;
760  int odd = n&1;
761  bool cancelEven = odd ^ Quad;
762  double dt = dT*(1./this->LWDM);
763  if(m==0 || m==M)if(cancelEven) return 0; // DOESN'T WORK FOR ODD M!!
764  DataType_t* TFMap = this->pWWS;
765  if(Quad) TFMap += this->nWWS/2;
766 
767  //double* tfc = TFMap[m]+n;
768  DataType_t* tfc = TFMap + n*M1 + m;
769  double res, res1;
770 
771  //SymmArray<double>& sa = T0[dT];
772  SymmArraySSE<float>& sa = T0[dT];
773  double sumOdd = 0;
774  double sumEven = tfc[0]*sa[0];
775  for(int i=2; i<=sa.Last(); i+=2) sumEven += (tfc[i*M1]*sa[i] +tfc[-i*M1]*sa[-i]);
776  for(int i=1; i<=sa.Last(); i+=2) sumOdd += (tfc[i*M1]*sa[i] +tfc[-i*M1]*sa[-i]);
777  if(m==0 || m==M) {
778  if(cancelEven)sumEven = 0;
779  else sumOdd = 0;
780  }
781 
782  //res = 3.2*sumEven + sumOdd;
783  if(odd) res = sumEven*cos((m*Pi*dt)/M) - sumOdd*sin((m*Pi*dt)/M);
784  else res = sumEven*cos((m*Pi*dt)/M) + sumOdd*sin((m*Pi*dt)/M);
785 
786  //SymmArray<double>& sax = Tx[dT];
787  SymmArraySSE<float>& sax = Tx[dT];
788  if(m>0){ // m-1 case
789  // tfc = TFMap[m-1]+n;
790  tfc--;
791  sumOdd = 0;
792  sumEven = tfc[0]*sax[0];
793  for(int i=2; i<=sax.Last(); i+=2) sumEven += (tfc[i*M1]*sax[i] +tfc[-i*M1]*sax[-i]);
794  for(int i=1; i<=sax.Last(); i+=2) sumOdd += (tfc[i*M1]*sax[i] + tfc[-i*M1]*sax[-i]);
795 
796  if(m==1) {
797  if(cancelEven)sumEven = 0;
798  else sumOdd = 0;
799  }
800  //res1 = 2.1*sumEven + sumOdd;
801  if(odd) res1 = (sumEven - sumOdd)*sin((2*m-1)*dt*Pi/(2.*M));
802  else res1 = (sumEven + sumOdd)*sin((2*m-1)*dt*Pi/(2.*M));
803  if(m==1 || m == M) res1 *= sqrt(2);
804  if((m+n)&1) res -= res1;
805  else res += res1;
806  }
807 
808  if(m<M){ // m+1 case
809  //tfc = TFMap[m+1]+n;
810  tfc = TFMap + n*M1 + m + 1;
811  sumOdd = 0;
812  sumEven = tfc[0]*sax[0];
813  for(int i=2; i<=sax.Last(); i+=2) sumEven += (tfc[i*M1]*sax[i] + tfc[-i*M1]*sax[-i]);
814  for(int i=1; i<=sax.Last(); i+=2) sumOdd += (tfc[i*M1]*sax[i] + tfc[-i*M1]*sax[-i]);
815 
816  if(m==M-1) {
817  if(cancelEven)sumEven = 0;
818  else sumOdd = 0;
819  }
820  //res1 = 2.3*sumEven + sumOdd;
821  if(odd) res1 = (sumEven + sumOdd)*sin((2*m+1)*dt*Pi/(2.*M));
822  else res1 = (sumEven - sumOdd)*sin((2*m+1)*dt*Pi/(2.*M)); //two (-1)^l cancel
823  if(m==0 || m==M-1) res1 *= sqrt(2);
824  if((m+n)&1) res -= res1;
825  else res += res1;
826  }
827  return res;
828 }
829 
830 
831 // Quadrature: equivalent to multiplying C_l by -i (plus shifting the 0,M layers)
832 template<class DataType_t>
833 double WDM<DataType_t>::getPixelAmplitudeSSEOld(int m, int n, int dT, bool Quad)
834 {
835 // computes one time dalay for one pixel using SSE instructions
836 // m : frequency index (0...M)
837 // n : time index (in pixels)
838 // dT : time delay (in units of the sampling period or fractions of it,
839 // depending on the TD Filters setting
840 // Quad: true to compute the delayed amplitude for the quadrature [m,n] pixel
841 
842 
843  static const double sqrt2 = sqrt(2);
844  int M = this->m_Layer;
845  int M1 = M+1;
846  int L = this->LWDM;
847  int odd = n&1;
848  bool cancelEven = odd ^ Quad;
849 
850  if(m==0 || m==M)if(cancelEven) return 0; // DOESN'T WORK FOR ODD M!!
851  DataType_t* TFMap = this->pWWS;
852  if(Quad) TFMap += this->nWWS/2;
853  DataType_t* tfc = TFMap + n*M1 + m;
854 
855  double res, res1;
856 
857  SymmArraySSE<float>& sa = T0[dT];
858  watasm_data = td_buffer + 4;
859  watasm_filter = sa.SSEPointer();
860 
861  const int last = sa.Last();
862  for(int i=-last, j = -last*M1; i<=last; ++i){td_data[i] = tfc[j] ; j+=M1;}
863  SSE_TDF();
864 
865  float sumEven = watasm_xmm0[0] + watasm_xmm0[2];
866  float sumOdd = watasm_xmm0[1] + watasm_xmm0[3];
867 
868  if(m==0 || m==M) {
869  if(cancelEven)sumEven = 0;
870  else sumOdd = 0;
871  }
872 
873  //res = sumEven + sumOdd;
874  int index = (m*dT)%(2*M*L);
875  if(index<0) index +=2*M*L;
876 
877  if(odd) res = sumEven*cosTD[index] - sumOdd*sinTD[index];
878  else res = sumEven*cosTD[index] + sumOdd*sinTD[index];
879 
880  SymmArraySSE<float>& sax = Tx[dT];
881  watasm_filter = sax.SSEPointer();
882 
883  if(m>0){ // m-1 case
884  tfc--;
885  for(int i=-last, j = -last*M1; i<=last; ++i){td_data[i] = tfc[j] ; j+=M1;}
886  SSE_TDF();
887  sumEven = watasm_xmm0[0] + watasm_xmm0[2];
888  sumOdd = watasm_xmm0[1] + watasm_xmm0[3];
889 
890  if(m==1) {
891  if(cancelEven)sumEven = 0;
892  else sumOdd = 0;
893  }
894 
895  index = (2*m-1)*dT%(4*M*L);
896  if(index<0)index += 4*M*L;
897  if(odd) res1 = (sumEven - sumOdd)*sinTDx[index];
898  else res1 = (sumEven + sumOdd)*sinTDx[index];
899 
900  if(m==1 || m == M) res1 *= sqrt2;
901  if((m+n)&1) res -= res1;
902  else res += res1;
903  }
904 
905  if(m<M){ // m+1 case
906  tfc = TFMap + n*M1 + m + 1;
907  for(int i=-last, j = -last*M1; i<=last; ++i){td_data[i] = tfc[j] ; j+=M1;}
908  SSE_TDF();
909 
910  sumEven = watasm_xmm0[0] + watasm_xmm0[2];
911  sumOdd = watasm_xmm0[1] + watasm_xmm0[3];
912 
913  if(m==M-1) {
914  if(cancelEven)sumEven = 0;
915  else sumOdd = 0;
916  }
917  index = (2*m+1)*dT%(4*M*L);
918  if(index<0)index += 4*M*L;
919  if(odd) res1 = (sumEven + sumOdd)*sinTDx[index];
920  else res1 = (sumEven - sumOdd)*sinTDx[index]; //two (-1)^l cancel
921  if(m==0 || m==M-1) res1 *= sqrt2;
922  if((m+n)&1) res -= res1;
923  else res += res1;
924  }
925  return res;
926 }
927 
928 template<class DataType_t>
929 float WDM<DataType_t>::getPixelAmplitudeSSE(int m, int n, int dT, bool Quad)
930 {
931 // computes one time dalay for one pixel using SSE instructions; requires td_halo be initialized
932 // m : frequency index (0...M)
933 // n : time index (in pixels)
934 // dT : time delay (in units of the sampling period or fractions of it,
935 // depending on the TD Filters setting
936 // Quad: true to compute the delayed amplitude for the quadrature [m,n] pixel
937 
938  static const double sqrt2 = sqrt(2);
939  int M = this->m_Layer;
940  int L = this->LWDM;
941  int odd = n&1;
942  bool cancelEven = odd ^ Quad;
943  if(m==0 || m==M)if(cancelEven)return 0; // DOESN'T WORK FOR ODD M!!
944 
945  int QuadShift = Quad? 3:0;
946  double res, res1;
947 
948  SymmArraySSE<float>& sa = T0[dT];
949  watasm_data = td_halo[QuadShift + 1].SSEPointer();
950  watasm_filter = sa.SSEPointer();
951  SSE_TDF();
952  float sumEven = watasm_xmm0[0] + watasm_xmm0[2];
953  float sumOdd = watasm_xmm0[1] + watasm_xmm0[3];
954 
955  if(m==0 || m==M) {
956  if(cancelEven)sumEven = 0;
957  else sumOdd = 0;
958  }
959 
960  int index = (m*dT)%(2*M*L);
961  if(index<0) index +=2*M*L;
962  if(odd) res = sumEven*cosTD[index] - sumOdd*sinTD[index];
963  else res = sumEven*cosTD[index] + sumOdd*sinTD[index];
964 
965  SymmArraySSE<float>& sax = Tx[dT];
966  watasm_filter = sax.SSEPointer();
967 
968  if(m>0){ // m-1 case
969  watasm_data = td_halo[QuadShift + 0].SSEPointer();
970  SSE_TDF();
971  sumEven = watasm_xmm0[0] + watasm_xmm0[2];
972  sumOdd = watasm_xmm0[1] + watasm_xmm0[3];
973 
974  if(m==1) {
975  if(cancelEven)sumEven = 0;
976  else sumOdd = 0;
977  }
978 
979  index = (2*m-1)*dT%(4*M*L);
980  if(index<0)index += 4*M*L;
981  if(odd) res1 = (sumEven - sumOdd)*sinTDx[index];
982  else res1 = (sumEven + sumOdd)*sinTDx[index];
983 
984  if(m==1 || m == M) res1 *= sqrt2;
985  if((m+n)&1) res -= res1;
986  else res += res1;
987  }
988 
989  if(m<M){ // m+1 case
990  watasm_data = td_halo[QuadShift + 2].SSEPointer();
991  SSE_TDF();
992  sumEven = watasm_xmm0[0] + watasm_xmm0[2];
993  sumOdd = watasm_xmm0[1] + watasm_xmm0[3];
994 
995  if(m==M-1) {
996  if(cancelEven)sumEven = 0;
997  else sumOdd = 0;
998  }
999  index = (2*m+1)*dT%(4*M*L);
1000  if(index<0)index += 4*M*L;
1001 
1002  if(odd) res1 = (sumEven + sumOdd)*sinTDx[index];
1003  else res1 = (sumEven - sumOdd)*sinTDx[index]; //two (-1)^l cancel
1004  if(m==0 || m==M-1) res1 *= sqrt2;
1005  if((m+n)&1) res -= res1;
1006  else res += res1;
1007  }
1008  return res;
1009 }
1010 
1011 template<class DataType_t>
1012 void WDM<DataType_t>::getPixelAmplitudeSSE(int m, int n, int t1, int t2, float* r, bool Quad)
1013 {
1014 // computes time dalays for one pixel
1015 // m : frequency index (0...M)
1016 // n : time index (in pixels)
1017 // t1, t2 : time delays range [t1, t1+1, .. t2]
1018 // r : where to store the time delayed amplitude for pixel [m,n]
1019 // Quad : true to request time delays for the quadrature amplitude
1020 
1021  static const double sqrt2 = sqrt(2);
1022  int M = this->m_Layer;
1023  int L = this->LWDM;
1024  int odd = n&1;
1025  bool cancelEven = odd ^ Quad;
1026  if(m==0 || m==M)if(cancelEven){ // DOESN'T WORK FOR ODD M!!
1027  for(int dT = t1 ; dT<=t2; ++dT)*r++ = 0;
1028  return;
1029  }
1030 
1031  double res, res1;
1032  int QuadShift = Quad? 3:0;
1033  int _2ml = 2*M*L;
1034  int _4ml = 2*_2ml;
1035 
1036  for(int dT = t1 ; dT<=t2; ++dT){
1037 
1038  watasm_data = td_halo[QuadShift + 1].SSEPointer();
1039  watasm_filter = T0[dT].SSEPointer();
1040  SSE_TDF();
1041  float sumEven = watasm_xmm0[0] + watasm_xmm0[2];
1042  float sumOdd = watasm_xmm0[1] + watasm_xmm0[3];
1043 
1044  if(m==0 || m==M) {
1045  if(cancelEven)sumEven = 0;
1046  else sumOdd = 0;
1047  }
1048 
1049  int index = (m*dT)%_2ml;
1050  if(index<0) index +=_2ml;
1051  if(odd) res = sumEven*cosTD[index] - sumOdd*sinTD[index];
1052  else res = sumEven*cosTD[index] + sumOdd*sinTD[index];
1053 
1054  SymmArraySSE<float>& sax = Tx[dT];
1055  watasm_filter = sax.SSEPointer();
1056 
1057  if(m>0){ // m-1 case
1058  watasm_data = td_halo[QuadShift + 0].SSEPointer();
1059  SSE_TDF();
1060  sumEven = watasm_xmm0[0] + watasm_xmm0[2];
1061  sumOdd = watasm_xmm0[1] + watasm_xmm0[3];
1062 
1063  if(m==1) {
1064  if(cancelEven)sumEven = 0;
1065  else sumOdd = 0;
1066  }
1067 
1068  index = (2*m-1)*dT%_4ml;
1069  if(index<0)index += _4ml;
1070  if(odd) res1 = (sumEven - sumOdd)*sinTDx[index];
1071  else res1 = (sumEven + sumOdd)*sinTDx[index];
1072 
1073  if(m==1 || m == M) res1 *= sqrt2;
1074  if((m+n)&1) res -= res1;
1075  else res += res1;
1076  }
1077 
1078  if(m<M){ // m+1 case
1079  watasm_data = td_halo[QuadShift + 2].SSEPointer();
1080  SSE_TDF();
1081  sumEven = watasm_xmm0[0] + watasm_xmm0[2];
1082  sumOdd = watasm_xmm0[1] + watasm_xmm0[3];
1083 
1084  if(m==M-1) {
1085  if(cancelEven)sumEven = 0;
1086  else sumOdd = 0;
1087  }
1088 
1089  index = (2*m+1)*dT%_4ml;
1090  if(index<0)index += _4ml;
1091  if(odd) res1 = (sumEven + sumOdd)*sinTDx[index];
1092  else res1 = (sumEven - sumOdd)*sinTDx[index]; //two (-1)^l cancel
1093  if(m==0 || m==M-1) res1 *= sqrt2;
1094  if((m+n)&1) res -= res1;
1095  else res += res1;
1096  }
1097  *r++ = res;
1098  }
1099 }
1100 
1101 
1102 template<class DataType_t>
1103 float WDM<DataType_t>::getTDamp(int j, int k, char c)
1104 {
1105 // override getTDamp from WaveDWT
1106 // return time-delayed amplitude for delay index k:
1107 // j - global index in the TF map
1108 // k - delay index assuming time delay step 1/(LWDM*rate)
1109 // c - mode: 'a'/'A' - returns 00/90 amplitude, 'p'or'P' - returns power
1110 
1111  int N = (int)this->nWWS/2;
1112  if(this->nWWS/this->nSTS!=2) {
1113  printf("WDM:getTDamp() - time delays can not be produced with this TF data.");
1114  exit(0);
1115  }
1116  if(!this->Last()) {
1117  printf("WDM:getTDamp() - time delay filter is not set");
1118  exit(0);
1119  }
1120  if(j>=N) j -= N;
1121 
1122  int M = this->m_Layer;
1123  int J = M*LWDM; // number of delays in one pixel
1124  int n = j/(M+1); // time index
1125  int m = j%(M+1); // frequency index
1126  int wdmShift = k/J;
1127  k %= J;
1128 
1129  if(n<0 || n>N/(M+1)) {
1130  cout<<"WDM::getTDamp(): index outside TF map"<<endl; exit(1);
1131  }
1132 
1133 
1134  if(c=='a'){
1135  if(wdmShift%2){ // not working well for odd M when m=0/M !!
1136  if((n+m)%2) return -getPixelAmplitude(m, n-wdmShift, k, true);
1137  else return getPixelAmplitude(m, n-wdmShift, k, true);
1138  }
1139  return getPixelAmplitude(m, n-wdmShift, k, false);
1140  }
1141 
1142  if(c=='A'){
1143  if(wdmShift%2){ // not working well for odd M when m=0/M !!
1144  if((n+m)%2) return getPixelAmplitude(m, n-wdmShift, k, false);
1145  else return -getPixelAmplitude(m, n-wdmShift, k, false);
1146  }
1147  return getPixelAmplitude(m, n-wdmShift, k, true);
1148  }
1149 
1150  double a00 = getPixelAmplitude(m, n-wdmShift, k, false);
1151  double a90 = getPixelAmplitude(m, n-wdmShift, k, true);
1152  return float(a00*a00+a90*a90)/2;
1153 }
1154 
1155 
1156 
1157 // set array of time-delay amplitudes/energy
1158 template<class DataType_t>
1160 {
1161 // override getTDAmp from WaveDWT
1162 // return array of time-delayed amplitudes in the format:
1163 // [-K*dt , 0, K*dt, -K*dt, 0, K*dt] - total 2(2K+1) amplitudes
1164 // or array of time-delayed power in the format:
1165 // [-K*dt , 0, K*dt] - total (2K+1) amplitudes
1166 // j - index in TF map
1167 // K - range of unit delays dt
1168 // c - mode: 'a','A' - amplitude, 'p','P' - power
1169 
1170  if(this->nWWS/this->nSTS!=2) {
1171  printf("WDM:getTDAmp() - time delays can not be produced with this TF data.");
1172  exit(0);
1173  }
1174  if(!this->Last()) {
1175  printf("WDM:getTDAmp() - time delay filter is not set");
1176  exit(0);
1177  }
1178  if(j>=(int)this->nWWS/2) j -= this->nWWS/2;
1179 
1180  int mode = (c=='a' || c=='A') ? 2 : 1;
1181  wavearray<float> amp(mode*(2*K+1));
1182 
1183  float* p00 = amp.data;
1184  float* p90 = amp.data + (mode-1)*(2*K+1);
1185  int i = 0;
1186 
1187  for(int k=-K; k<=K; k++) {
1188  if(mode==2) {
1189  p00[i] = float(getTDamp(j, k, 'a'));
1190  p90[i] = float(getTDamp(j, k, 'A'));
1191  }
1192  else p00[i] = getTDamp(j, k, 'p');
1193  i++;
1194  }
1195  return amp;
1196 }
1197 
1198 
1199 // set array of time-delay amplitudes/energy
1200 template<class DataType_t>
1202 {
1203 // return array of time-delayed amplitudes in the format:
1204 // [-K*dt , 0, K*dt, -K*dt, 0, K*dt] - total 2(2K+1) amplitudes
1205 // or array of time-delayed power in the format:
1206 // [-K*dt , 0, K*dt] - total (2K+1) amplitudes
1207 // j - index in TF map
1208 // K - range of unit delays dt
1209 // c - mode: 'a'/'A' - amplitude, 'p','P' - power
1210 
1211  if(!this->Last()) {
1212  printf("WDM:getTDAmp() - time delay filter is not set");
1213  exit(0);
1214  }
1215 
1216  int M = this->m_Layer;
1217  int n = j/(M+1); // time index
1218  int m = j%(M+1); // frequency index
1219 
1220  int J = M*LWDM;
1221  int max_wdmShift = ((K + J)/(2*J))*2;
1222  int max_k = (K + J)%(2*J) - J;
1223  int min_wdmShift = (-K + J)/(2*J);
1224  int min_k = (-K + J)%(2*J);
1225  if(min_k<0){
1226  min_k += 2*J;
1227  --min_wdmShift;
1228  }
1229  min_wdmShift *=2;
1230  min_k -= J;
1231 
1232  //printf("minShift = %d maxShift = %d mink = %d maxk = %d\n",
1233  // min_wdmShift, max_wdmShift, min_k, max_k);
1234 
1235  int mode = (c=='a' || c=='A') ? 2 : 1;
1236  wavearray<float> amp(mode*(2*K+1));
1238  if(mode==1) aux.resize(2*K+1);
1239 
1240  float* p00 = amp.data;
1241  float* p90 = aux.data;
1242  if(mode==2) p90 = amp.data + (mode-1)*(2*K+1);
1243 
1244  if(min_wdmShift==max_wdmShift){
1245  pSS->GetSTFdata(j, td_halo);
1246  getPixelAmplitudeSSE(m, n, min_k, max_k, p00, false);
1247  getPixelAmplitudeSSE(m, n, min_k, max_k, p90, true);
1248  }
1249  else{
1250  pSS->GetSTFdata(j - min_wdmShift*(M+1), td_halo);
1251  getPixelAmplitudeSSE(m, n-min_wdmShift, min_k, J-1, p00, false);
1252  getPixelAmplitudeSSE(m, n-min_wdmShift, min_k, J-1, p90, true);
1253  p00 += J - min_k;
1254  p90 += J - min_k;
1255 
1256  for(int i=min_wdmShift+2; i<max_wdmShift; i+=2){
1257  pSS->GetSTFdata(j - i*(M+1), td_halo);
1258  getPixelAmplitudeSSE(m, n-i, -J, J-1, p00, false);
1259  getPixelAmplitudeSSE(m, n-i, -J, J-1, p90, true);
1260  p00 += 2*J;
1261  p90 += 2*J;
1262  }
1263  pSS->GetSTFdata(j - max_wdmShift*(M+1), td_halo);
1264  getPixelAmplitudeSSE(m, n-max_wdmShift, -J, max_k, p00, false);
1265  getPixelAmplitudeSSE(m, n-max_wdmShift, -J, max_k, p90, true);
1266  }
1267 
1268  if(mode==1){
1269  p00 = amp.data;
1270  p90 = aux.data;
1271  for(int i=0; i<=2*K+1; ++i)p00[i] = (p00[i]*p00[i] + p90[i]*p90[i])/2;
1272  }
1273  return amp;
1274 }
1275 
1276 template<class DataType_t>
1278 {
1279 // fill r with amplitudes needed to compute time-delays (for both quadratures)
1280 // j - global index in the TF map
1281 
1282  if(this->nWWS/this->nSTS!=2) {
1283  printf("WDM:getTDAmp() - time delays can not be produced with this TF data.");
1284  exit(0);
1285  }
1286  if(!this->Last()) {
1287  printf("WDM:getTDAmp() - time delay filter is not set");
1288  exit(0);
1289  }
1290  DataType_t* TFMap = this->pWWS;
1291  if(j>=(int)this->nWWS/2) j -= this->nWWS/2;
1292 
1293  int M = this->m_Layer;
1294  int M1 = M+1;
1295  int L = this->LWDM;
1296  int BlockSize = T0[0].SSESize()/4;
1297  r.resize(6*BlockSize);
1298  float* data = r.data;
1299  data += BlockSize/2;
1300  for(int i=-L; i<=L; ++i)data[i] = TFMap[j + i*M1];
1301  data += BlockSize;
1302  --j;
1303  for(int i=-L; i<=L; ++i)data[i] = TFMap[j + i*M1];
1304  data += BlockSize;
1305  j+=2;
1306  for(int i=-L; i<=L; ++i)data[i] = TFMap[j + i*M1];
1307 
1308  j+= this->nWWS/2 - 1; // 90 degree phase
1309  data += BlockSize;
1310  for(int i=-L; i<=L; ++i)data[i] = TFMap[j + i*M1];
1311  data += BlockSize;
1312  --j;
1313  for(int i=-L; i<=L; ++i)data[i] = TFMap[j + i*M1];
1314  data += BlockSize;
1315  j+=2;
1316  for(int i=-L; i<=L; ++i)data[i] = TFMap[j + i*M1];
1317 }
1318 
1319 
1320 template<class DataType_t>
1322 {
1323 // a rough CPU time test for time delays
1324  double res = 0;
1325  int M = this->m_Layer;
1326  int N = this->nWWS/2/(M+1);
1327  int nCoeff = T0[0].Last();
1328  for(int k=0; k<100; ++k)
1329  for(int i=nCoeff; i<N-nCoeff; ++i)
1330  for(int j=0; j<=M; ++j) res += getPixelAmplitude(j, i, dt);
1331  return res;
1332 }
1333 
1334 
1335 template<class DataType_t>
1337 {
1338 // a rough CPU time test for time delays using SSE instructions
1339 
1340  double res;
1341  int M = this->m_Layer;
1342  int N = this->nWWS/2/(M+1);
1343  int nCoeff = T0[0].Last();
1344  for(int k=0; k<100; ++k)
1345  for(int i=nCoeff; i<N-nCoeff; ++i)
1346  for(int j=0; j<=M; ++j) res += getPixelAmplitudeSSEOld(j, i, dt);
1347  return res;
1348 }
1349 
1350 
1351 template<class DataType_t>
1353 {
1354 // legacy function
1355 
1356  return this->m_Level;
1357 }
1358 
1359 template<class DataType_t>
1360 int WDM<DataType_t>::getOffset(int level, int layer)
1361 {
1362 // legacy function
1363 
1364  return layer;
1365 }
1366 
1367 template<class DataType_t>
1368 void WDM<DataType_t>::forward(int level, int layer)
1369 {
1370 // legacy function, should not be called
1371 
1372  printf("ERROR, WDM::forward called\n");
1373 }
1374 
1375 template<class DataType_t>
1376 void WDM<DataType_t>::inverse(int level,int layer)
1377 {
1378 // legacy function, should not be called
1379 
1380  printf("ERROR, WDM::inverse called\n");
1381 }
1382 
1383 
1384 template<class DataType_t>
1386 {
1387 // direct transform
1388 // MM = 0 requests power map of combined quadratures (not amplitudes for both)
1389 
1390  //pWWS contains the TD data, nWWS the size;
1391  static const double sqrt2 = sqrt(2);
1392 
1393  int n, m, j, J;
1394  int M = this->m_Layer; // max layer
1395  int M1 = M+1;
1396  int M2 = M*2;
1397  int nWDM = this->m_H;
1398  int nTS = this->nWWS;
1399  int KK = MM;
1400 
1401  if(MM<=0) MM = M;
1402 
1403  // adjust nWWS
1404  this->nWWS += this->nWWS%MM ? MM-this->nWWS%MM : 0;
1405 
1406  // initialize time series with boundary conditions (mirror)
1407  m = this->nWWS+2*nWDM;
1408  double* ts = (double*) _mm_malloc(m*sizeof(double), 16); // align to 16-byte for SSE
1409  for(n=0; n<=nWDM; n++) ts[nWDM-n] = (double)(this->pWWS[n]);
1410  for(n=0; n < nTS; n++) ts[nWDM+n] = (double)(this->pWWS[n]);
1411  for(n=0; n<int(m - nWDM-nTS); n++) ts[n+nWDM+nTS] = (double)(this->pWWS[nTS-n-1]);
1412 
1413  double* pTS = ts + nWDM;
1414 
1415 
1416  // create aligned symmetric arrays
1417  double* wdm = (double*) _mm_malloc((nWDM)*sizeof(double),16); // align to 16-byte for SSE
1418  double* INV = (double*) _mm_malloc((nWDM)*sizeof(double),16); // align to 16-byte for SSE
1419  for(n=0; n<nWDM; n++) {
1420  wdm[n] = wdmFilter.data[n];
1421  INV[n] = wdmFilter.data[nWDM-n-1];
1422  }
1423  double* WDM = INV+nWDM-1;
1424 
1425  // reallocate TF array
1426  int N = int(this->nWWS/MM);
1427  int L = KK<0 ? 2*N*M1 : N*M1;
1428  this->m_L = KK<0 ? this->m_H : 0;
1429  DataType_t* pWDM = (DataType_t *)realloc(this->pWWS,L*sizeof(DataType_t));
1430  this->release();
1431  this->allocate(size_t(L), pWDM);
1432 
1433  double* re = new double[M2];
1434  double* im = new double[M2];
1435 
1436  double* reX = (double*) _mm_malloc(M2 * sizeof(double), 16); // align to 16-byte for SSE
1437  double* imX = (double*) _mm_malloc(M2 * sizeof(double), 16); // align to 16-byte for SSE
1438 
1439  DataType_t* map00 = pWDM;
1440  DataType_t* map90 = pWDM + N*M1;
1441 
1442  int odd = 0;
1443  int sign,kind; sign=0; // dummy parameters
1444  TFFTRealComplex fft(1,&M2,false);
1445  fft.Init("ES",sign, &kind); // EX - optimized, ES - least optimized
1446 
1447  for(n=0; n<N; n++){
1448  /*
1449  for(m=0; m<M2; m++) {re[m] = 0; im[m] = 0.;}
1450  for(j=0; j<nWDM-1; j+=M2) {
1451  J = M2 + j;
1452  for(m=0; m<M2; m++)
1453  re[m] += *(pTS+j+m)*wdm[j+m] + *(pTS-J+m)*wdm[J-m];
1454  }
1455  re[0] += wdm[J]*pTS[J];
1456  printf("%d %.16f %.16f /",n,re[0],re[1]);
1457  */
1458 
1459  for(m=0; m<M2; m++) {reX[m] = 0; imX[m] = 0.;}
1460 
1461  for(j=0; j<nWDM-1; j+=M2) {
1462  J = M2 + j;
1463 
1464  __m128d* PR = (__m128d*) reX;
1465  for(m=0; m<M2; m+=2) {
1466  __m128d ptj = _mm_loadu_pd(pTS+j+m); // non-aligned pTS array
1467  __m128d pTJ = _mm_loadu_pd(pTS-J+m);
1468  __m128d pwj = _mm_load_pd(wdm+j+m); // aligned wdm and WDM arrays
1469  __m128d pWJ = _mm_load_pd(WDM-J+m);
1470  __m128d pjj = _mm_mul_pd(ptj,pwj);
1471  __m128d pJJ = _mm_mul_pd(pTJ,pWJ);
1472  *PR = _mm_add_pd(*PR, _mm_add_pd(pjj,pJJ));
1473  PR++;
1474  }
1475 
1476  }
1477 
1478  reX[0] += wdm[J]*pTS[J];
1479 // printf("%.16f %.16f \n",reX[0],reX[1]);
1480 
1481  fft.SetPoints(reX);
1482  fft.Transform();
1483  fft.GetPointsComplex(reX, imX);
1484 /*
1485  wavefft(re,im, M2,M2,M2, 1); // InvFFT(reX, imX, M2);
1486 
1487  printf("%d %.16f %.16f /",n,re[0],im[0]);
1488  printf("%.16f %.16f \n",reX[0],imX[0]);
1489  printf("%d %.16f %.16f /",n,re[1],im[1]);
1490  printf("%.16f %.16f \n",reX[1],imX[1]);
1491  printf("%d %.16f %.16f /",n,re[M],im[M]);
1492  printf("%.16f %.16f \n",reX[M],imX[M]);
1493 */
1494  reX[0] = imX[0] = reX[0]/sqrt2;
1495  reX[M] = imX[M] = reX[M]/sqrt2;
1496 
1497 // with fftw replace im with -im
1498 
1499  if(KK<0) {
1500  for(int m=0; m<=M; ++m) {
1501  if((m+odd) & 1){
1502  map00[m] = sqrt2*imX[m];
1503  map90[m] = sqrt2*reX[m];
1504  }
1505  else{
1506  map00[m] = sqrt2*reX[m];
1507  map90[m] = -sqrt2*imX[m];
1508  }
1509  }
1510  }
1511 
1512  else // power map
1513  for(int m=0; m<=M; m++)
1514  map00[m] = reX[m]*reX[m]+imX[m]*imX[m];
1515 
1516 
1517  odd = 1 - odd;
1518  pTS += MM;
1519  map00 += M1;
1520  map90 += M1;
1521  }
1522 
1523  this->m_Level = this->m_Layer;
1524 
1525  _mm_free(reX);
1526  _mm_free(imX);
1527  _mm_free(wdm);
1528  _mm_free(INV);
1529  _mm_free(ts);
1530 
1531  delete [] re;
1532  delete [] im;
1533 
1534 }
1535 
1536 
1537 
1538 template<class DataType_t>
1539 void WDM<DataType_t>::w2t(int flag)
1540 {
1541 // inverse transform
1542 // flag = -2 indicates that the Quadrature coefficients will be used
1543 
1544  if(flag==-2){
1545  w2tQ(0);
1546  return;
1547  }
1548  int n, j;
1549  int M = this->m_Layer;
1550  int M2 = M*2;
1551  int N = this->nWWS/(M2+2);
1552  int nWDM = this->m_H;
1553  int nTS = M*N + 2*nWDM;
1554 
1555  double sqrt2 = sqrt(2.);
1556  double* reX = new double[M2];
1557  double* imX = new double[M2];
1558  double* wdm = wdmFilter.data;
1559  DataType_t* TFmap = this->pWWS;
1560 
1561  if(M*N/this->nSTS!=1) {
1562  printf("WDM<DataType_t>::w2t: Inverse is not defined for the up-sampled map\n");
1563  exit(0);
1564  }
1565 
1567  DataType_t* tmp = ts.data + nWDM;
1568 
1569  int odd = 0;
1570  int sign,kind; sign=0; // dummy parameters
1571  TFFTComplexReal ifft(1,&M2,false);
1572  ifft.Init("ES",sign, &kind); // EX - optimized, ES - least optimized
1573 
1574  // Warning : if wavefft is used than reX/imX must be multiplied by sqrt(2)
1575  // if fftw is used than reX/imX must be divided by sqrt(2)
1576  for(n=0; n<N; ++n){
1577  for(j = 0; j<M2; ++j) reX[j] = imX[j] = 0;
1578  for(j=1; j<M; ++j)
1579  if( (n+j) & 1 ) imX[j] = TFmap[j]/sqrt2;
1580  else if(j&1) reX[j] = - TFmap[j]/sqrt2;
1581  else reX[j] = TFmap[j]/sqrt2;
1582 
1583  /* works only for EVEN M now
1584  if(n & 1);
1585  else{
1586  reX[0] = TFmap[0];
1587  reX[M] = TFmap[M];
1588  }
1589  */
1590 
1591  if( (n & 1) == 0 )reX[0] = TFmap[0];
1592  if( ((n+M) & 1) == 0 )
1593  if(M&1) reX[M] = -TFmap[M];
1594  else reX[M] = TFmap[M];
1595 
1596  ifft.SetPointsComplex(reX,imX);
1597  ifft.Transform();
1598  ifft.GetPoints(reX);
1599 
1600 // InvFFT(reX, imX, M2);
1601  DataType_t* x = tmp + n*M;
1602  int m = M*(n&1);
1603  int mm = m;
1604 
1605  x[0] += wdm[0]*reX[m];
1606  for(j = 1; j<nWDM; ++j){
1607  ++m;
1608  if(m==2*M)m=0;
1609  x[j] += wdm[j]*reX[m];
1610  if(mm==0)mm = 2*M - 1;
1611  else --mm;
1612  x[-j] += wdm[j]*reX[mm];
1613  }
1614  TFmap += M+1;
1615  }
1616 
1617  DataType_t* pTS = (DataType_t *)realloc(this->pWWS,this->nSTS*sizeof(DataType_t));
1618  this->release();
1619  this->allocate(size_t(this->nSTS), pTS);
1620  for(n=0; n<int(this->nSTS); n++) pTS[n] = tmp[n];
1621  this->m_Level = 0; // used in getSlice() function
1622 
1623  delete [] reX;
1624  delete [] imX;
1625 }
1626 
1627 
1628 // extra -i in C_k...
1629 
1630 template<class DataType_t>
1632 {
1633 // inverse transform using the quadrature coefficients
1634 
1635  int n, j;
1636  int M = this->m_Layer;
1637  int M1 = M+1;
1638  int M2 = M*2;
1639  int N = this->nWWS/(M2+2);
1640  int nWDM = this->m_H;
1641  int nTS = M*N + 2*nWDM;
1642 
1643  double sqrt2 = sqrt(2.);
1644  double* reX = new double[M2];
1645  double* imX = new double[M2];
1646  double* wdm = wdmFilter.data;
1647  //Data ype_t* TFmap = this->pWWS;
1648 
1649  DataType_t* map90 = this->pWWS + N*M1;
1650 
1651  if(M*N/this->nSTS!=1) {
1652  printf("WDM<DataType_t>::w2tQ: Inverse is not defined for the up-sampled map\n");
1653  exit(0);
1654  }
1655 
1657  DataType_t* tmp = ts.data + nWDM;
1658 
1659  int odd = 0;
1660  int sign,kind; sign=0; // dummy parameters
1661  TFFTComplexReal ifft(1,&M2,false);
1662  ifft.Init("ES",sign, &kind); // EX - optimized, ES - least optimized
1663 
1664  // Warning : if wavefft is used than reX/imX must be multiplied by sqrt(2)
1665  // if fftw is used than reX/imX must be divided by sqrt(2)
1666  for(n=0; n<N; ++n){
1667  for(j = 0; j<M2; ++j) reX[j] = imX[j] = 0;
1668  for(j=1; j<M; ++j)
1669  if( (n+j) & 1 ) reX[j] = map90[j]/sqrt2;
1670  else if(j&1) imX[j] = map90[j]/sqrt2;
1671  else imX[j] = -map90[j]/sqrt2;
1672 
1673  /* works only for EVEN M now
1674  if(n & 1){
1675  reX[0] = map90[0];
1676  reX[M] = map90[M];
1677  }
1678  */
1679 
1680  if((n+M) & 1)reX[M] = map90[M];
1681  if(n & 1)reX[0] = map90[0];
1682 
1683  ifft.SetPointsComplex(reX,imX);
1684  ifft.Transform();
1685  ifft.GetPoints(reX);
1686 
1687 // InvFFT(reX, imX, M2);
1688  DataType_t* x = tmp + n*M;
1689  int m = M*(n&1);
1690  int mm = m;
1691 
1692  x[0] += wdm[0]*reX[m];
1693  for(j = 1; j<nWDM; ++j){
1694  ++m;
1695  if(m==2*M)m=0;
1696  x[j] += wdm[j]*reX[m];
1697  if(mm==0)mm = 2*M - 1;
1698  else --mm;
1699  x[-j] += wdm[j]*reX[mm];
1700  }
1701 
1702  map90 += M+1;
1703  }
1704 
1705  DataType_t* pTS = (DataType_t *)realloc(this->pWWS,this->nSTS*sizeof(DataType_t));
1706  this->release();
1707  this->allocate(size_t(this->nSTS), pTS);
1708  for(n=0; n<int(this->nSTS); n++) pTS[n] = tmp[n];
1709  this->m_Level = 0; // used in getSlice() function
1710 
1711  delete [] reX;
1712  delete [] imX;
1713 }
1714 
1715 
1716 
1717 #define CLASS_INSTANTIATION(class_) template class WDM< class_ >;
1718 
1719 CLASS_INSTANTIATION(float)
1720 CLASS_INSTANTIATION(double)
1721 
1722 #undef CLASS_INSTANTIATION
void sse_dp5()
DataType_t ** TFMap00
Definition: WDM.hh:174
int getBaseWave(int m, int n, SymmArray< double > &w)
Definition: WDM.cc:323
double aux
Definition: testWDM_5.C:14
void w2t(int flag)
Definition: WDM.cc:1539
wavearray< float > cosTD
Definition: WDM.hh:172
static double SinCosSize[MAXBETA]
Definition: WDM.hh:138
#define MAXBETA
Definition: WDM.hh:37
TH1 * t1
double M
Definition: DrawEBHH.C:13
virtual WDM * Init() const
Definition: WDM.cc:301
static double CosSize[MAXBETA]
Definition: WDM.hh:138
wavearray< float > getTDvecSSE(int j, int K, char c, SSeries< double > *pss)
Definition: WDM.cc:1201
unsigned int Last()
Definition: SymmObjArray.hh:40
float M2
wavearray< double > a(hp.size())
#define B
int n
Definition: cwb_net.C:28
void initFourier()
Definition: WDM.cc:80
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
void sse_dp9()
WDM< double > wdm(nLAYERS, nLAYERS, 6, 10)
Record * SSEPointer()
Definition: SymmArraySSE.hh:41
float * watasm_data
STL namespace.
static double * SinCos[MAXBETA]
Definition: WDM.hh:137
static double * Cos2[MAXBETA]
Definition: WDM.hh:137
int m
Definition: cwb_net.C:28
float watasm_xmm0[4]
DataType_t * pWWS
Definition: WaveDWT.hh:141
int j
Definition: cwb_net.C:28
int getMaxLevel()
Definition: WDM.cc:1352
i drho i
int BetaOrder
Definition: WDM.hh:164
void setTDFilter(int nCoeffs, int L=1)
Definition: WDM.cc:639
SymmObjArray< SymmArraySSE< float > > Tx
Definition: WDM.hh:171
void sse_dp10()
void w2tQ(int)
Definition: WDM.cc:1631
#define N
void Resize(int sz)
Definition: SymmArray.cc:56
double getPixelAmplitude(int, int, int, bool=false)
Definition: WDM.cc:748
void inverse(int, int)
Definition: WDM.cc:1376
void release()
Definition: WaveDWT.cc:223
void sse_dp6()
int getOffset(int, int)
Definition: WDM.cc:1360
static int objCounter
Definition: WDM.hh:139
WDM()
Definition: WDM.cc:86
void getTFvec(int j, wavearray< float > &r)
Definition: WDM.cc:1277
wavearray< float > sinTDx
Definition: WDM.hh:172
size_t mode
wavearray< double > w
Definition: Test1.C:27
#define CLASS_INSTANTIATION(class_)
Definition: WDM.cc:1717
virtual size_t size() const
Definition: wavearray.hh:145
wavearray< float > getTDvec(int j, int K, char c='p')
Definition: WDM.cc:1159
double getPixelAmplitudeSSEOld(int, int, int, bool=false)
Definition: WDM.cc:833
void wavefft(double a[], double b[], int ntot, int n, int nspan, int isn)
Definition: wavefft.cc:41
int getBaseWaveQ(int m, int n, SymmArray< double > &w)
Definition: WDM.cc:374
int m_L
Definition: Wavelet.hh:124
static double Cos2Size[MAXBETA]
Definition: WDM.hh:138
bool GetSTFdata(int index, SymmArraySSE< float > *pS)
Definition: sseries.cc:362
unsigned long nSTS
Definition: WaveDWT.hh:143
wavearray< float > sinTD
Definition: WDM.hh:172
static double * Cos[MAXBETA]
Definition: WDM.hh:137
void Resize(unsigned int sz)
Definition: SymmObjArray.cc:54
i() int(T_cor *100))
void sse_dp8()
void SetTFMap()
Definition: WDM.cc:447
double * tmp
Definition: testWDM_5.C:31
void FFT(double *a, double *b, int n)
Definition: WDM.cc:68
wavearray< double > getFilter(int n)
Definition: WDM.cc:498
virtual WDM * Clone() const
return: Wavelet* - duplicate of *this, allocated on heap
Definition: WDM.cc:292
printf("total live time: non-zero lags = %10.1f \, liveTot)
wavearray< double > getTDFilter2(int n, int L)
Definition: WDM.cc:597
void sse_dp11()
Definition: Wavelet.hh:49
int LWDM
Definition: WDM.hh:167
float * td_data
Definition: WDM.hh:178
int KWDM
Definition: WDM.hh:166
void forward(int, int)
Definition: WDM.cc:1368
int k
static double A
Definition: geodesics.cc:26
DataType_t ** TFMap90
pointer to 0-phase data, by default not initialized
Definition: WDM.hh:175
int m_Level
Definition: Wavelet.hh:115
void InvFFT(double *a, double *b, int n)
Definition: WDM.cc:73
virtual ~WDM()
Definition: WDM.cc:273
int Last()
Definition: SymmArray.hh:41
void t2w(int)
Definition: WDM.cc:1385
bool allocate()
Definition: WaveDWT.cc:216
bool m_Heterodine
Definition: Wavelet.hh:126
void(* SSE_TDF)()
pointer to 90-phase data, by default not initialized
Definition: WDM.hh:176
float getPixelAmplitudeSSE(int m, int n, int dT, bool Quad)
Definition: WDM.cc:929
double dt
void Resize(int nn)
Definition: SymmArraySSE.cc:53
regression r
Definition: Regression_H1.C:44
s s
Definition: cwb_net.C:155
char filter[1024]
unsigned long nWWS
pointer to wavelet work space
Definition: WaveDWT.hh:142
size_t Last(int n=0)
Definition: WDM.hh:159
double TimeShiftTest(int dt)
Definition: WDM.cc:1321
float getTDamp(int n, int m, char c='p')
Definition: WDM.cc:1103
wavearray< double > ts(N)
wavearray< double > wdmFilter
Definition: WDM.hh:168
CosSize[2]
wavearray< int > index
WDM< double > * pWDM
Definition: testWDM_5.C:28
double fabs(const Complex &x)
Definition: numpy.cc:55
SymmArraySSE< float > td_halo[6]
Definition: WDM.hh:179
float M1
float * td_buffer
Definition: WDM.hh:177
std::slice getSlice(double n)
Definition: WDM.cc:463
wavearray< double > getTDFilter1(int n, int L)
Definition: WDM.cc:551
double TimeShiftTestSSE(int dt)
Definition: WDM.cc:1336
DataType_t * data
Definition: wavearray.hh:319
void sse_dp7()
enum WAVETYPE m_WaveType
Definition: Wavelet.hh:106
double dT
Definition: testWDM_5.C:12
static const double Pi
Definition: WDM.cc:57
int m_Layer
Definition: Wavelet.hh:118
void sse_dp4()
double shift[NIFO_MAX]
virtual void resize(unsigned int)
Definition: wavearray.cc:463
int precision
Definition: WDM.hh:165
void initSSEPointers()
Definition: WDM.cc:701
void ZeroExtraElements()
Definition: SymmArraySSE.cc:74
float * watasm_filter
int m_H
Definition: Wavelet.hh:121
SymmObjArray< SymmArraySSE< float > > T0
Definition: WDM.hh:170
exit(0)
void cpf(const wavearray< DataType_t > &, int=0, int=0, int=0)
Definition: wavearray.cc:717