魔法师の高塔

快速傅里叶变换

快速傅里叶变换具体原理参阅『算法导论』。

用rust实现如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
use std::ops::{Add, AddAssign, Mul, Div};
use std::f32;
static PI: f32 = 3.1415926535;
#[derive(Debug,Clone,Copy)]
struct Complex {
real: f32,
image: f32,
}
impl Complex {
fn new(real: f32, image: f32) -> Self {
Complex {
real: real,
image: image,
}
}
fn quick_pow(&self, n: usize) -> Self {
match n {
0 => Complex::new(1f32, 0f32),
1 => *self,
_ if n % 2 == 0 => {
let a = self.quick_pow(n >> 1);
a * a
}
_ => {
let a = self.quick_pow(n - 1);
(*self) * a
}
}
}
}
impl Add for Complex {
type Output = Complex;
fn add(self, comp: Complex) -> Complex {
Complex {
real: self.real + comp.real,
image: self.image + comp.image,
}
}
}
impl AddAssign for Complex {
fn add_assign(&mut self, comp: Complex) {
self.real = self.real + comp.real;
self.image = self.image + comp.image;
}
}
impl Mul for Complex {
type Output = Complex;
fn mul(self, comp: Complex) -> Complex {
Complex {
real: self.real * comp.real - self.image * comp.image,
image: self.real * comp.image + self.image * comp.real,
}
}
}
impl Div for Complex {
type Output = Complex;
fn div(self, comp: Complex) -> Complex {
let distance = comp.real * comp.real + comp.image * comp.image;
Complex {
real: (self.real * comp.real + self.image * comp.image) / distance,
image: (self.image * comp.real - self.real * comp.image) / distance,
}
}
}
fn dft(num: Vec<Complex>, res: bool) -> Vec<Complex> {
let mut root_vec: Vec<Complex> = Vec::new();
for i in 0..num.len() {
let value: f32 = PI / 4f32 * (i as f32);
root_vec.push(Complex::new(value.cos(), value.sin()))
}
let mut result: Vec<Complex> = Vec::new();
for i in 0..num.len() {
let mut x: Complex = Complex::new(0f32, 0f32);
for j in 0..num.len() {
if res == true {
x += num[j] / root_vec[j].quick_pow(i) / Complex::new(num.len() as f32, 0f32)
} else {
x += num[j] * root_vec[j].quick_pow(i)
}
}
result.push(x)
}
result
}
fn vecf322veccomplex(num: Vec<f32>) -> Vec<Complex> {
let mut result: Vec<Complex> = Vec::new();
for i in 0..num.len() {
result.push(Complex::new(num[i], 0f32))
}
result
}
fn mul_vec_complex(num1: Vec<Complex>, num2: Vec<Complex>) -> Vec<Complex> {
let mut result = Vec::new();
for i in 0..num1.len() {
result.push(num1[i] * num2[i])
}
result
}
fn add_vec_complex(num1: Vec<Complex>, num2: Vec<Complex>) -> Vec<Complex> {
let mut result = Vec::new();
for i in 0..num1.len() {
result.push(num1[i] + num2[i])
}
result
}
fn dft2vec(num: Vec<Complex>) -> Vec<f32> {
let mut result1 = Vec::new();
for value in num {
result1.push(value.real.round())
}
let mut result = Vec::new();
for i in 0..result1.len() {
if result1[i] > 10.0 {
let a = (result1[i] / 10.0).trunc();
result.push(result1[i] - a * 10.0);
result1[i + 1] += a;
} else {
result.push(result1[i])
}
}
result
}
fn main() {
let do1 = vec![1.0, 3.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0];
let do2 = vec![3.0, 2.0, 4.0, 3.0, 2.0, 0.0, 0.0, 0.0];
let dft1 = dft(vecf322veccomplex(do1), false);
let dft2 = dft(vecf322veccomplex(do2), false);
let reserve_dft = dft(mul_vec_complex(dft1, dft2), true);
println!("{:?}", 331 * 23423);
println!("{:?}", reserve_dft);
println!("{:?}", dft2vec(reserve_dft));
}

快速傅里叶变换

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
fn fft(num: Vec<Complex>, res: bool) -> Vec<Complex> {
if num.len() == 1 {
num
} else {
let w: f32 = 2f32 * PI / (num.len() as f32);
let root = Complex::new(w.cos(), w.sin());
let mut a0 = Vec::new();
let mut a1 = Vec::new();
let mut result = Vec::new();
for _ in 0..num.len() {
result.push(Complex::new(0f32, 0f32));
}
for i in 0..num.len() {
if i % 2 == 0 {
a0.push(num[i]);
} else {
a1.push(num[i]);
}
}
a0 = fft(a0, res);
a1 = fft(a1, res);
for i in 0..num.len() / 2 {
if res == true {
let val = Complex::new(0.5, 0f32);
result[i] += val * (a0[i] + a1[i] / root.quick_pow(i));
result[i + num.len() / 2] += val * (a0[i] - a1[i] / root.quick_pow(i));
} else {
result[i] += a0[i] + a1[i] * root.quick_pow(i);
result[i + num.len() / 2] += a0[i] - a1[i] * root.quick_pow(i);
}
}
result
}
}