Commit 30c887d1 authored by Martin Reinecke's avatar Martin Reinecke

tweak trigonometric transforms

parent e0126108
......@@ -235,6 +235,8 @@ template<typename T> struct cmplx {
};
template<typename T> inline void PMINPLACE(T &a, T &b)
{ T t = a; a+=b; b=t-b; }
template<typename T> inline void MPINPLACE(T &a, T &b)
{ T t = a; a-=b; b=t+b; }
template<typename T> void PMC(T &a, T &b, const T &c, const T &d)
{ a = c+d; b = c-d; }
template<typename T> cmplx<T> conj(const cmplx<T> &a)
......@@ -2373,23 +2375,21 @@ template<typename T0> class T_dcst23
if (!cosine)
for (size_t k=1; k<N; k+=2)
c[k] = -c[k];
for (size_t i=2; i<N; i+=2)
for (size_t k=1; k<N-1; k+=2)
{
T xim1 = T0(0.5)*(c[i-1]+c[i]);
c[i] = T0(0.5)*(c[i]-c[i-1]);
c[i-1] = xim1;
T tmp = T0(0.5)*(c[k]+c[k+1]);
c[k+1] = T0(0.5)*(c[k+1]-c[k]);
c[k] = tmp;
}
fftplan.backward(c, fct);
for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
{
T tmp = twiddle[k-1]*c[kc]+twiddle[kc-1]*c[k];
c[kc] = twiddle[k-1]*c[k]-twiddle[kc-1]*c[kc];
c[k] = tmp;
T t1 = twiddle[k-1]*c[kc]+twiddle[kc-1]*c[k];
T t2 = twiddle[k-1]*c[k]-twiddle[kc-1]*c[kc];
c[k] = t1+t2; c[kc]=t1-t2;
}
if ((N&1)==0)
c[NS2] = twiddle[NS2-1]*(c[NS2]+c[NS2]);
for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
PMINPLACE(c[k], c[kc]);
c[0] *= 2;
if (!cosine)
for (size_t k=0, kc=N-1; k<kc; ++k, --kc)
......@@ -2402,30 +2402,22 @@ template<typename T0> class T_dcst23
{
constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
size_t N=length();
if (ortho) c[0]*=sqrt2;
size_t NS2 = (N+1)/2;
if (ortho) c[0]*=sqrt2;
if (!cosine)
for (size_t k=0, kc=N-1; k<NS2; ++k, --kc)
swap(c[k], c[kc]);
for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
PMINPLACE(c[k], c[kc]);
if ((N&1)==0)
c[NS2] = c[NS2]+c[NS2];
for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
{
T tmp = twiddle[k-1]*c[k]-twiddle[kc-1]*c[kc];
c[k] = twiddle[k-1]*c[kc]+twiddle[kc-1]*c[k];
c[kc] = tmp;
T t1=c[k]+c[kc], t2=c[k]-c[kc];
c[k] = twiddle[k-1]*t2+twiddle[kc-1]*t1;
c[kc]= twiddle[k-1]*t1-twiddle[kc-1]*t2;
}
if ((N&1)==0)
c[NS2] = twiddle[NS2-1]*c[NS2];
c[NS2] = 2*twiddle[NS2-1]*c[NS2];
fftplan.forward(c, fct);
for (size_t i=2; i<N; i+=2)
{
T xim1 = c[i-1]-c[i];
c[i] += c[i-1];
c[i-1] = xim1;
}
for (size_t k=1; k<N-1; k+=2)
MPINPLACE(c[k], c[k+1]);
if (!cosine)
for (size_t k=1; k<N; k+=2)
c[k] = -c[k];
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment