Commit 4525275e authored by Martin Reinecke's avatar Martin Reinecke

accurate DCT/DST IV for even lengths

parent 7e77dfec
......@@ -196,6 +196,13 @@ template<typename T> struct cmplx {
{ r+=other.r; i+=other.i; return *this; }
template<typename T2>cmplx &operator*= (T2 other)
{ r*=other; i*=other; return *this; }
template<typename T2>cmplx &operator*= (const cmplx<T2> &other)
{
T tmp = r*other.r - i*other.i;
i = r*other.i + i*other.r;
r = tmp;
return *this;
}
cmplx operator+ (const cmplx &other) const
{ return cmplx(r+other.r, i+other.i); }
cmplx operator- (const cmplx &other) const
......@@ -2221,10 +2228,10 @@ template<typename T0> class pocketfft_c
packplan=unique_ptr<cfftp<T0>>(new cfftp<T0>(length));
}
template<typename T> POCKETFFT_NOINLINE void backward(cmplx<T> c[], T0 fct)
template<typename T> POCKETFFT_NOINLINE void backward(cmplx<T> c[], T0 fct) const
{ packplan ? packplan->backward(c,fct) : blueplan->backward(c,fct); }
template<typename T> POCKETFFT_NOINLINE void forward(cmplx<T> c[], T0 fct)
template<typename T> POCKETFFT_NOINLINE void forward(cmplx<T> c[], T0 fct) const
{ packplan ? packplan->forward(c,fct) : blueplan->forward(c,fct); }
size_t length() const { return len; }
......@@ -2419,31 +2426,65 @@ template<typename T0> class T_dct3
template<typename T0> class T_dct4
{
private:
size_t N;
T_dct2<T0> dct2;
arr<T0> C;
pocketfft_c<T0> fft;
arr<cmplx<T0>> C2;
public:
POCKETFFT_NOINLINE T_dct4(size_t length)
: dct2(length), C(length)
: N(length), dct2((N&1)?N:1), C((N&1)?N:0),
fft((N&1)?1:N/2), C2((N&1) ? 0:N/2)
{
constexpr T0 pi = T0(3.141592653589793238462643383279502884197L);
for (size_t i=0; i<length; ++i)
C[i] = cos(T0(0.5)*pi*(T0(i)+T0(0.5))/T0(length));
if (N&1)
{
for (size_t i=0; i<N; ++i)
C[i] = cos(T0(0.5)*pi*(T0(i)+T0(0.5))/T0(N));
}
else
{
for (size_t i=0; i<N/2; ++i)
{
T0 ang = -pi/T0(N)*(T0(i)+T0(0.125));
C2[i].Set(cos(ang), sin(ang));
}
}
}
template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct) const
{
// FIXME: due to the recurrence below, the algorithm is not very accurate
// Trying to find a better one...
size_t N = length();
for (size_t i=0; i<N; ++i)
c[i] *= C[i];
dct2.exec(c, fct);
for (size_t i=1; i<N; ++i)
c[i] = 2*c[i] - c[i-1];
if (N&1)
{
for (size_t i=0; i<N; ++i)
c[i] *= C[i];
dct2.exec(c, fct);
for (size_t i=1; i<N; ++i)
c[i] = 2*c[i] - c[i-1];
}
else
{
arr<cmplx<T>> y(N/2);
for(size_t i=0; i<N/2; ++i)
{
y[i].Set(c[2*i],c[N-1-2*i]);
y[i] *= C2[i];
}
fft.forward(y.data(), fct);
for(size_t i=0; i<N/2; ++i)
y[i] *= C2[i];
for(size_t i=0; i<N/2; ++i)
{
c[2*i] = 2*y[i].r;
c[2*i+1] = -2*y[N/2-1-i].i;
}
}
}
size_t length() const { return dct2.length(); }
size_t length() const { return N; }
};
template<typename T0> class T_dst1
......
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