本文不保证读者能都读懂。写本文的主要目的是测试数学公式功能,文字讲解基本没有,请谨慎阅读。
傅里叶变换
连续傅里叶变换是一个
它的作用是对于一个时域的复数函数,求出其频谱。
当自变量x表示时间(以秒为单位),变换变量ξ表示频率(以赫兹为单位)。在适当条件下,
它的作用是对于一个频域的复数函数,求出其时域表示。
离散傅里叶变换(DFT)
考虑将上述
以上两个变换有一种美妙的性质:都符合卷积定理,即可以将某一个域内的卷积对偶于其对应域内的乘法。
为了简化起见,可以将
为了方便起见,设符合方程
单位复根有以下性质:
那么,从矩阵角度考察
显然,以上两矩阵互逆。
考察矩阵1,
那么考虑,
对
即,原矩阵的这些位置
等于一个
这样就得到了两个
不断分治,就得到了一个
C++代码
typedef complex<double> T;
T dwg(int n, int flag)
{
return T(cos(pi * 2 / n), flag * sin(pi * 2 / n));
}
void fft(T *a, int m, int flag) // 求2^m点FFT
{
int nn = 1 << m;
while (m) {
int n2 = 1 << (m - 1);
int n = n2 << 1;
T wn(dwg(n, -flag));
for (int j = 0; j < nn; j += n) {
= 1;
T w *l = a + j, *r = a + n2 + j;
T *rb = a + n + j;
T ;
T tfor (; r < rb; ++l, ++r, w *= wn) {
= *l - *r;
t *l += *r;
*r = t * w;
}
}
--m;
}
for (int i = 0; i < n; ++i) { // 位反转
if (R[i] > i) {
(a[i], a[R[i]]);
swap}
}
}
对于正向FFT,取
例题
#include "bits/stdc++.h"
using namespace std;
const int BUF_SIZE = (int)1e6 + 10;
struct fastIO
{
char buf[BUF_SIZE];
int cur;
FILE *in, *out;
()
fastIO{
= BUF_SIZE;
cur = stdin;
in = stdout;
out }
inline char nextChar()
{
if (cur == BUF_SIZE) {
(buf, BUF_SIZE, 1, in);
fread= 0;
cur }
return buf[cur++];
}
inline int nextInt()
{
int x = 0;
char c = nextChar();
while (!('0' <= c && c <= '9')) c = nextChar();
while ('0' <= c && c <= '9') {
= x * 10 + c - '0';
x = nextChar();
c }
return x;
}
inline void printChar(char ch)
{
if (cur == BUF_SIZE) {
(buf, BUF_SIZE, 1, out);
fwrite= 0;
cur }
[cur++] = ch;
buf}
inline void printInt(int x)
{
if (x >= 10) printInt(x / 10);
(x % 10 + '0');
printChar}
inline void close()
{
if (cur > 0) {
(buf, cur, 1, out);
fwrite}
= 0;
cur }
} fin, fout;
const int lg2maxn = 18;
const int maxn = (1 << lg2maxn);
const double pi = acos(-1);
namespace nscpx
{
typedef double T;
struct cpx
{
, y;
T x(T a = 0, T b = 0):x(a), y(b) {}
cpxoperator-() const
cpx {
return cpx(-x, -y);
}
&operator+=(const cpx &a)
cpx {
+= a.x;
x += a.y;
y return *this;
}
&operator-=(const cpx &a)
cpx {
-= a.x;
x -= a.y;
y return *this;
}
&operator*=(const cpx &a)
cpx {
;
T nx= x * a.x - y * a.y;
nx = y * a.x + x * a.y;
y = nx;
x return *this;
}
operator+(const cpx &a) const
cpx {
return cpx(x + a.x, y + a.y);
}
operator-(const cpx &a) const
cpx {
return cpx(x - a.x, y - a.y);
}
operator*(const cpx &a) const
cpx {
return cpx(x * a.x - y * a.y, x * a.y + y * a.x);
}
};
}
using nscpx::cpx;
int R[maxn];
typedef cpx T;
T dwg(int n, int flag)
{
return T(cos(pi * 2 / n), flag * sin(pi * 2 / n));
}
void ffft(T *a, int m, int flag)
{
int nn = 1 << m;
while (m) {
int n2 = 1 << (m - 1);
int n = n2 << 1;
T wn(dwg(n, -flag));
for (int j = 0; j < nn; j += n) {
= 1;
T w *l = a + j, *r = a + n2 + j;
T *rb = a + n + j;
T ;
T tfor (; r < rb; ++l, ++r, w *= wn) {
= *l - *r;
t *l += *r;
*r = t * w;
}
}
--m;
}
}
void rev(T *a, int m)
{
int n = 1 << m;
for (int i = 0; i < n; ++i) {
if (R[i] > i) {
(a[i], a[R[i]]);
swap}
}
}
int n, m;
[maxn], b[maxn];
T a#define lowbit(x) ((x)&-(x))
int main()
{
= fin.nextInt();
n = fin.nextInt();
m int maxnn = m + n + 1;
while (maxnn != lowbit(maxnn)) maxnn += lowbit(maxnn);
int shit = 0;
while (maxnn > (1 << shit)) ++shit;
assert(maxnn == (1 << shit));
for (int i = 0; i < maxnn; ++i) {
[i] = (R[i >> 1] >> 1) | ((i & 1) << (shit - 1));
R}
for (int i = 0; i <= n; ++i) {
[i].x = fin.nextInt();
a}
for (int i = 0; i <= m; ++i) {
[i].x = fin.nextInt();
b}
(a, shit, 1);
ffft(b, shit, 1);
ffftfor (int i = 0; i < maxnn; ++i) {
[i] *= b[i];
a}
(a, shit);
rev(a, shit, -1);
ffft(a, shit);
rev.cur = 0;
foutfor (int i = 0; i <= (n + m); ++i) {
.printInt(a[i].x / maxnn + 0.5);
fout.printChar(' ');
fout}
.printChar('\n');
fout.close();
foutreturn 0;
}
钦此。