| 2158 | |
| 2159 | |
| 2160 | static void |
| 2161 | DCTInit( int n, int elem_size, void* _wave, int inv ) |
| 2162 | { |
| 2163 | static const double DctScale[] = |
| 2164 | { |
| 2165 | 0.707106781186547570, 0.500000000000000000, 0.353553390593273790, |
| 2166 | 0.250000000000000000, 0.176776695296636890, 0.125000000000000000, |
| 2167 | 0.088388347648318447, 0.062500000000000000, 0.044194173824159223, |
| 2168 | 0.031250000000000000, 0.022097086912079612, 0.015625000000000000, |
| 2169 | 0.011048543456039806, 0.007812500000000000, 0.005524271728019903, |
| 2170 | 0.003906250000000000, 0.002762135864009952, 0.001953125000000000, |
| 2171 | 0.001381067932004976, 0.000976562500000000, 0.000690533966002488, |
| 2172 | 0.000488281250000000, 0.000345266983001244, 0.000244140625000000, |
| 2173 | 0.000172633491500622, 0.000122070312500000, 0.000086316745750311, |
| 2174 | 0.000061035156250000, 0.000043158372875155, 0.000030517578125000 |
| 2175 | }; |
| 2176 | |
| 2177 | int i; |
| 2178 | Complex<double> w, w1; |
| 2179 | double t, scale; |
| 2180 | |
| 2181 | if( n == 1 ) |
| 2182 | return; |
| 2183 | |
| 2184 | assert( (n&1) == 0 ); |
| 2185 | |
| 2186 | if( (n & (n - 1)) == 0 ) |
| 2187 | { |
| 2188 | int m; |
| 2189 | for( m = 0; (unsigned)(1 << m) < (unsigned)n; m++ ) |
| 2190 | ; |
| 2191 | scale = (!inv ? 2 : 1)*DctScale[m]; |
| 2192 | w1.re = DFTTab[m+2][0]; |
| 2193 | w1.im = -DFTTab[m+2][1]; |
| 2194 | } |
| 2195 | else |
| 2196 | { |
| 2197 | t = 1./(2*n); |
| 2198 | scale = (!inv ? 2 : 1)*std::sqrt(t); |
| 2199 | w1.im = sin(-CV_PI*t); |
| 2200 | w1.re = std::sqrt(1. - w1.im*w1.im); |
| 2201 | } |
| 2202 | n >>= 1; |
| 2203 | |
| 2204 | if( elem_size == sizeof(Complex<double>) ) |
| 2205 | { |
| 2206 | Complex<double>* wave = (Complex<double>*)_wave; |
| 2207 | |
| 2208 | w.re = scale; |
| 2209 | w.im = 0.; |
| 2210 | |
| 2211 | for( i = 0; i <= n; i++ ) |
| 2212 | { |
| 2213 | wave[i] = w; |
| 2214 | t = w.re*w1.re - w.im*w1.im; |
| 2215 | w.im = w.re*w1.im + w.im*w1.re; |
| 2216 | w.re = t; |
| 2217 | } |