1- import numpy as np
1+ """
2+ Fourier algorithm related utils
3+ """
24
5+ import numpy as np
36
4- def DFT_slow (x ):
5- """ Compute the discrete Fourier Transform of the one-dimensional array x"""
6- x = np .asarray (x )
7- N = x .shape [0 ]
8- n = np .arange (N )
9- k = n .reshape ((N , 1 ))
10- M = np .exp (- 2j * np .pi * k * n / N )
11- return np .dot (M , x )
127
8+ def dft_slow (input_data ):
9+ """
10+ Compute the discrete Fourier Transform of the one-dimensional array x
11+ """
12+ input_data = np .asarray (input_data )
13+ data_length = input_data .shape [0 ]
14+ data_sorted = np .arange (data_length )
15+ data_array = data_sorted .reshape ((data_length , 1 ))
16+ exp_data = np .exp (- 2j * np .pi * data_array * data_sorted / data_length )
17+ return np .dot (exp_data , input_data )
1318
14- def FFT_vectorized (x ):
15- """A vectorized, non-recursive version of the Cooley-Tukey FFT"""
16- x = np .asarray (x )
17- N = x .shape [0 ]
1819
19- if np .log2 (N ) % 1 > 0 :
20- raise ValueError ("size of x must be a power of 2" )
20+ def fft_vectorized (input_data ):
21+ """
22+ A vectorized, non-recursive version of the Cooley-Tukey FFT
23+ """
24+ input_data = np .asarray (input_data )
25+ data_length = input_data .shape [0 ]
26+
27+ if np .log2 (data_length ) % 1 > 0 :
28+ raise ValueError ("size of input data must be a power of 2" )
2129
2230 # N_min here is equivalent to the stopping condition above,
2331 # and should be a power of 2
24- N_min = min (N , 32 )
25-
26- # Perform an O[N^2] DFT on all length-N_min sub-problems at once
27- n = np .arange (N_min )
28- k = n [:, None ]
29- M = np .exp (- 2j * np .pi * n * k / N_min )
30- X = np .dot (M , x .reshape ((N_min , - 1 )))
32+ min_data = min (data_length , 32 )
33+
34+ # Perform an O[N^2] DFT on all length-min_data sub-problems at once
35+ data_sorted = np .arange (min_data )
36+ data_matrix = data_sorted [:, None ]
37+ exp_data = np .exp (- 2j * np .pi * data_sorted * data_matrix / min_data )
38+ dot_product = np .dot (exp_data , input_data .reshape ((min_data , - 1 )))
3139
3240 # build-up each level of the recursive calculation all at once
33- while X .shape [0 ] < N :
34- X_even = X [:, :X .shape [1 ] // 2 ]
35- X_odd = X [:, X .shape [1 ] // 2 :]
36- factor = np .exp (- 1j * np .pi * np .arange (X .shape [0 ])
37- / X .shape [0 ])[:, None ]
38- X = np .vstack ([X_even + factor * X_odd ,
39- X_even - factor * X_odd ])
41+ while dot_product .shape [0 ] < data_length :
42+ data_even = dot_product [:, :dot_product .shape [1 ] // 2 ]
43+ data_odd = dot_product [:, dot_product .shape [1 ] // 2 :]
44+ factor = np .exp (- 1j * np .pi * np .arange (dot_product .shape [0 ])
45+ / dot_product .shape [0 ])[:, None ]
46+ dot_product = np .vstack ([data_even + factor * data_odd ,
47+ data_even - factor * data_odd ])
4048
41- return X .ravel ()
49+ return dot_product .ravel ()
4250
4351
44- def fftfreq ( n , d = 1.0 ):
52+ def fft_freq ( window_len , spacing = 1.0 ):
4553 """
4654 Return the Discrete Fourier Transform sample frequencies.
4755
@@ -75,14 +83,13 @@ def fftfreq(n, d=1.0):
7583 >>> freq = np.fft.fftfreq(n, d=timestep)
7684 >>> freq
7785 array([ 0. , 1.25, 2.5 , 3.75, -5. , -3.75, -2.5 , -1.25])
78-
7986 """
80-
81- val = 1.0 / (n * d )
82- results = np .empty (n , int )
83- N = (n - 1 )// 2 + 1
84- p1 = np .arange (0 , N , dtype = int )
85- results [:N ] = p1
86- p2 = np .arange (- (n // 2 ), 0 , dtype = int )
87- results [N :] = p2
88- return results * val
87+
88+ val = 1.0 / (window_len * spacing )
89+ results = np .empty (window_len , int )
90+ window_half = (window_len - 1 )// 2 + 1
91+ window_p1 = np .arange (0 , window_half , dtype = int )
92+ results [:window_half ] = window_p1
93+ window_p2 = np .arange (- (window_len // 2 ), 0 , dtype = int )
94+ results [window_half :] = window_p2
95+ return results * val
0 commit comments