Commit d556ae21 authored by Martin Reinecke's avatar Martin Reinecke

cleanup DCT/DST, part 1

parent 656b2373
...@@ -2342,40 +2342,29 @@ template<typename T0> class T_dct2 ...@@ -2342,40 +2342,29 @@ template<typename T0> class T_dct2
{ {
constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L); constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
size_t N=length(); size_t N=length();
if (N==1) size_t NS2 = (N+1)/2;
c[0]*=2*fct; for (size_t i=2; i<N; i+=2)
else if (N==2)
{ {
T x1 = 2*fct*(c[0]+c[1]); T xim1 = T0(0.5)*(c[i-1]+c[i]);
c[1] = sqrt2*fct*(c[0]-c[1]); c[i] = T0(0.5)*(c[i]-c[i-1]);
c[0] = x1; c[i-1] = xim1;
} }
else fftplan.backward(c, fct);
for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
{ {
size_t NS2 = (N+1)/2; T tmp = twiddle[k-1]*c[kc]+twiddle[kc-1]*c[k];
for (size_t i=2; i<N; i+=2) c[kc] = twiddle[k-1]*c[k]-twiddle[kc-1]*c[kc];
{ c[k] = tmp;
T xim1 = T0(0.5)*(c[i-1]+c[i]); }
c[i] = T0(0.5)*(c[i]-c[i-1]); if ((N&1)==0)
c[i-1] = xim1; c[NS2] = twiddle[NS2-1]*(c[NS2]+c[NS2]);
} for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
fftplan.backward(c, fct); {
for (size_t k=1, kc=N-1; k<NS2; ++k, --kc) T tmp = c[k]+c[kc];
{ c[kc] = c[k]-c[kc];
T tmp = twiddle[k-1]*c[kc]+twiddle[kc-1]*c[k]; c[k] = tmp;
c[kc] = twiddle[k-1]*c[k]-twiddle[kc-1]*c[kc];
c[k] = tmp;
}
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)
{
T tmp = c[k]+c[kc];
c[kc] = c[k]-c[kc];
c[k] = tmp;
}
c[0] *= 2;
} }
c[0] *= 2;
if (ortho) c[0]/=sqrt2; if (ortho) c[0]/=sqrt2;
} }
...@@ -2402,40 +2391,29 @@ template<typename T0> class T_dct3 ...@@ -2402,40 +2391,29 @@ template<typename T0> class T_dct3
constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L); constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
size_t N=length(); size_t N=length();
if (ortho) c[0]*=sqrt2; if (ortho) c[0]*=sqrt2;
if (N==1) size_t NS2 = (N+1)/2;
c[0]*=fct; for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
else if (N==2)
{ {
T TSQX = sqrt2*c[1]; T tmp = c[k]-c[kc];
c[1] = fct*(c[0]-TSQX); c[k] = c[k]+c[kc];
c[0] = fct*(c[0]+TSQX); c[kc] = tmp;
} }
else if ((N&1)==0)
c[NS2] = c[NS2]+c[NS2];
for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
{ {
size_t NS2 = (N+1)/2; T tmp = twiddle[k-1]*c[k]-twiddle[kc-1]*c[kc];
for (size_t k=1, kc=N-1; k<NS2; ++k, --kc) c[k] = twiddle[k-1]*c[kc]+twiddle[kc-1]*c[k];
{ c[kc] = tmp;
T tmp = c[k]-c[kc]; }
c[k] = c[k]+c[kc]; if ((N&1)==0)
c[kc] = tmp; c[NS2] = twiddle[NS2-1]*c[NS2];
} fftplan.forward(c, fct);
if ((N&1)==0) for (size_t i=2; i<N; i+=2)
c[NS2] = c[NS2]+c[NS2]; {
for (size_t k=1, kc=N-1; k<NS2; ++k, --kc) T xim1 = c[i-1]-c[i];
{ c[i] += c[i-1];
T tmp = twiddle[k-1]*c[k]-twiddle[kc-1]*c[kc]; c[i-1] = xim1;
c[k] = twiddle[k-1]*c[kc]+twiddle[kc-1]*c[k];
c[kc] = tmp;
}
if ((N&1)==0)
c[NS2] = 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;
}
} }
} }
...@@ -2577,16 +2555,11 @@ template<typename T0> class T_dst2 ...@@ -2577,16 +2555,11 @@ template<typename T0> class T_dst2
{ {
constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L); constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
size_t N=length(); size_t N=length();
if (N==1) for (size_t k=1; k<N; k+=2)
c[0]*=2*fct; c[k] = -c[k];
else dct.exec(c, fct, false);
{ for (size_t k=0, kc=N-1; k<kc; ++k, --kc)
for (size_t k=1; k<N; k+=2) swap(c[k], c[kc]);
c[k] = -c[k];
dct.exec(c, fct, false);
for (size_t k=0, kc=N-1; k<kc; ++k, --kc)
swap(c[k], c[kc]);
}
if (ortho) c[0]/=sqrt2; if (ortho) c[0]/=sqrt2;
} }
...@@ -2607,17 +2580,12 @@ template<typename T0> class T_dst3 ...@@ -2607,17 +2580,12 @@ template<typename T0> class T_dst3
constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L); constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
size_t N=length(); size_t N=length();
if (ortho) c[0]*=sqrt2; if (ortho) c[0]*=sqrt2;
if (N==1) size_t NS2 = N/2;
c[0]*=fct; for (size_t k=0, kc=N-1; k<NS2; ++k, --kc)
else swap(c[k], c[kc]);
{ dct.exec(c, fct, false);
size_t NS2 = N/2; for (size_t k=1; k<N; k+=2)
for (size_t k=0, kc=N-1; k<NS2; ++k, --kc) c[k] = -c[k];
swap(c[k], c[kc]);
dct.exec(c, fct, false);
for (size_t k=1; k<N; k+=2)
c[k] = -c[k];
}
} }
size_t length() const { return dct.length(); } size_t length() const { return dct.length(); }
...@@ -2635,7 +2603,6 @@ template<typename T0> class T_dst4 ...@@ -2635,7 +2603,6 @@ template<typename T0> class T_dst4
template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool /*ortho*/) template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool /*ortho*/)
{ {
size_t N=length(); size_t N=length();
//if (N==1) { c[0]*=fct; return; }
size_t NS2 = N/2; size_t NS2 = N/2;
for (size_t k=0, kc=N-1; k<NS2; ++k, --kc) for (size_t k=0, kc=N-1; k<NS2; ++k, --kc)
swap(c[k], c[kc]); swap(c[k], c[kc]);
......
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