FFT

Posted __AiR_H

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了FFT相关的知识,希望对你有一定的参考价值。


 

1 HDU 1402 A * B Problem Plus

技术分享
 1 #include <iostream>
 2 #include <cstdio>
 3 #include <cstring>
 4 #include <cstdlib>
 5 #include <algorithm>
 6 #include <queue>
 7 #include <vector>
 8 #include <stack>
 9 #include <map>
10 #include <set>
11 #include <cmath>
12 #include <cctype>
13 #include <ctime>
14 #include <ext/rope>
15 
16 using namespace std;
17 using namespace __gnu_cxx;
18 
19 #define REP(i, n) for (int i = 0; i < (n); ++i)
20 #define eps 1e-9
21 #define SZ(x) ((int)x.size())
22 
23 typedef long long ll;
24 typedef pair<int, int> pii;
25 
26 #define PI acos(-1.0)
27 const int maxn = 1e5 + 10;
28 struct Complex {
29     double x, y;
30     Complex(double _x = 0.0, double _y = 0.0) { x = _x; y = _y; }
31     Complex operator - (const Complex &b) const { return Complex(x - b.x, y - b.y); }
32     Complex operator + (const Complex &b) const { return Complex(x + b.x, y + b.y); }
33     Complex operator * (const Complex &b) const { return Complex(x * b.x - y * b.y, x * b.y + y * b.x); }
34 };
35 char str_A[maxn], str_B[maxn];
36 Complex A[maxn * 4], B[maxn * 4], ans[maxn * 4];
37 int ans_final[maxn * 4], len_final = 0;
38 int N, M, n = 1;
39 
40 void fft(Complex a[], int len, int oper) {
41     for (int i = 1, j = 0; i < len - 1; ++i) {
42         for (int s = len; ~j & s;) { j ^= s >>= 1; }
43         if (i < j) { swap(a[i], a[j]); }
44     }
45     Complex omega_m, omega, t, u;
46     int m, m2;
47     for (int d = 0; (1 << d) < len; ++d) {
48         m = (1 << d); m2 = m << 1;
49         omega_m = Complex(cos(PI / m * oper), sin(PI / m * oper));
50         for (int i = 0; i < len; i += m2) {
51             omega = Complex(1, 0);
52             for (int j = 0; j < m; ++j) {
53                 u = a[i + j]; t = omega * a[i + j + m];
54                 a[i + j] = u + t; a[i + j + m] = u - t;
55                 omega = omega * omega_m;
56             }
57         }
58     }
59     if (oper == -1) { REP(i, N + M - 1) { a[i].x /= len; } }
60 }
61 
62 
63 int main() {
64 #ifdef __AiR_H
65     freopen("in.txt", "r", stdin);
66 //    freopen("out.txt", "w", stdout);
67 #endif // __AiR_H
68     while (scanf("%s %s", str_A, str_B) != EOF) {
69         memset(ans, 0, sizeof(ans)); memset(ans_final, 0, sizeof(ans_final));
70         N = strlen(str_A); M = strlen(str_B); n = 1;
71         while (n < N + M - 1) { n <<= 1; }
72         REP(i, N) { A[i] = Complex(str_A[N - 1 - i] - 0, 0); }
73         for (int i = N; i < n; ++i) { A[i] = Complex(0, 0); }
74         REP(i, M) { B[i] = Complex(str_B[M - 1 - i] - 0, 0); }
75         for (int i = M; i < n; ++i) { B[i] = Complex(0, 0); }
76         fft(A, n, 1); fft(B, n, 1);
77         REP(i, n) { ans[i] = A[i] * B[i]; }
78         fft(ans, n, -1);
79         REP(i, N + M) { ans_final[i] = ans[i].x + 0.5; }
80         REP(i, N + M - 1) { ans_final[i + 1] += ans_final[i] / 10; ans_final[i] %= 10; }
81         len_final = N + M - 1;
82         if (ans_final[len_final]) { ++len_final; }
83         while (ans_final[len_final - 1] == 0 && len_final > 1) { --len_final; }
84         for (int i = len_final - 1; i > 0; --i) { printf("%d", ans_final[i]); } printf("%d\n", ans_final[0]);
85     }
86     return 0;
87 }
View Code

 

 

 

 


 

以上是关于FFT的主要内容,如果未能解决你的问题,请参考以下文章

为啥不使用 Scipy 的 FFT 代码中的结果与 Scipy FFT 不相似?

基数 32 FFT 实现

FFT的大小实际上是啥意思

数字信号处理3: 快速傅里叶变换(FFT)(含代码)

数字信号处理3: 快速傅里叶变换(FFT)(含代码)

FFT频谱分析(补零频谱泄露栅栏效应加窗细化频谱混叠),MatlabC语言代码