当前位置：首页 → 问答吧 → 成功了

成功了

时间：2014-03-22

来源：互联网

Discrete cosine transform type II~IV, Discrete sine transform type II~IV
快速算法(全部都是使用FFT)

复制内容到剪贴板代码:template <typename type, size_t Rank, typename Allocator>
multiarray<type, Rank, Allocator> &DCTII(const multiarray<type, Rank, Allocator> &buffer, multiarray<type, Rank, Allocator> &result)
{
typedef complex_type<type>::type complex_type;
typedef std::allocator_traits<Allocator>::rebind_alloc<complex_type>::other complex_alloc_type;
size_t N = buffer.size();
multiarray<complex_type, 1, complex_alloc_type> temp(N);
for(size_t i = 0; i < N / 2; ++i)
{
temp[i] = buffer[2 * i];
temp[N - i - 1] = buffer[2 * i + 1];
}
if(N % 2)
temp[N / 2] = buffer[N - 1];
Fourier(temp, temp);
result.reshape(N);
result[0] = temp[0].real();
for(size_t i = 1; i < N; ++i)
result[i] = std::real(temp[i] * std::polar<type>((type)M_SQRT2, -(type)M_PI_2 * (type)i / (type)N));
return result;
}
template <typename type, size_t Rank, typename Allocator>
multiarray<type, Rank, Allocator> &DCTIII(const multiarray<type, Rank, Allocator> &buffer, multiarray<type, Rank, Allocator> &result)
{
typedef complex_type<type>::type complex_type;
typedef std::allocator_traits<Allocator>::rebind_alloc<complex_type>::other complex_alloc_type;
size_t N = buffer.size();
multiarray<complex_type, 1, complex_alloc_type> temp(N);
temp[0] = (complex_type)buffer[0];
for(size_t i = 1; i < N; ++i)
temp[i] = std::polar<type>(buffer[i] * (type)M_SQRT2, -(type)M_PI_2 * (type)i / (type)N);
Fourier(temp, temp);
result.reshape(N);
for(size_t i = 0; i < N / 2; ++i)
{
result[2 * i] = temp[i].real();
result[2 * i + 1] = temp[N - i - 1].real();
}
if(N % 2)
result[N - 1] = temp[N / 2].real();
return result;
}
template <typename type, size_t Rank, typename Allocator>
multiarray<type, Rank, Allocator> &DCTIV(const multiarray<type, Rank, Allocator> &buffer, multiarray<type, Rank, Allocator> &result)
{
typedef complex_type<type>::type complex_type;
typedef std::allocator_traits<Allocator>::rebind_alloc<complex_type>::other complex_alloc_type;
size_t N = buffer.size();
multiarray<complex_type, 1, complex_alloc_type> temp(N);
for(size_t i = 0; i < N / 2; ++i)
{
temp[i] = buffer[2 * i];
temp[N - i - 1] = -buffer[2 * i + 1];
}
if(N % 2)
temp[N / 2] = buffer[N - 1];
for(size_t i = 1; i < N; ++i)
temp[i] *= std::polar<type>(1, -(type)M_PI * (type)i / (type)N);
Fourier(temp, temp);
result.reshape(N);
for(size_t i = 0; i < N; ++i)
result[i] = std::real(temp[i] * std::polar<type>((type)M_SQRT2, -(type)M_PI_4 * (type)(2 * i + 1) / (type)N));
return result;
}
template <typename type, size_t Rank, typename Allocator>
multiarray<type, Rank, Allocator> &DSTII(const multiarray<type, Rank, Allocator> &buffer, multiarray<type, Rank, Allocator> &result)
{
typedef complex_type<type>::type complex_type;
typedef std::allocator_traits<Allocator>::rebind_alloc<complex_type>::other complex_alloc_type;
size_t N = buffer.size();
multiarray<complex_type, 1, complex_alloc_type> temp(N);
for(size_t i = 0; i < N / 2; ++i)
{
temp[i] = -buffer[2 * i];
temp[N - i - 1] = buffer[2 * i + 1];
}
if(N % 2)
temp[N / 2] = -buffer[N - 1];
for(size_t i = 1; i < N; ++i)
temp[i] *= std::polar<type>(1, -(type)2 * (type)M_PI * (type)i / (type)N);
Fourier(temp, temp);
result.reshape(N);
for(size_t i = 0; i < N - 1; ++i)
result[i] = std::imag(temp[i] * std::polar<type>((type)M_SQRT2, -(type)M_PI_2 * (type)(i + 1) / (type)N));
result[N - 1] = std::imag(temp[N - 1] * std::polar<type>(1, -(type)M_PI_2));
return result;
}
template <typename type, size_t Rank, typename Allocator>
multiarray<type, Rank, Allocator> &DSTIII(const multiarray<type, Rank, Allocator> &buffer, multiarray<type, Rank, Allocator> &result)
{
typedef complex_type<type>::type complex_type;
typedef std::allocator_traits<Allocator>::rebind_alloc<complex_type>::other complex_alloc_type;
size_t N = buffer.size();
multiarray<complex_type, 1, complex_alloc_type> temp(N);
for(size_t i = 0; i < N - 1; ++i)
temp[i] = std::polar<type>(buffer[i] * (type)M_SQRT2, -(type)M_PI_2 * (type)(i + 1) / (type)N);
temp[N - 1] = std::polar<type>(buffer[N - 1], -(type)M_PI_2);
Fourier(temp, temp);
for(size_t i = 1; i < N; ++i)
temp[i] *= std::polar<type>(1, -(type)2 * (type)M_PI * (type)i / (type)N);
result.reshape(N);
for(size_t i = 0; i < N / 2; ++i)
{
result[2 * i] = -temp[i].imag();
result[2 * i + 1] = temp[N - i - 1].imag();
}
if(N % 2)
result[N - 1] = -temp[N / 2].imag();
return result;
}
template <typename type, size_t Rank, typename Allocator>
multiarray<type, Rank, Allocator> &DSTIV(const multiarray<type, Rank, Allocator> &buffer, multiarray<type, Rank, Allocator> &result)
{
typedef complex_type<type>::type complex_type;
typedef std::allocator_traits<Allocator>::rebind_alloc<complex_type>::other complex_alloc_type;
size_t N = buffer.size();
multiarray<complex_type, 1, complex_alloc_type> temp(N);
for(size_t i = 0; i < N / 2; ++i)
{
temp[i] = -buffer[2 * i];
temp[N - i - 1] = -buffer[2 * i + 1];
}
if(N % 2)
temp[N / 2] = -buffer[N - 1];
for(size_t i = 1; i < result.size(0); ++i)
temp[i] *= std::polar<type>(1, -(type)M_PI * (type)i / (type)N);
Fourier(temp, temp);
result.reshape(N);
for(size_t i = 0; i < N; ++i)
result[i] = std::imag(temp[i] * std::polar<type>((type)M_SQRT2, -(type)M_PI_4 * (type)(2 * i + 1) / (type)N));
return result;
}

作者: Susan﹏汪汪发布时间: 2014-03-22

复制内容到剪贴板代码:multiarray<double> a(GetRandom());

generate(a.begin(), a.end(), GetRandom);
show(a);
cout << endl;

cout << "DFT: " << endl;
for(size_t i = 0; i < a.size(); ++i)
cout << round(abs(DiscreteFourier(a, i)) * 100) / 100 << '\t';
cout << endl;
for(size_t i = 0; i < a.size(); ++i)
cout << round(arg(DiscreteFourier(a, i)) * 100) / 100 << '\t';
cout << endl << endl;

cout << "FFT: " << endl;
multiarray<complex<double>> temp;
Fourier(a, temp);
for(size_t i = 0; i < temp.size(); ++i)
cout << round(abs(temp.at(i)) * 100) / 100 << '\t';
cout << endl;
for(size_t i = 0; i < temp.size(); ++i)
cout << round(arg(temp.at(i)) * 100) / 100 << '\t';
cout << endl << endl;

cout << "IFFT: " << endl;
Flip(Fourier(temp, temp), temp);
for(size_t i = 0; i < temp.size(); ++i)
cout << round(abs(temp.at(i)) * 100) / 100 << '\t';
cout << endl;
for(size_t i = 0; i < temp.size(); ++i)
cout << round(arg(temp.at(i)) * 100) / 100 << '\t';
cout << endl << endl;

cout << "DCTII: " << endl;
show(DCTII(a, a));
cout << "DCTIII: " << endl;
show(DCTIII(a, a));
cout << "DCTIV: " << endl;
show(DCTIV(a, a));
cout << "DCTIV: " << endl;
show(DCTIV(a, a));

cout << "DSTII: " << endl;
show(DSTII(a, a));
cout << "DSTIII: " << endl;
show(DSTIII(a, a));
cout << "DSTIV: " << endl;
show(DSTIV(a, a));
cout << "DSTIV: " << endl;
show(DSTIV(a, a));

未命名.png (11.14 KB)

2014-1-11 12:52 AM