本文实例讲述了java实现的傅里叶变化算法。分享给大家供大家参考，具体如下：
用java实现傅里叶变化 结果为复数形式 a+bi
废话不多说，实现代码如下，共两个class
fft.class 傅里叶变化功能实现代码
?
|
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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
|
package fft.test;
/*************************************************************************
* compilation: javac fft.java execution: java fft n dependencies: complex.java
*
* compute the fft and inverse fft of a length n complex sequence. bare bones
* implementation that runs in o(n log n) time. our goal is to optimize the
* clarity of the code, rather than performance.
*
* limitations ----------- - assumes n is a power of 2
*
* - not the most memory efficient algorithm (because it uses an object type for
* representing complex numbers and because it re-allocates memory for the
* subarray, instead of doing in-place or reusing a single temporary array)
*
*************************************************************************/
public class fft {
// compute the fft of x[], assuming its length is a power of 2
public static complex[] fft(complex[] x) {
int n = x.length;
// base case
if (n == 1)
return new complex[] { x[0] };
// radix 2 cooley-tukey fft
if (n % 2 != 0) {
throw new runtimeexception("n is not a power of 2");
}
// fft of even terms
complex[] even = new complex[n / 2];
for (int k = 0; k < n / 2; k++) {
even[k] = x[2 * k];
}
complex[] q = fft(even);
// fft of odd terms
complex[] odd = even; // reuse the array
for (int k = 0; k < n / 2; k++) {
odd[k] = x[2 * k + 1];
}
complex[] r = fft(odd);
// combine
complex[] y = new complex[n];
for (int k = 0; k < n / 2; k++) {
double kth = -2 * k * math.pi / n;
complex wk = new complex(math.cos(kth), math.sin(kth));
y[k] = q[k].plus(wk.times(r[k]));
y[k + n / 2] = q[k].minus(wk.times(r[k]));
}
return y;
}
// compute the inverse fft of x[], assuming its length is a power of 2
public static complex[] ifft(complex[] x) {
int n = x.length;
complex[] y = new complex[n];
// take conjugate
for (int i = 0; i < n; i++) {
y[i] = x[i].conjugate();
}
// compute forward fft
y = fft(y);
// take conjugate again
for (int i = 0; i < n; i++) {
y[i] = y[i].conjugate();
}
// divide by n
for (int i = 0; i < n; i++) {
y[i] = y[i].scale(1.0 / n);
}
return y;
}
// compute the circular convolution of x and y
public static complex[] cconvolve(complex[] x, complex[] y) {
// should probably pad x and y with 0s so that they have same length
// and are powers of 2
if (x.length != y.length) {
throw new runtimeexception("dimensions don't agree");
}
int n = x.length;
// compute fft of each sequence，求值
complex[] a = fft(x);
complex[] b = fft(y);
// point-wise multiply，点值乘法
complex[] c = new complex[n];
for (int i = 0; i < n; i++) {
c[i] = a[i].times(b[i]);
}
// compute inverse fft，插值
return ifft(c);
}
// compute the linear convolution of x and y
public static complex[] convolve(complex[] x, complex[] y) {
complex zero = new complex(0, 0);
complex[] a = new complex[2 * x.length];// 2n次数界，高阶系数为0.
for (int i = 0; i < x.length; i++)
a[i] = x[i];
for (int i = x.length; i < 2 * x.length; i++)
a[i] = zero;
complex[] b = new complex[2 * y.length];
for (int i = 0; i < y.length; i++)
b[i] = y[i];
for (int i = y.length; i < 2 * y.length; i++)
b[i] = zero;
return cconvolve(a, b);
}
// display an array of complex numbers to standard output
public static void show(complex[] x, string title) {
system.out.println(title);
system.out.println("-------------------");
int complexlength = x.length;
for (int i = 0; i < complexlength; i++) {
// 输出复数
// system.out.println(x[i]);
// 输出幅值需要 * 2 / length
system.out.println(x[i].abs() * 2 / complexlength);
}
system.out.println();
}
/**
* 将数组数据重组成2的幂次方输出
*
* @param data
* @return
*/
public static double[] pow2doublearr(double[] data) {
// 创建新数组
double[] newdata = null;
int datalength = data.length;
int sumnum = 2;
while (sumnum < datalength) {
sumnum = sumnum * 2;
}
int addlength = sumnum - datalength;
if (addlength != 0) {
newdata = new double[sumnum];
system.arraycopy(data, 0, newdata, 0, datalength);
for (int i = datalength; i < sumnum; i++) {
newdata[i] = 0d;
}
} else {
newdata = data;
}
return newdata;
}
/**
* 去偏移量
*
* @param originalarr
* 原数组
* @return 目标数组
*/
public static double[] deskew(double[] originalarr) {
// 过滤不正确的参数
if (originalarr == null || originalarr.length <= 0) {
return null;
}
// 定义目标数组
double[] resarr = new double[originalarr.length];
// 求数组总和
double sum = 0d;
for (int i = 0; i < originalarr.length; i++) {
sum += originalarr[i];
}
// 求数组平均值
double aver = sum / originalarr.length;
// 去除偏移值
for (int i = 0; i < originalarr.length; i++) {
resarr[i] = originalarr[i] - aver;
}
return resarr;
}
public static void main(string[] args) {
// int n = integer.parseint(args[0]);
double[] data = { -0.35668879080953375, -0.6118094913035987, 0.8534269560320435, -0.6699697478438837, 0.35425500561437717,
0.8910250650549392, -0.025718699518642918, 0.07649691490732002 };
// 去除偏移量
data = deskew(data);
// 个数为2的幂次方
data = pow2doublearr(data);
int n = data.length;
system.out.println(n + "数组n中数量....");
complex[] x = new complex[n];
// original data
for (int i = 0; i < n; i++) {
// x[i] = new complex(i, 0);
// x[i] = new complex(-2 * math.random() + 1, 0);
x[i] = new complex(data[i], 0);
}
show(x, "x");
// fft of original data
complex[] y = fft(x);
show(y, "y = fft(x)");
// take inverse fft
complex[] z = ifft(y);
show(z, "z = ifft(y)");
// circular convolution of x with itself
complex[] c = cconvolve(x, x);
show(c, "c = cconvolve(x, x)");
// linear convolution of x with itself
complex[] d = convolve(x, x);
show(d, "d = convolve(x, x)");
}
}
/*********************************************************************
* % java fft 8 x ------------------- -0.35668879080953375 -0.6118094913035987
* 0.8534269560320435 -0.6699697478438837 0.35425500561437717 0.8910250650549392
* -0.025718699518642918 0.07649691490732002
*
* y = fft(x) ------------------- 0.5110172121330208 -1.245776663065442 +
* 0.7113504894129803i -0.8301420417085572 - 0.8726884066879042i
* -0.17611092978238008 + 2.4696418005143532i 1.1395317305034673
* -0.17611092978237974 - 2.4696418005143532i -0.8301420417085572 +
* 0.8726884066879042i -1.2457766630654419 - 0.7113504894129803i
*
* z = ifft(y) ------------------- -0.35668879080953375 -0.6118094913035987 +
* 4.2151962932466006e-17i 0.8534269560320435 - 2.691607282636124e-17i
* -0.6699697478438837 + 4.1114763914420734e-17i 0.35425500561437717
* 0.8910250650549392 - 6.887033953004965e-17i -0.025718699518642918 +
* 2.691607282636124e-17i 0.07649691490732002 - 1.4396387316837096e-17i
*
* c = cconvolve(x, x) ------------------- -1.0786973139009466 -
* 2.636779683484747e-16i 1.2327819138980782 + 2.2180047699856214e-17i
* 0.4386976685553382 - 1.3815636262919812e-17i -0.5579612069781844 +
* 1.9986455722517509e-16i 1.432390480003344 + 2.636779683484747e-16i
* -2.2165857430333684 + 2.2180047699856214e-17i -0.01255525669751989 +
* 1.3815636262919812e-17i 1.0230680492494633 - 2.4422465262488753e-16i
*
* d = convolve(x, x) ------------------- 0.12722689348916738 +
* 3.469446951953614e-17i 0.43645117531775324 - 2.78776395788635e-18i
* -0.2345048043334932 - 6.907818131459906e-18i -0.5663280251946803 +
* 5.829891518914417e-17i 1.2954076913348198 + 1.518836016779236e-16i
* -2.212650940696159 + 1.1090023849928107e-17i -0.018407034687857718 -
* 1.1306778366296569e-17i 1.023068049249463 - 9.435675069681485e-17i
* -1.205924207390114 - 2.983724378680108e-16i 0.796330738580325 +
* 2.4967811657742562e-17i 0.6732024728888314 - 6.907818131459906e-18i
* 0.00836681821649593 + 1.4156564203603091e-16i 0.1369827886685242 +
* 1.1179436667055108e-16i -0.00393480233720922 + 1.1090023849928107e-17i
* 0.005851777990337828 + 2.512241462921638e-17i 1.1102230246251565e-16 -
|