|
2 | 2 | from cmath import exp, pi
|
3 | 3 | from math import log2
|
4 | 4 |
|
| 5 | + |
5 | 6 | def dft(X):
|
6 | 7 | N = len(X)
|
7 |
| - temp = [0]*N |
| 8 | + temp = [0] * N |
8 | 9 | for i in range(N):
|
9 | 10 | for k in range(N):
|
10 |
| - temp[i] += X[k] * exp(-2.0j*pi*i*k/N) |
| 11 | + temp[i] += X[k] * exp(-2.0j * pi * i * k / N) |
11 | 12 | return temp
|
12 | 13 |
|
| 14 | + |
13 | 15 | def cooley_tukey(X):
|
14 |
| - N = len(X) |
15 |
| - if N <= 1: |
16 |
| - return X |
17 |
| - even = cooley_tukey(X[0::2]) |
18 |
| - odd = cooley_tukey(X[1::2]) |
19 |
| - |
20 |
| - temp = [i for i in range(N)] |
21 |
| - for k in range(N//2): |
22 |
| - temp[k] = even[k] + exp(-2j*pi*k/N) * odd[k] |
23 |
| - temp[k+N//2] = even[k] - exp(-2j*pi*k/N)*odd[k] |
24 |
| - return temp |
25 |
| - |
26 |
| -def bitReverse(X): |
27 |
| - N = len(X) |
28 |
| - temp = [i for i in range(N)] |
29 |
| - for k in range(N): |
30 |
| - b = sum(1<<(int(log2(N))-1-i) for i in range(int(log2(N))) if k>>i&1) |
31 |
| - temp[k] = X[b] |
32 |
| - temp[b] = X[k] |
33 |
| - return temp |
| 16 | + N = len(X) |
| 17 | + if N <= 1: |
| 18 | + return X |
| 19 | + even = cooley_tukey(X[0::2]) |
| 20 | + odd = cooley_tukey(X[1::2]) |
| 21 | + |
| 22 | + temp = [i for i in range(N)] |
| 23 | + for k in range(N // 2): |
| 24 | + temp[k] = even[k] + exp(-2.0j * pi * k / N) * odd[k] |
| 25 | + temp[k + N // 2] = even[k] - exp(-2.0j * pi * k / N) * odd[k] |
| 26 | + return temp |
| 27 | + |
| 28 | + |
| 29 | +def bit_reverse(X): |
| 30 | + N = len(X) |
| 31 | + temp = [i for i in range(N)] |
| 32 | + for k in range(N): |
| 33 | + b = sum(1 << int(log2(N)) - 1 - |
| 34 | + i for i in range(int(log2(N))) if k >> i & 1) |
| 35 | + temp[k] = X[b] |
| 36 | + temp[b] = X[k] |
| 37 | + return temp |
| 38 | + |
34 | 39 |
|
35 | 40 | def iterative_cooley_tukey(X):
|
36 |
| - N = len(X) |
| 41 | + N = len(X) |
| 42 | + |
| 43 | + X = bit_reverse(X) |
37 | 44 |
|
38 |
| - X = bitReverse(X) |
| 45 | + for i in range(1, int(log2(N)) + 1): |
| 46 | + stride = 2 ** i |
| 47 | + w = exp(-2.0j * pi / stride) |
| 48 | + for j in range(0, N, stride): |
| 49 | + v = 1 |
| 50 | + for k in range(stride // 2): |
| 51 | + X[k + j + stride // 2] = X[k + j] - v * X[k + j + stride // 2] |
| 52 | + X[k + j] -= X[k + j + stride // 2] - X[k + j] |
| 53 | + v *= w |
| 54 | + return X |
39 | 55 |
|
40 |
| - for i in range(1, int(log2(N)) + 1): |
41 |
| - stride = 2**i |
42 |
| - w = exp(-2j*pi/stride) |
43 |
| - for j in range(0, N, stride): |
44 |
| - v = 1 |
45 |
| - for k in range(stride//2): |
46 |
| - X[k + j + stride//2] = X[k + j] - v*X[k + j + stride//2]; |
47 |
| - X[k + j] -= (X[k + j + stride//2] - X[k + j]); |
48 |
| - v *= w; |
49 |
| - return X |
50 | 56 |
|
51 | 57 | X = []
|
52 | 58 |
|
53 | 59 | for i in range(64):
|
54 |
| - X.append(random()) |
| 60 | + X.append(random()) |
55 | 61 |
|
56 | 62 | Y = cooley_tukey(X)
|
57 | 63 | Z = iterative_cooley_tukey(X)
|
|
0 commit comments