【注意】最后更新于 December 9, 2015，文中内容可能已过时，请谨慎使用。
昨天晚上数字信号处理粗略的复习了下，感觉快速傅里叶变换对于使用挺重要的，就想着编程实现它。
今天上午没课，想着一上午实现吧，结果到晚上才编出来。
反正是出结果了，我也不管结果对不对了，反正原理算是差不多了解了。
首先，什么是傅里叶变换，傅里叶变换最基本的就不赘述了。
主要说下傅里叶变换的速度以及快速傅里叶变换的原理。
我用的方法是时间抽选基2的算法（库里-图基算法）。
这个遵循时间偶奇分解，频率前后分解的规则。
核心思想是利用旋转因子的对称性和周期性，将其变成小的DFT来计算，降低运算复杂度。
最容易理解的就是蝶形图，下面的是蝶形图.
嗯，大体上就是上面这个图。
由于在进行运算的时候要进行偶奇分开，就得把正常的时间序列转换下。
主要方法是变址计算。比如8进制的如下表
| 自然顺序 |
二进制 |
码位倒置 |
倒置顺序 |
| 0 |
000 |
000 |
0 |
| 1 |
001 |
100 |
4 |
| 2 |
010 |
010 |
2 |
| 3 |
011 |
110 |
6 |
| 4 |
100 |
001 |
1 |
| ... |
... |
... |
... |
就是这个意思，不打啦。变址之后就进行正常的傅里叶计算，从最小的开始，逐步扩大。
这个地方有点像递归计算。不过递归的效率只能呵呵了。
在不同级进行傅里叶变换的长度间隔是一个具体的表达式，不愿意写啦，有个参考的写得挺详细。
再说下，代码我也没检查是不是执行正确的，就是这么个思想。
实在是不想再改了，数字信号处理复习的一塌糊涂，得复习了。
这个是主要的代码
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
|
public class FFT {
public FFT() {
}
public static final int N = 8;
public static final double PI = Math.PI;
public static void main(String[] args) {
FFT fft = new FFT();
Complex time[] = new Complex[N];
//创建个时域的序列
for (int i = 0; i < time.length; i++) {
Complex com = new Complex();
com.re = Math.random() * 10;
com.im = Math.random() * 10;
time[i] = com;
}
fft.listNum(time);
//进行二进制数倒置
fft.changeOrder(time);
fft.listNum(time);
//进行快速傅里叶变换
fft.trans(time);
fft.listNum(time);
}
private void trans(Complex[] time) {
int length = time.length;
int m = log2(length);
for (int i = 0; i < m; i++) {
//间隔是2的i次幂,也代表了有几个数参与运算
int l = (int) Math.pow(2, i);
//在第i次有几组进行运算
int g = length / (l * 2), t = 0;
for (int j = 0; j < g; j++) {
//计算每一组的变换值
for (; t < 2 * l * (j + 1); t++) {
//前l个是加运算，后l是减运算
if (t < l * (j + 1)) {
Complex com = time[t];
Complex product = comMultiply(time[t + l], getWN(t));
time[t] = comPlus(com, product);
} else {
Complex com = time[t - l];
Complex product = comMultiply(time[t], getWN(t));
time[t] = comMinus(com, product);
}
}
}
}
}
public void changeOrder(Complex time[]) {
int length = time.length;
//经过m次分解
int m = log2(length);
int k, tempNum;
for (int i = 0; i < length; i++) {
Complex com = new Complex();
k = i;
tempNum = 0;
//进行二进制倒置,m就代表倒置几次
for (int t = 0; t < m; t++) {
//将低位的乘到高位
tempNum *= 2;
//取出k的低位，通过循环回来乘到高位
tempNum += k % 2;
//将已经取出来的数删除
k /= 2;
}
//完成调换
if (tempNum > i) {
com = time[i];
time[i] = time[tempNum];
time[tempNum] = com;
}
}
}
public void listNum(Complex time[]) {
System.out.println("结果如下");
for (Complex com : time) {
System.out.println(com.re + "+ j" + com.im);
}
}
public int log2(int value) {
return (int) (Math.log(value) / Math.log(2));
}
public Complex comPlus(Complex com1, Complex com2) {
Complex com = new Complex();
com.re = com1.re + com2.re;
com.im = com1.im + com2.im;
return com;
}
public Complex comMinus(Complex com1, Complex com2) {
Complex com = new Complex();
com.re = com1.re - com2.re;
com.im = com1.im - com2.im;
return com;
}
public Complex comMultiply(Complex com1, Complex com2) {
Complex com = new Complex();
com.re = com1.re * com2.re - com1.im * com2.im;
com.im = com1.re * com2.im + com2.re * com1.im;
return com;
}
public Complex getWN(int n) {
Complex com = new Complex();
com.re = Math.cos(2 * PI / N * n);
com.im = Math.sin(2 * PI / N * n) * (-1);
return com;
}
}
|
我把复数写成一个单独的类了
1
2
3
4
5
6
7
8
9
|
public class Complex {
double re,im;
public Complex() {
}
public Complex(double re, double im) {
this.re = re;
this.im = im;
}
}
|
关于代码主要比较晦涩的地方我也注了注释，希望对理解有帮助。
如果还是看不懂，那就没法了，反正我写这段代码的时候也是晕头转向的。
说实话用这种面向对象的语言来编写这面向过程的程序，除了真爱想不到其他的了。
参考文献
1.学步园
2.百度文库
3.最重要的还是课本啦！
