diff --git a/pocketfft_hdronly.h b/pocketfft_hdronly.h index ac2ac7f85dd449bfb1aaa2ec408b522f1fd7b3d3..96cdf05f482b036e7372c8ebb749608ecd199326 100644 --- a/pocketfft_hdronly.h +++ b/pocketfft_hdronly.h @@ -2492,10 +2492,9 @@ template class T_dcst4 auto SGN_SET = [](T x, size_t i) {return (i%2) ? -x : x;}; arr y(N); size_t n2 = N/2; - size_t i; { - size_t m; - for (i=0, m=n2; m class T_dcst4 y[i] = c[m]; } rfft->forward(y.data(), fct); - for (i=0; i+i+1 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> y(N/2); for(size_t i=0; i void dct(const shape_t &shape, ndarr aout(data_out, shape, stride_out); if (type==1) general_dcst>(ain, aout, axes, fct, ortho, type, true, nthreads); - else if (type==2) - general_dcst>(ain, aout, axes, fct, ortho, type, true, nthreads); - else if (type==3) - general_dcst>(ain, aout, axes, fct, ortho, type, true, nthreads); else if (type==4) general_dcst>(ain, aout, axes, fct, ortho, type, true, nthreads); else - throw runtime_error("unsupported DCT type"); + general_dcst>(ain, aout, axes, fct, ortho, type, true, nthreads); } template void dst(const shape_t &shape, @@ -3362,14 +3362,10 @@ template void dst(const shape_t &shape, ndarr aout(data_out, shape, stride_out); if (type==1) general_dcst>(ain, aout, axes, fct, ortho, type, false, nthreads); - else if (type==2) - general_dcst>(ain, aout, axes, fct, ortho, type, false, nthreads); - else if (type==3) - general_dcst>(ain, aout, axes, fct, ortho, type, false, nthreads); else if (type==4) general_dcst>(ain, aout, axes, fct, ortho, type, false, nthreads); else - throw runtime_error("unsupported DST type"); + general_dcst>(ain, aout, axes, fct, ortho, type, false, nthreads); } template void r2c(const shape_t &shape_in,