Chiebot-Mirror
/
opencv
8
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Merge pull request #18158 from legrosbuffle:3.4-vectorize-dft-radix

Browse Source
pull/18724/head
Alexander Alekhin 4 years ago
parent 97eaddd93b da555a2c9b
commit 716450ceb5
1 changed files with 280 additions and 109 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 389
      modules/core/src/dxt.cpp

389
modules/core/src/dxt.cpp
Unescape View File

@ -122,6 +122,33 @@ static const double DFTTab[][2] =
{ 1.00000000000000000, 0.00000000292583616 }
};
namespace {
template <typename T>
struct Constants {
static const T sin_120;
static const T fft5_2;
static const T fft5_3;
static const T fft5_4;
static const T fft5_5;
};
template <typename T>
const T Constants<T>::sin_120 = (T)0.86602540378443864676372317075294;
template <typename T>
const T Constants<T>::fft5_2 = (T)0.559016994374947424102293417182819;
template <typename T>
const T Constants<T>::fft5_3 = (T)-0.951056516295153572116439333379382;
template <typename T>
const T Constants<T>::fft5_4 = (T)-1.538841768587626701285145288018455;
template <typename T>
const T Constants<T>::fft5_5 = (T)0.363271264002680442947733378740309;
} //namespace
#define BitRev(i,shift) \
((int)((((unsigned)bitrevTab[(i)&255] << 24)+ \
((unsigned)bitrevTab[((i)>> 8)&255] << 16)+ \
@ -372,6 +399,149 @@ DFTInit( int n0, int nf, const int* factors, int* itab, int elem_size, void* _wa
}
}
// Reference radix-2 implementation.
template<typename T> struct DFT_R2
{
void operator()(Complex<T>* dst, const int c_n, const int n, const int dw0, const Complex<T>* wave) const {
const int nx = n/2;
for(int i = 0 ; i < c_n; i += n)
{
Complex<T>* v = dst + i;
T r0 = v[0].re + v[nx].re;
T i0 = v[0].im + v[nx].im;
T r1 = v[0].re - v[nx].re;
T i1 = v[0].im - v[nx].im;
v[0].re = r0; v[0].im = i0;
v[nx].re = r1; v[nx].im = i1;
for( int j = 1, dw = dw0; j < nx; j++, dw += dw0 )
{
v = dst + i + j;
r1 = v[nx].re*wave[dw].re - v[nx].im*wave[dw].im;
i1 = v[nx].im*wave[dw].re + v[nx].re*wave[dw].im;
r0 = v[0].re; i0 = v[0].im;
v[0].re = r0 + r1; v[0].im = i0 + i1;
v[nx].re = r0 - r1; v[nx].im = i0 - i1;
}
}
}
};
// Reference radix-3 implementation.
template<typename T> struct DFT_R3
{
void operator()(Complex<T>* dst, const int c_n, const int n, const int dw0, const Complex<T>* wave) const {
const int nx = n / 3;
for(int i = 0; i < c_n; i += n )
{
{
Complex<T>* v = dst + i;
T r1 = v[nx].re + v[nx*2].re;
T i1 = v[nx].im + v[nx*2].im;
T r0 = v[0].re;
T i0 = v[0].im;
T r2 = Constants<T>::sin_120*(v[nx].im - v[nx*2].im);
T i2 = Constants<T>::sin_120*(v[nx*2].re - v[nx].re);
v[0].re = r0 + r1; v[0].im = i0 + i1;
r0 -= (T)0.5*r1; i0 -= (T)0.5*i1;
v[nx].re = r0 + r2; v[nx].im = i0 + i2;
v[nx*2].re = r0 - r2; v[nx*2].im = i0 - i2;
}
for(int j = 1, dw = dw0; j < nx; j++, dw += dw0 )
{
Complex<T>* v = dst + i + j;
T r0 = v[nx].re*wave[dw].re - v[nx].im*wave[dw].im;
T i0 = v[nx].re*wave[dw].im + v[nx].im*wave[dw].re;
T i2 = v[nx*2].re*wave[dw*2].re - v[nx*2].im*wave[dw*2].im;
T r2 = v[nx*2].re*wave[dw*2].im + v[nx*2].im*wave[dw*2].re;
T r1 = r0 + i2; T i1 = i0 + r2;
r2 = Constants<T>::sin_120*(i0 - r2); i2 = Constants<T>::sin_120*(i2 - r0);
r0 = v[0].re; i0 = v[0].im;
v[0].re = r0 + r1; v[0].im = i0 + i1;
r0 -= (T)0.5*r1; i0 -= (T)0.5*i1;
v[nx].re = r0 + r2; v[nx].im = i0 + i2;
v[nx*2].re = r0 - r2; v[nx*2].im = i0 - i2;
}
}
}
};
// Reference radix-5 implementation.
template<typename T> struct DFT_R5
{
void operator()(Complex<T>* dst, const int c_n, const int n, const int dw0, const Complex<T>* wave) const {
const int nx = n / 5;
for(int i = 0; i < c_n; i += n )
{
for(int j = 0, dw = 0; j < nx; j++, dw += dw0 )
{
Complex<T>* v0 = dst + i + j;
Complex<T>* v1 = v0 + nx*2;
Complex<T>* v2 = v1 + nx*2;
T r0, i0, r1, i1, r2, i2, r3, i3, r4, i4, r5, i5;
r3 = v0[nx].re*wave[dw].re - v0[nx].im*wave[dw].im;
i3 = v0[nx].re*wave[dw].im + v0[nx].im*wave[dw].re;
r2 = v2[0].re*wave[dw*4].re - v2[0].im*wave[dw*4].im;
i2 = v2[0].re*wave[dw*4].im + v2[0].im*wave[dw*4].re;
r1 = r3 + r2; i1 = i3 + i2;
r3 -= r2; i3 -= i2;
r4 = v1[nx].re*wave[dw*3].re - v1[nx].im*wave[dw*3].im;
i4 = v1[nx].re*wave[dw*3].im + v1[nx].im*wave[dw*3].re;
r0 = v1[0].re*wave[dw*2].re - v1[0].im*wave[dw*2].im;
i0 = v1[0].re*wave[dw*2].im + v1[0].im*wave[dw*2].re;
r2 = r4 + r0; i2 = i4 + i0;
r4 -= r0; i4 -= i0;
r0 = v0[0].re; i0 = v0[0].im;
r5 = r1 + r2; i5 = i1 + i2;
v0[0].re = r0 + r5; v0[0].im = i0 + i5;
r0 -= (T)0.25*r5; i0 -= (T)0.25*i5;
r1 = Constants<T>::fft5_2*(r1 - r2); i1 = Constants<T>::fft5_2*(i1 - i2);
r2 = -Constants<T>::fft5_3*(i3 + i4); i2 = Constants<T>::fft5_3*(r3 + r4);
i3 *= -Constants<T>::fft5_5; r3 *= Constants<T>::fft5_5;
i4 *= -Constants<T>::fft5_4; r4 *= Constants<T>::fft5_4;
r5 = r2 + i3; i5 = i2 + r3;
r2 -= i4; i2 -= r4;
r3 = r0 + r1; i3 = i0 + i1;
r0 -= r1; i0 -= i1;
v0[nx].re = r3 + r2; v0[nx].im = i3 + i2;
v2[0].re = r3 - r2; v2[0].im = i3 - i2;
v1[0].re = r0 + r5; v1[0].im = i0 + i5;
v1[nx].re = r0 - r5; v1[nx].im = i0 - i5;
}
}
}
};
template<typename T> struct DFT_VecR2
{
void operator()(Complex<T>* dst, const int c_n, const int n, const int dw0, const Complex<T>* wave) const {
return DFT_R2<T>()(dst, c_n, n, dw0, wave);
}
};
template<typename T> struct DFT_VecR3
{
void operator()(Complex<T>* dst, const int c_n, const int n, const int dw0, const Complex<T>* wave) const {
return DFT_R3<T>()(dst, c_n, n, dw0, wave);
}
};
template<typename T> struct DFT_VecR4
{
int operator()(Complex<T>*, int, int, int&, const Complex<T>*) const { return 1; }
@ -379,6 +549,98 @@ template<typename T> struct DFT_VecR4
#if CV_SSE3
// multiplies *a and *b:
// r_re + i*r_im = (a_re + i*a_im)*(b_re + i*b_im)
// r_re and r_im are placed respectively in bits 31:0 and 63:32 of the resulting
// vector register.
inline __m128 complexMul(const Complex<float>* const a, const Complex<float>* const b) {
const __m128 z = _mm_setzero_ps();
const __m128 neg_elem0 = _mm_set_ps(0.0f,0.0f,0.0f,-0.0f);
// v_a[31:0] is a->re and v_a[63:32] is a->im.
const __m128 v_a = _mm_loadl_pi(z, (const __m64*)a);
const __m128 v_b = _mm_loadl_pi(z, (const __m64*)b);
// x_1 = v[nx] * wave[dw].
const __m128 v_a_riri = _mm_shuffle_ps(v_a, v_a, _MM_SHUFFLE(0, 1, 0, 1));
const __m128 v_b_irri = _mm_shuffle_ps(v_b, v_b, _MM_SHUFFLE(1, 0, 0, 1));
const __m128 mul = _mm_mul_ps(v_a_riri, v_b_irri);
const __m128 xored = _mm_xor_ps(mul, neg_elem0);
return _mm_hadd_ps(xored, z);
}
// optimized radix-2 transform
template<> struct DFT_VecR2<float> {
void operator()(Complex<float>* dst, const int c_n, const int n, const int dw0, const Complex<float>* wave) const {
const __m128 z = _mm_setzero_ps();
const int nx = n/2;
for(int i = 0 ; i < c_n; i += n)
{
{
Complex<float>* v = dst + i;
float r0 = v[0].re + v[nx].re;
float i0 = v[0].im + v[nx].im;
float r1 = v[0].re - v[nx].re;
float i1 = v[0].im - v[nx].im;
v[0].re = r0; v[0].im = i0;
v[nx].re = r1; v[nx].im = i1;
}
for( int j = 1, dw = dw0; j < nx; j++, dw += dw0 )
{
Complex<float>* v = dst + i + j;
const __m128 x_1 = complexMul(&v[nx], &wave[dw]);
const __m128 v_0 = _mm_loadl_pi(z, (const __m64*)&v[0]);
_mm_storel_pi((__m64*)&v[0], _mm_add_ps(v_0, x_1));
_mm_storel_pi((__m64*)&v[nx], _mm_sub_ps(v_0, x_1));
}
}
}
};
// Optimized radix-3 implementation.
template<> struct DFT_VecR3<float> {
void operator()(Complex<float>* dst, const int c_n, const int n, const int dw0, const Complex<float>* wave) const {
const int nx = n / 3;
const __m128 z = _mm_setzero_ps();
const __m128 neg_elem1 = _mm_set_ps(0.0f,0.0f,-0.0f,0.0f);
const __m128 sin_120 = _mm_set1_ps(Constants<float>::sin_120);
const __m128 one_half = _mm_set1_ps(0.5f);
for(int i = 0; i < c_n; i += n )
{
{
Complex<float>* v = dst + i;
float r1 = v[nx].re + v[nx*2].re;
float i1 = v[nx].im + v[nx*2].im;
float r0 = v[0].re;
float i0 = v[0].im;
float r2 = Constants<float>::sin_120*(v[nx].im - v[nx*2].im);
float i2 = Constants<float>::sin_120*(v[nx*2].re - v[nx].re);
v[0].re = r0 + r1; v[0].im = i0 + i1;
r0 -= (float)0.5*r1; i0 -= (float)0.5*i1;
v[nx].re = r0 + r2; v[nx].im = i0 + i2;
v[nx*2].re = r0 - r2; v[nx*2].im = i0 - i2;
}
for(int j = 1, dw = dw0; j < nx; j++, dw += dw0 )
{
Complex<float>* v = dst + i + j;
const __m128 x_0 = complexMul(&v[nx], &wave[dw]);
const __m128 x_2 = complexMul(&v[nx*2], &wave[dw*2]);
const __m128 x_1 = _mm_add_ps(x_0, x_2);
const __m128 v_0 = _mm_loadl_pi(z, (const __m64*)&v[0]);
_mm_storel_pi((__m64*)&v[0], _mm_add_ps(v_0, x_1));
const __m128 x_3 = _mm_mul_ps(sin_120, _mm_xor_ps(_mm_sub_ps(x_2, x_0), neg_elem1));
const __m128 x_3s = _mm_shuffle_ps(x_3, x_3, _MM_SHUFFLE(0, 1, 0, 1));
const __m128 x_4 = _mm_sub_ps(v_0, _mm_mul_ps(one_half, x_1));
_mm_storel_pi((__m64*)&v[nx], _mm_add_ps(x_4, x_3s));
_mm_storel_pi((__m64*)&v[nx*2], _mm_sub_ps(x_4, x_3s));
}
}
}
};
// optimized radix-4 transform
template<> struct DFT_VecR4<float>
{
@ -573,12 +835,6 @@ struct OcvDftOptions {
template<typename T> static void
DFT(const OcvDftOptions & c, const Complex<T>* src, Complex<T>* dst)
{
static const T sin_120 = (T)0.86602540378443864676372317075294;
static const T fft5_2 = (T)0.559016994374947424102293417182819;
static const T fft5_3 = (T)-0.951056516295153572116439333379382;
static const T fft5_4 = (T)-1.538841768587626701285145288018455;
static const T fft5_5 = (T)0.363271264002680442947733378740309;
const Complex<T>* wave = (Complex<T>*)c.wave;
const int * itab = c.itab;
@ -775,30 +1031,18 @@ DFT(const OcvDftOptions & c, const Complex<T>* src, Complex<T>* dst)
for( ; n < c.factors[0]; )
{
// do the remaining radix-2 transform
nx = n;
n *= 2;
dw0 /= 2;
for( i = 0; i < c.n; i += n )
if(c.haveSSE3)
{
Complex<T>* v = dst + i;
T r0 = v[0].re + v[nx].re;
T i0 = v[0].im + v[nx].im;
T r1 = v[0].re - v[nx].re;
T i1 = v[0].im - v[nx].im;
v[0].re = r0; v[0].im = i0;
v[nx].re = r1; v[nx].im = i1;
for( j = 1, dw = dw0; j < nx; j++, dw += dw0 )
{
v = dst + i + j;
r1 = v[nx].re*wave[dw].re - v[nx].im*wave[dw].im;
i1 = v[nx].im*wave[dw].re + v[nx].re*wave[dw].im;
r0 = v[0].re; i0 = v[0].im;
v[0].re = r0 + r1; v[0].im = i0 + i1;
v[nx].re = r0 - r1; v[nx].im = i0 - i1;
}
DFT_VecR2<T> vr2;
vr2(dst, c.n, n, dw0, wave);
}
else
{
DFT_R2<T> vr2;
vr2(dst, c.n, n, dw0, wave);
}
}
}
@ -813,94 +1057,21 @@ DFT(const OcvDftOptions & c, const Complex<T>* src, Complex<T>* dst)
if( factor == 3 )
{
// radix-3
for( i = 0; i < c.n; i += n )
if(c.haveSSE3)
{
Complex<T>* v = dst + i;
T r1 = v[nx].re + v[nx*2].re;
T i1 = v[nx].im + v[nx*2].im;
T r0 = v[0].re;
T i0 = v[0].im;
T r2 = sin_120*(v[nx].im - v[nx*2].im);
T i2 = sin_120*(v[nx*2].re - v[nx].re);
v[0].re = r0 + r1; v[0].im = i0 + i1;
r0 -= (T)0.5*r1; i0 -= (T)0.5*i1;
v[nx].re = r0 + r2; v[nx].im = i0 + i2;
v[nx*2].re = r0 - r2; v[nx*2].im = i0 - i2;
for( j = 1, dw = dw0; j < nx; j++, dw += dw0 )
{
v = dst + i + j;
r0 = v[nx].re*wave[dw].re - v[nx].im*wave[dw].im;
i0 = v[nx].re*wave[dw].im + v[nx].im*wave[dw].re;
i2 = v[nx*2].re*wave[dw*2].re - v[nx*2].im*wave[dw*2].im;
r2 = v[nx*2].re*wave[dw*2].im + v[nx*2].im*wave[dw*2].re;
r1 = r0 + i2; i1 = i0 + r2;
r2 = sin_120*(i0 - r2); i2 = sin_120*(i2 - r0);
r0 = v[0].re; i0 = v[0].im;
v[0].re = r0 + r1; v[0].im = i0 + i1;
r0 -= (T)0.5*r1; i0 -= (T)0.5*i1;
v[nx].re = r0 + r2; v[nx].im = i0 + i2;
v[nx*2].re = r0 - r2; v[nx*2].im = i0 - i2;
}
DFT_VecR3<T> vr3;
vr3(dst, c.n, n, dw0, wave);
}
else
{
DFT_R3<T> vr3;
vr3(dst, c.n, n, dw0, wave);
}
}
else if( factor == 5 )
{
// radix-5
for( i = 0; i < c.n; i += n )
{
for( j = 0, dw = 0; j < nx; j++, dw += dw0 )
{
Complex<T>* v0 = dst + i + j;
Complex<T>* v1 = v0 + nx*2;
Complex<T>* v2 = v1 + nx*2;
T r0, i0, r1, i1, r2, i2, r3, i3, r4, i4, r5, i5;
r3 = v0[nx].re*wave[dw].re - v0[nx].im*wave[dw].im;
i3 = v0[nx].re*wave[dw].im + v0[nx].im*wave[dw].re;
r2 = v2[0].re*wave[dw*4].re - v2[0].im*wave[dw*4].im;
i2 = v2[0].re*wave[dw*4].im + v2[0].im*wave[dw*4].re;
r1 = r3 + r2; i1 = i3 + i2;
r3 -= r2; i3 -= i2;
r4 = v1[nx].re*wave[dw*3].re - v1[nx].im*wave[dw*3].im;
i4 = v1[nx].re*wave[dw*3].im + v1[nx].im*wave[dw*3].re;
r0 = v1[0].re*wave[dw*2].re - v1[0].im*wave[dw*2].im;
i0 = v1[0].re*wave[dw*2].im + v1[0].im*wave[dw*2].re;
r2 = r4 + r0; i2 = i4 + i0;
r4 -= r0; i4 -= i0;
r0 = v0[0].re; i0 = v0[0].im;
r5 = r1 + r2; i5 = i1 + i2;
v0[0].re = r0 + r5; v0[0].im = i0 + i5;
r0 -= (T)0.25*r5; i0 -= (T)0.25*i5;
r1 = fft5_2*(r1 - r2); i1 = fft5_2*(i1 - i2);
r2 = -fft5_3*(i3 + i4); i2 = fft5_3*(r3 + r4);
i3 *= -fft5_5; r3 *= fft5_5;
i4 *= -fft5_4; r4 *= fft5_4;
r5 = r2 + i3; i5 = i2 + r3;
r2 -= i4; i2 -= r4;
r3 = r0 + r1; i3 = i0 + i1;
r0 -= r1; i0 -= i1;
v0[nx].re = r3 + r2; v0[nx].im = i3 + i2;
v2[0].re = r3 - r2; v2[0].im = i3 - i2;
v1[0].re = r0 + r5; v1[0].im = i0 + i5;
v1[nx].re = r0 - r5; v1[nx].im = i0 - i5;
}
}
DFT_R5<T> vr5;
vr5(dst, c.n, n, dw0, wave);
}
else
{

Write Preview
Loading…
Cancel
Save