Commit c9dc81da authored by Martin Reinecke's avatar Martin Reinecke

more polishing

parent c2da68d4
......@@ -2492,10 +2492,9 @@ template<typename T0> class T_dcst4
auto SGN_SET = [](T x, size_t i) {return (i%2) ? -x : x;};
arr<T> y(N);
size_t n2 = N/2;
size_t i;
{
size_t m;
for (i=0, m=n2; m<N; ++i, m+=4)
size_t i=0, m=n2;
for (; m<N; ++i, m+=4)
y[i] = c[m];
for (; m<2*N; ++i, m+=4)
y[i] = -c[2*N-m-1];
......@@ -2508,7 +2507,9 @@ template<typename T0> class T_dcst4
y[i] = c[m];
}
rfft->forward(y.data(), fct);
for (i=0; i+i+1<n2; ++i)
{
size_t i=0;
for (; i+i+1<n2; ++i)
{
size_t k = i+i+1;
T c1=y[2*k-1], s1=y[2*k], c2=y[2*k+1], s2=y[2*k+2];
......@@ -2523,12 +2524,15 @@ template<typename T0> class T_dcst4
c[i] = sqrt2 * (SGN_SET(cx, (i+1)/2) + SGN_SET(sx, i/2));
c[N-(i+1)] = sqrt2 * (SGN_SET(cx, (i+2)/2) + SGN_SET(sx, (i+1)/2));
}
}
c[n2] = sqrt2 * SGN_SET(y[0], (n2+1)/2);
// FFTW-derived code ends here
}
else
{
// even length algorithm from
// https://www.appletonaudio.com/blog/2013/derivation-of-fast-dct-4-algorithm-based-on-dft/
arr<cmplx<T>> y(N/2);
for(size_t i=0; i<N/2; ++i)
{
......@@ -3341,14 +3345,10 @@ template<typename T> void dct(const shape_t &shape,
ndarr<T> aout(data_out, shape, stride_out);
if (type==1)
general_dcst<T_dct1<T>>(ain, aout, axes, fct, ortho, type, true, nthreads);
else if (type==2)
general_dcst<T_dcst23<T>>(ain, aout, axes, fct, ortho, type, true, nthreads);
else if (type==3)
general_dcst<T_dcst23<T>>(ain, aout, axes, fct, ortho, type, true, nthreads);
else if (type==4)
general_dcst<T_dcst4<T>>(ain, aout, axes, fct, ortho, type, true, nthreads);
else
throw runtime_error("unsupported DCT type");
general_dcst<T_dcst23<T>>(ain, aout, axes, fct, ortho, type, true, nthreads);
}
template<typename T> void dst(const shape_t &shape,
......@@ -3362,14 +3362,10 @@ template<typename T> void dst(const shape_t &shape,
ndarr<T> aout(data_out, shape, stride_out);
if (type==1)
general_dcst<T_dst1<T>>(ain, aout, axes, fct, ortho, type, false, nthreads);
else if (type==2)
general_dcst<T_dcst23<T>>(ain, aout, axes, fct, ortho, type, false, nthreads);
else if (type==3)
general_dcst<T_dcst23<T>>(ain, aout, axes, fct, ortho, type, false, nthreads);
else if (type==4)
general_dcst<T_dcst4<T>>(ain, aout, axes, fct, ortho, type, false, nthreads);
else
throw runtime_error("unsupported DST type");
general_dcst<T_dcst23<T>>(ain, aout, axes, fct, ortho, type, false, nthreads);
}
template<typename T> void r2c(const shape_t &shape_in,
......
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