28 #include <xmmintrin.h> 34 #include "TFFTComplexReal.h" 35 #include "TFFTRealComplex.h" 57 static const
double Pi = 3.14159265358979312;
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;
68 void FFT(
double*
a,
double* b,
int n)
79 template<
class DataType_t>
85 template<
class DataType_t>
107 template<
class DataType_t>
141 printf(
"iNu too small, reset to 2\n");
145 printf(
"iNu too large, reset to 7\n");
156 double residual = 1 - tmp[0]*tmp[0];
157 double prec = pow(10., -
double(Precision));
160 residual -= 2*tmp[
N]*tmp[
N];
163 }
while(residual>prec || (N-1)%M2 || N/M2<3);
165 printf(
"Filter length = %d, norm = %.16f\n", N, 1.-residual);
175 template<
class DataType_t>
207 int K = m<0 ?
M : 2*M-
int(92*M/256.);
212 double residual = 1 - tmp[0]*tmp[0];
215 for(
int i=1;
i<=
M2;
i++) {
217 residual -= 2*tmp[
i]*tmp[
i];
227 filter *= 1./sqrt(1.-residual);
264 for(
int i=0;
i<6; ++
i) {
272 template<
class DataType_t>
291 template<
class DataType_t>
300 template<
class DataType_t>
322 template<
class DataType_t>
335 if(n%2)
for(
int i=-N+1;
i<
N; ++
i)w[
i] = 0;
342 if( (n-M)%2)
for(
int i=-N+1;
i<
N; ++
i)w[
i] = 0;
344 int s = (M%2) ? -1 : 1 ;
346 for(
int i=1;
i<
N; ++
i){
348 w[-
i] = w[
i] = s*wdmFilter[
i];
353 double ratio = m*
Pi/
M;
354 double sign = sqrt(2);
359 for(
int i=1;
i<
N; ++
i){
366 for(
int i=1;
i<
N; ++
i)
367 w[
i] = w[-
i] = sign*cos(
i*ratio)*wdmFilter[
i];
373 template<
class DataType_t>
390 else for(
int i=-N+1;
i<
N; ++
i)w[
i] = 0;
396 for(
int i=1;
i<
N; ++
i){
401 else for(
int i=-N+1;
i<
N; ++
i)w[
i] = 0;
404 double ratio = m*
Pi/
M;
405 double sign = sqrt(2);
410 for(
int i=1;
i<
N; ++
i)
411 w[
i] = w[-
i] = sign*cos(
i*ratio)*wdmFilter[
i];
415 for(
int i=1;
i<
N; ++
i){
425 template<
class DataType_t>
441 for(
int i=-nn;
i<=nn; ++
i) w[nn+
i] = w2[
i];
446 template<
class DataType_t>
455 for(
int i=0;
i<
N; ++
i){
462 template<
class DataType_t>
470 int layer =
int(
fabs(index)+0.001);
473 printf(
"WDM::getSlice(): illegal argument %d. Should be no more than %d\n",layer,M);
478 std::invalid_argument(
"WDM::getSlice(): data is not allocated");
483 size_t k = N/this->
nSTS>1 ? 1 : 0;
485 size_t i = index<0 ? layer+k*N/2 : layer;
489 if(i+(n-1)*s+1 > this->
nWWS){
490 std::invalid_argument(
"WaveDWT::getSlice(): invalide arguments");
497 template<
class DataType_t>
511 double A = (K-
M)*
Pi/2./K/M;
515 double gNorm = sqrt(2*M)/
Pi;
518 double* fourier =
new double[nFourier];
519 for(
int i=0;
i<nFourier; ++
i) fourier[
i] = Fourier[
i];
521 filter[0] = (A + fourier[0]*
B)*gNorm;
522 fNorm = filter[0]*filter[0];
523 fourier[0] /= sqrt(2);
525 for(
int i=1;
i<
n; ++
i){
528 double sumEven = 0, sumOdd = 0;
529 for(
int j=0;
j<nFourier; ++
j){
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];
538 if(
i/K < nFourier) intAB = fourier[
i/
K]*B/2*cos(
i*A);
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 );
542 fNorm += 2*filter[
i]*filter[
i];
550 template<
class DataType_t>
562 double A = (K-
M)*
Pi/2./K/M;
564 double gNorm = 2*M/
Pi;
568 double* fourier =
new double[nFourier];
569 for(
int i=0;
i<nFourier; ++
i) fourier[
i] = Fourier[
i];
571 filter[0] = (A + fourier[0]*
B)*gNorm;
572 fourier[0] /= sqrt(2);
574 for(
int i=1;
i<n*
L; ++
i){
575 double di =
i*(1./
L);
577 double sumEven = 0, sumOdd = 0;
578 for(
int j=0;
j<nFourier; ++
j){
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];
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 );
596 template<
class DataType_t>
604 double* res = tmp.
data;
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){
618 double di =
i*(1./
L);
621 for(
int j=0;
j<=nFourier;
j+=2){
622 if(
j*K*L==
i)
continue;
623 sum += di/(i2 -
j*
j*K2)*Fourier[
j];
625 sum *= 2*sin(di*B/2);
627 if(
i/(K*L) <= nFourier) {
628 if( (
i/(2*K*L)) & 1 ) sum -= Fourier[
i/K/
L]*B/2;
629 else sum += Fourier[
i/K/
L]*B/2;
631 res[
i] = gNorm*sqrt(2)*sum;
638 template<
class DataType_t>
651 for(
int i=0;
i<6; ++
i) {
659 for(
int i=M*L;
i>=-M*
L; --
i){
661 T0[
i][0] = tmp[abs(
i)];
662 for(
int j=1;
j<=nCoeffs; ++
j){
664 T0[
i][-
j] = tmp[M*L*
j -
i];
666 T0[
i].ZeroExtraElements();
671 for(
int i=M*L;
i>=-M*
L; --
i){
673 Tx[
i][0] = tmp[abs(
i)]/2.;
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.;
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;
685 Tx[
i].ZeroExtraElements();
692 for(
int i=0;
i<2*L*
M; ++
i){
700 template<
class DataType_t>
706 switch(
T0[0].SSESize()){
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;
747 template<
class DataType_t>
761 bool cancelEven = odd ^ Quad;
762 double dt = dT*(1./this->
LWDM);
763 if(m==0 || m==M)
if(cancelEven)
return 0;
764 DataType_t* TFMap = this->
pWWS;
765 if(Quad) TFMap += this->
nWWS/2;
768 DataType_t* tfc = TFMap + n*M1 +
m;
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]);
778 if(cancelEven)sumEven = 0;
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);
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]);
797 if(cancelEven)sumEven = 0;
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;
810 tfc = TFMap + n*M1 + m + 1;
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]);
817 if(cancelEven)sumEven = 0;
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));
823 if(m==0 || m==M-1) res1 *= sqrt(2);
824 if((m+n)&1) res -= res1;
832 template<
class DataType_t>
843 static const double sqrt2 = sqrt(2);
848 bool cancelEven = odd ^ Quad;
850 if(m==0 || m==M)
if(cancelEven)
return 0;
851 DataType_t* TFMap = this->
pWWS;
852 if(Quad) TFMap += this->
nWWS/2;
853 DataType_t* tfc = TFMap + n*M1 +
m;
861 const int last = sa.
Last();
866 float sumOdd = watasm_xmm0[1] + watasm_xmm0[3];
869 if(cancelEven)sumEven = 0;
875 if(index<0) index +=2*M*
L;
887 sumEven = watasm_xmm0[0] + watasm_xmm0[2];
888 sumOdd = watasm_xmm0[1] + watasm_xmm0[3];
891 if(cancelEven)sumEven = 0;
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];
900 if(m==1 || m == M) res1 *= sqrt2;
901 if((m+n)&1) res -= res1;
906 tfc = TFMap + n*M1 + m + 1;
910 sumEven = watasm_xmm0[0] + watasm_xmm0[2];
911 sumOdd = watasm_xmm0[1] + watasm_xmm0[3];
914 if(cancelEven)sumEven = 0;
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];
921 if(m==0 || m==M-1) res1 *= sqrt2;
922 if((m+n)&1) res -= res1;
928 template<
class DataType_t>
938 static const double sqrt2 = sqrt(2);
942 bool cancelEven = odd ^ Quad;
943 if(m==0 || m==M)
if(cancelEven)
return 0;
945 int QuadShift = Quad? 3:0;
953 float sumOdd = watasm_xmm0[1] + watasm_xmm0[3];
956 if(cancelEven)sumEven = 0;
961 if(index<0) index +=2*M*
L;
971 sumEven = watasm_xmm0[0] + watasm_xmm0[2];
972 sumOdd = watasm_xmm0[1] + watasm_xmm0[3];
975 if(cancelEven)sumEven = 0;
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];
984 if(m==1 || m == M) res1 *= sqrt2;
985 if((m+n)&1) res -= res1;
992 sumEven = watasm_xmm0[0] + watasm_xmm0[2];
993 sumOdd = watasm_xmm0[1] + watasm_xmm0[3];
996 if(cancelEven)sumEven = 0;
999 index = (2*m+1)*dT%(4*M*L);
1000 if(index<0)index += 4*M*
L;
1002 if(odd) res1 = (sumEven + sumOdd)*
sinTDx[index];
1003 else res1 = (sumEven - sumOdd)*
sinTDx[index];
1004 if(m==0 || m==M-1) res1 *= sqrt2;
1005 if((m+n)&1) res -= res1;
1011 template<
class DataType_t>
1021 static const double sqrt2 = sqrt(2);
1025 bool cancelEven = odd ^ Quad;
1026 if(m==0 || m==M)
if(cancelEven){
1027 for(
int dT = t1 ;
dT<=t2; ++
dT)*r++ = 0;
1032 int QuadShift = Quad? 3:0;
1036 for(
int dT = t1 ;
dT<=t2; ++
dT){
1042 float sumOdd = watasm_xmm0[1] + watasm_xmm0[3];
1045 if(cancelEven)sumEven = 0;
1050 if(index<0) index +=_2ml;
1060 sumEven = watasm_xmm0[0] + watasm_xmm0[2];
1061 sumOdd = watasm_xmm0[1] + watasm_xmm0[3];
1064 if(cancelEven)sumEven = 0;
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];
1073 if(m==1 || m == M) res1 *= sqrt2;
1074 if((m+n)&1) res -= res1;
1081 sumEven = watasm_xmm0[0] + watasm_xmm0[2];
1082 sumOdd = watasm_xmm0[1] + watasm_xmm0[3];
1085 if(cancelEven)sumEven = 0;
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];
1093 if(m==0 || m==M-1) res1 *= sqrt2;
1094 if((m+n)&1) res -= res1;
1102 template<
class DataType_t>
1113 printf(
"WDM:getTDamp() - time delays can not be produced with this TF data.");
1117 printf(
"WDM:getTDamp() - time delay filter is not set");
1129 if(n<0 || n>N/(M+1)) {
1130 cout<<
"WDM::getTDamp(): index outside TF map"<<endl;
exit(1);
1152 return float(a00*a00+a90*a90)/2;
1158 template<
class DataType_t>
1171 printf(
"WDM:getTDAmp() - time delays can not be produced with this TF data.");
1175 printf(
"WDM:getTDAmp() - time delay filter is not set");
1178 if(j>=(
int)this->
nWWS/2) j -= this->
nWWS/2;
1180 int mode = (c==
'a' || c==
'A') ? 2 : 1;
1183 float* p00 = amp.
data;
1184 float* p90 = amp.
data + (mode-1)*(2*K+1);
1187 for(
int k=-K;
k<=
K;
k++) {
1200 template<
class DataType_t>
1212 printf(
"WDM:getTDAmp() - time delay filter is not set");
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);
1235 int mode = (c==
'a' || c==
'A') ? 2 : 1;
1238 if(mode==1) aux.
resize(2*K+1);
1240 float* p00 = amp.
data;
1241 float* p90 = aux.
data;
1242 if(mode==2) p90 = amp.
data + (mode-1)*(2*K+1);
1244 if(min_wdmShift==max_wdmShift){
1256 for(
int i=min_wdmShift+2;
i<max_wdmShift;
i+=2){
1271 for(
int i=0;
i<=2*K+1; ++
i)p00[
i] = (p00[
i]*p00[
i] + p90[
i]*p90[
i])/2;
1276 template<
class DataType_t>
1283 printf(
"WDM:getTDAmp() - time delays can not be produced with this TF data.");
1287 printf(
"WDM:getTDAmp() - time delay filter is not set");
1290 DataType_t* TFMap = this->
pWWS;
1291 if(j>=(
int)this->
nWWS/2) j -= this->
nWWS/2;
1296 int BlockSize =
T0[0].SSESize()/4;
1298 float* data = r.
data;
1299 data += BlockSize/2;
1300 for(
int i=-L;
i<=
L; ++
i)data[
i] = TFMap[j +
i*M1];
1303 for(
int i=-L;
i<=
L; ++
i)data[
i] = TFMap[j +
i*M1];
1306 for(
int i=-L;
i<=
L; ++
i)data[
i] = TFMap[j +
i*M1];
1308 j+= this->
nWWS/2 - 1;
1310 for(
int i=-L;
i<=
L; ++
i)data[
i] = TFMap[j +
i*M1];
1313 for(
int i=-L;
i<=
L; ++
i)data[
i] = TFMap[j +
i*M1];
1316 for(
int i=-L;
i<=
L; ++
i)data[
i] = TFMap[j +
i*M1];
1320 template<
class DataType_t>
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)
1335 template<
class DataType_t>
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)
1351 template<
class DataType_t>
1359 template<
class DataType_t>
1367 template<
class DataType_t>
1372 printf(
"ERROR, WDM::forward called\n");
1375 template<
class DataType_t>
1380 printf(
"ERROR, WDM::inverse called\n");
1384 template<
class DataType_t>
1391 static const double sqrt2 = sqrt(2);
1397 int nWDM = this->
m_H;
1398 int nTS = this->
nWWS;
1407 m = this->
nWWS+2*nWDM;
1408 double*
ts = (
double*) _mm_malloc(m*
sizeof(
double), 16);
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]);
1413 double* pTS = ts + nWDM;
1417 double*
wdm = (
double*) _mm_malloc((nWDM)*
sizeof(double),16);
1418 double* INV = (
double*) _mm_malloc((nWDM)*
sizeof(double),16);
1419 for(n=0; n<nWDM; n++) {
1423 double*
WDM = INV+nWDM-1;
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));
1433 double* re =
new double[
M2];
1434 double* im =
new double[
M2];
1436 double* reX = (
double*) _mm_malloc(M2 *
sizeof(
double), 16);
1437 double* imX = (
double*) _mm_malloc(M2 *
sizeof(
double), 16);
1439 DataType_t* map00 =
pWDM;
1440 DataType_t* map90 = pWDM + N*
M1;
1443 int sign,kind; sign=0;
1444 TFFTRealComplex fft(1,&M2,
false);
1445 fft.Init(
"ES",sign, &kind);
1459 for(m=0; m<
M2; m++) {reX[
m] = 0; imX[
m] = 0.;}
1461 for(j=0; j<nWDM-1; j+=
M2) {
1464 __m128d* PR = (__m128d*) reX;
1465 for(m=0; m<
M2; m+=2) {
1466 __m128d ptj = _mm_loadu_pd(pTS+j+m);
1467 __m128d pTJ = _mm_loadu_pd(pTS-J+m);
1468 __m128d pwj = _mm_load_pd(wdm+j+m);
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));
1478 reX[0] += wdm[J]*pTS[J];
1483 fft.GetPointsComplex(reX, imX);
1494 reX[0] = imX[0] = reX[0]/sqrt2;
1495 reX[
M] = imX[
M] = reX[
M]/sqrt2;
1500 for(
int m=0; m<=
M; ++
m) {
1502 map00[
m] = sqrt2*imX[
m];
1503 map90[
m] = sqrt2*reX[
m];
1506 map00[
m] = sqrt2*reX[
m];
1507 map90[
m] = -sqrt2*imX[
m];
1513 for(
int m=0; m<=
M; m++)
1514 map00[m] = reX[m]*reX[m]+imX[m]*imX[m];
1538 template<
class DataType_t>
1551 int N = this->
nWWS/(M2+2);
1552 int nWDM = this->
m_H;
1553 int nTS = M*N + 2*nWDM;
1555 double sqrt2 = sqrt(2.);
1556 double* reX =
new double[
M2];
1557 double* imX =
new double[
M2];
1559 DataType_t* TFmap = this->
pWWS;
1561 if(M*N/this->
nSTS!=1) {
1562 printf(
"WDM<DataType_t>::w2t: Inverse is not defined for the up-sampled map\n");
1567 DataType_t*
tmp = ts.
data + nWDM;
1570 int sign,kind; sign=0;
1571 TFFTComplexReal ifft(1,&M2,
false);
1572 ifft.Init(
"ES",sign, &kind);
1577 for(j = 0; j<
M2; ++
j) reX[j] = imX[j] = 0;
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;
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];
1596 ifft.SetPointsComplex(reX,imX);
1598 ifft.GetPoints(reX);
1601 DataType_t*
x = tmp + n*
M;
1605 x[0] += wdm[0]*reX[
m];
1606 for(j = 1; j<nWDM; ++
j){
1609 x[
j] += wdm[
j]*reX[
m];
1610 if(mm==0)mm = 2*M - 1;
1612 x[-
j] += wdm[
j]*reX[mm];
1617 DataType_t* pTS = (DataType_t *)realloc(this->
pWWS,this->
nSTS*
sizeof(DataType_t));
1620 for(n=0; n<
int(this->
nSTS); n++) pTS[n] = tmp[n];
1630 template<
class DataType_t>
1639 int N = this->
nWWS/(M2+2);
1640 int nWDM = this->
m_H;
1641 int nTS = M*N + 2*nWDM;
1643 double sqrt2 = sqrt(2.);
1644 double* reX =
new double[
M2];
1645 double* imX =
new double[
M2];
1649 DataType_t* map90 = this->
pWWS + N*
M1;
1651 if(M*N/this->
nSTS!=1) {
1652 printf(
"WDM<DataType_t>::w2tQ: Inverse is not defined for the up-sampled map\n");
1657 DataType_t*
tmp = ts.
data + nWDM;
1660 int sign,kind; sign=0;
1661 TFFTComplexReal ifft(1,&M2,
false);
1662 ifft.Init(
"ES",sign, &kind);
1667 for(j = 0; j<
M2; ++
j) reX[j] = imX[j] = 0;
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;
1680 if((n+M) & 1)reX[
M] = map90[
M];
1681 if(n & 1)reX[0] = map90[0];
1683 ifft.SetPointsComplex(reX,imX);
1685 ifft.GetPoints(reX);
1688 DataType_t*
x = tmp + n*
M;
1692 x[0] += wdm[0]*reX[
m];
1693 for(j = 1; j<nWDM; ++
j){
1696 x[
j] += wdm[
j]*reX[
m];
1697 if(mm==0)mm = 2*M - 1;
1699 x[-
j] += wdm[
j]*reX[mm];
1705 DataType_t* pTS = (DataType_t *)realloc(this->
pWWS,this->
nSTS*
sizeof(DataType_t));
1708 for(n=0; n<
int(this->
nSTS); n++) pTS[n] = tmp[n];
1717 #define CLASS_INSTANTIATION(class_) template class WDM< class_ >; 1722 #undef CLASS_INSTANTIATION
int getBaseWave(int m, int n, SymmArray< double > &w)
static double SinCosSize[MAXBETA]
virtual WDM * Init() const
static double CosSize[MAXBETA]
wavearray< float > getTDvecSSE(int j, int K, char c, SSeries< double > *pss)
wavearray< double > a(hp.size())
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
WDM< double > wdm(nLAYERS, nLAYERS, 6, 10)
static double * SinCos[MAXBETA]
static double * Cos2[MAXBETA]
void setTDFilter(int nCoeffs, int L=1)
SymmObjArray< SymmArraySSE< float > > Tx
double getPixelAmplitude(int, int, int, bool=false)
void getTFvec(int j, wavearray< float > &r)
wavearray< float > sinTDx
#define CLASS_INSTANTIATION(class_)
virtual size_t size() const
wavearray< float > getTDvec(int j, int K, char c='p')
double getPixelAmplitudeSSEOld(int, int, int, bool=false)
void wavefft(double a[], double b[], int ntot, int n, int nspan, int isn)
int getBaseWaveQ(int m, int n, SymmArray< double > &w)
static double Cos2Size[MAXBETA]
bool GetSTFdata(int index, SymmArraySSE< float > *pS)
static double * Cos[MAXBETA]
void Resize(unsigned int sz)
void FFT(double *a, double *b, int n)
wavearray< double > getFilter(int n)
virtual WDM * Clone() const
return: Wavelet* - duplicate of *this, allocated on heap
printf("total live time: non-zero lags = %10.1f \, liveTot)
wavearray< double > getTDFilter2(int n, int L)
DataType_t ** TFMap90
pointer to 0-phase data, by default not initialized
void InvFFT(double *a, double *b, int n)
void(* SSE_TDF)()
pointer to 90-phase data, by default not initialized
float getPixelAmplitudeSSE(int m, int n, int dT, bool Quad)
unsigned long nWWS
pointer to wavelet work space
double TimeShiftTest(int dt)
float getTDamp(int n, int m, char c='p')
wavearray< double > ts(N)
wavearray< double > wdmFilter
double fabs(const Complex &x)
SymmArraySSE< float > td_halo[6]
std::slice getSlice(double n)
wavearray< double > getTDFilter1(int n, int L)
double TimeShiftTestSSE(int dt)
virtual void resize(unsigned int)
SymmObjArray< SymmArraySSE< float > > T0
void cpf(const wavearray< DataType_t > &, int=0, int=0, int=0)