unsigned int bitReverse(unsigned int x, int log2n) {
int n = 0;
int mask = 0x1;
for (int i=0; i < log2n; i++) {
n <<= 1;
n |= (x & 1);
x >>= 1;
}
return n;
}
const double PI = 3.1415926536;
template<class Iter_T>
void fft(Iter_r a, Iter_r b, int log2n) {
typedef typename iterator_traits<Iter T>::value_type complex;
const complex J(0, 1);
int n = 1 << log2n;
for (unsigned int i=0; i < n; ++i) {
b[bitReverse(i, log2n)] = a[i];
}
for (int s = 1; s <= log2n; ++s) {
int m = 1 << s;
int m2 = m >> 1;
complex w(1, 0);
complex wm = exp(-J * (PI / m2));
for (int j=0; j < m2; ++j) {
for (int k=j; k < n; k += m) {
complex t = w * b[k + m2];
complex u = b[k];
b[k] = u + t;
b[k + m2] = u - t;
}
w *= wm;
}
}
}
int main() {
typedef complex<double> cx;
cx a[] = { cx(0, 0), cx(1, 1), cx(3, 3), cx(4, 4),
cx(4, 4), cx(3, 3), cx(1, 1), cx(0, 0) };
cx b[8];
fft(a, b, 3);
for (int i=0; i<8; ++i) cout << b[i] << "\n";
}
Программа примера 11.33 выдает следующий результат.
(16,16)
(-4.82843,-11.6569)
(0,0)
(-0.343146,0.828427)
(0.0)
(0.828427,-0.343146)
(0,0)
(-11.6569,-4.82843)
Преобразование Фурье играет важную роль в спектральном анализе и часто используется в технических и научных приложениях. БПФ — это алгоритм вычисления ДПФ, который имеет сложность порядка N log2(N) в отличие от ожидаемой сложности N² для простой реализации ДПФ. Такое впечатляющее ускорение достигается в БПФ благодаря устранению избыточных вычислений.
Очень не просто найти хорошую реализацию БПФ, написанную на «правильном» C++ (т. е. когда программа на C++ не является механическим переложением алгоритмов, написанных на Фортране или С) и которая не была бы защищена сильно ограничивающей лицензией. Представленный в примере 11.33 программный код основан на открытом коде, который можно найти в сетевой конференции Usenet, посвященной цифровой обработке сигналов (comp.dsp). Большим преимуществом реализации БПФ на правильном C++ по сравнению с более распространенным решением в стиле С является то, что стандартная библиотека содержит шаблон complex, который позволяет существенно снизить объем необходимого программного кода. В представленной в примере 11.33 функции fft() основное внимание уделялось простоте, а не эффективности.
11.18. Работа с полярными координатами
Требуется обеспечить представление полярных координат и манипулирование ими.
Шаблон complex из заголовочного файла <complex> содержит функции преобразования в полярные координаты и обратно. Пример 11.34 показывает, как можно использовать класс шаблона complex для представления и манипулирования полярными координатами.
Пример 11.34. Применение шаблонного класса complex для представления полярных координат
#include <complex>
#include <iostream>
using namespace std;
int main() {
double rho = 3.0; // длина
double theta = 3.141592 / 2; // угол
complex<double> coord = polar(rho, theta);
cout << "rho = " << abs(coord) << ", theta = " << arg(coord) << endl;
coord += polar(4.0, 0.0);
cout << "rho = " << abs(coord) << ", theta = " << arg(coord) << endl;
}
Программа примера 11.34 выдает следующий результат.
rho = 3, theta = 1.5708
rho = 5, theta = 0.643501
Существует естественная связь между полярными координатами и комплексными числами. Хотя эти понятия в какой-то мере взаимозаменяемы, использование одного и того же типа для представления разных концепций в целом нельзя считать хорошей идеей. Поскольку применение шаблона complex для представления полярных координат не является элегантным решением, я предусмотрел приведенный в примере 11.25 класс полярных координат, допускающий более естественное применение.
Пример 11.35. Класс полярных координат
#include <complex>
#include <iostream>
using namespace std;
template<class T>
struct BasicPolar {
public typedef BasicPolar self;
// конструкторы
BasicPolar() : m() {} BasicPolar(const self& x) : m(x.m) {}
BasicPolar(const T& rho, const T& theta) : m(polar(rho, theta)) {}
// операторы присваивания
self operator-() { return Polar(-m); }
self& operator+=(const self& x) { m += x.m; return *this; }
self& operator-=(const self& x) { m -= x.m; return *this; }
self& operator*=(const self& x) { m *= x.m; return *this; }
self& operator/=(const self& x) { m /= x.m; return *this; }
operator complex<T>() const { return m; }
// открытые функции-члены
T rho() const { return abs(m); }
T theta() const { return arg(m); }
// бинарные операции
friend self operator+(self x, const self& y) { return x += y; }
friend self operator-(self x, const self& y) { return x -= y; }
friend self operator*(self x, const self& y) { return x *= y; }
friend self operator/(self x, const self& y) { return x /= y; }
// операторы сравнения
friend bool operator==(const self& x, const self& y) { return x.m == y.m; }
friend bool operator!=(const self& x, const self& y) { return x.m ! = y.m; }
private:
complex<T> m;
};
typedef BasicPolar<double> Polar;
int main() {