Chris@87
|
1 """
|
Chris@87
|
2 Discrete Fourier Transforms - helper.py
|
Chris@87
|
3
|
Chris@87
|
4 """
|
Chris@87
|
5 from __future__ import division, absolute_import, print_function
|
Chris@87
|
6
|
Chris@87
|
7 from numpy.compat import integer_types
|
Chris@87
|
8 from numpy.core import (
|
Chris@87
|
9 asarray, concatenate, arange, take, integer, empty
|
Chris@87
|
10 )
|
Chris@87
|
11
|
Chris@87
|
12 # Created by Pearu Peterson, September 2002
|
Chris@87
|
13
|
Chris@87
|
14 __all__ = ['fftshift', 'ifftshift', 'fftfreq', 'rfftfreq']
|
Chris@87
|
15
|
Chris@87
|
16 integer_types = integer_types + (integer,)
|
Chris@87
|
17
|
Chris@87
|
18
|
Chris@87
|
19 def fftshift(x, axes=None):
|
Chris@87
|
20 """
|
Chris@87
|
21 Shift the zero-frequency component to the center of the spectrum.
|
Chris@87
|
22
|
Chris@87
|
23 This function swaps half-spaces for all axes listed (defaults to all).
|
Chris@87
|
24 Note that ``y[0]`` is the Nyquist component only if ``len(x)`` is even.
|
Chris@87
|
25
|
Chris@87
|
26 Parameters
|
Chris@87
|
27 ----------
|
Chris@87
|
28 x : array_like
|
Chris@87
|
29 Input array.
|
Chris@87
|
30 axes : int or shape tuple, optional
|
Chris@87
|
31 Axes over which to shift. Default is None, which shifts all axes.
|
Chris@87
|
32
|
Chris@87
|
33 Returns
|
Chris@87
|
34 -------
|
Chris@87
|
35 y : ndarray
|
Chris@87
|
36 The shifted array.
|
Chris@87
|
37
|
Chris@87
|
38 See Also
|
Chris@87
|
39 --------
|
Chris@87
|
40 ifftshift : The inverse of `fftshift`.
|
Chris@87
|
41
|
Chris@87
|
42 Examples
|
Chris@87
|
43 --------
|
Chris@87
|
44 >>> freqs = np.fft.fftfreq(10, 0.1)
|
Chris@87
|
45 >>> freqs
|
Chris@87
|
46 array([ 0., 1., 2., 3., 4., -5., -4., -3., -2., -1.])
|
Chris@87
|
47 >>> np.fft.fftshift(freqs)
|
Chris@87
|
48 array([-5., -4., -3., -2., -1., 0., 1., 2., 3., 4.])
|
Chris@87
|
49
|
Chris@87
|
50 Shift the zero-frequency component only along the second axis:
|
Chris@87
|
51
|
Chris@87
|
52 >>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
|
Chris@87
|
53 >>> freqs
|
Chris@87
|
54 array([[ 0., 1., 2.],
|
Chris@87
|
55 [ 3., 4., -4.],
|
Chris@87
|
56 [-3., -2., -1.]])
|
Chris@87
|
57 >>> np.fft.fftshift(freqs, axes=(1,))
|
Chris@87
|
58 array([[ 2., 0., 1.],
|
Chris@87
|
59 [-4., 3., 4.],
|
Chris@87
|
60 [-1., -3., -2.]])
|
Chris@87
|
61
|
Chris@87
|
62 """
|
Chris@87
|
63 tmp = asarray(x)
|
Chris@87
|
64 ndim = len(tmp.shape)
|
Chris@87
|
65 if axes is None:
|
Chris@87
|
66 axes = list(range(ndim))
|
Chris@87
|
67 elif isinstance(axes, integer_types):
|
Chris@87
|
68 axes = (axes,)
|
Chris@87
|
69 y = tmp
|
Chris@87
|
70 for k in axes:
|
Chris@87
|
71 n = tmp.shape[k]
|
Chris@87
|
72 p2 = (n+1)//2
|
Chris@87
|
73 mylist = concatenate((arange(p2, n), arange(p2)))
|
Chris@87
|
74 y = take(y, mylist, k)
|
Chris@87
|
75 return y
|
Chris@87
|
76
|
Chris@87
|
77
|
Chris@87
|
78 def ifftshift(x, axes=None):
|
Chris@87
|
79 """
|
Chris@87
|
80 The inverse of `fftshift`. Although identical for even-length `x`, the
|
Chris@87
|
81 functions differ by one sample for odd-length `x`.
|
Chris@87
|
82
|
Chris@87
|
83 Parameters
|
Chris@87
|
84 ----------
|
Chris@87
|
85 x : array_like
|
Chris@87
|
86 Input array.
|
Chris@87
|
87 axes : int or shape tuple, optional
|
Chris@87
|
88 Axes over which to calculate. Defaults to None, which shifts all axes.
|
Chris@87
|
89
|
Chris@87
|
90 Returns
|
Chris@87
|
91 -------
|
Chris@87
|
92 y : ndarray
|
Chris@87
|
93 The shifted array.
|
Chris@87
|
94
|
Chris@87
|
95 See Also
|
Chris@87
|
96 --------
|
Chris@87
|
97 fftshift : Shift zero-frequency component to the center of the spectrum.
|
Chris@87
|
98
|
Chris@87
|
99 Examples
|
Chris@87
|
100 --------
|
Chris@87
|
101 >>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
|
Chris@87
|
102 >>> freqs
|
Chris@87
|
103 array([[ 0., 1., 2.],
|
Chris@87
|
104 [ 3., 4., -4.],
|
Chris@87
|
105 [-3., -2., -1.]])
|
Chris@87
|
106 >>> np.fft.ifftshift(np.fft.fftshift(freqs))
|
Chris@87
|
107 array([[ 0., 1., 2.],
|
Chris@87
|
108 [ 3., 4., -4.],
|
Chris@87
|
109 [-3., -2., -1.]])
|
Chris@87
|
110
|
Chris@87
|
111 """
|
Chris@87
|
112 tmp = asarray(x)
|
Chris@87
|
113 ndim = len(tmp.shape)
|
Chris@87
|
114 if axes is None:
|
Chris@87
|
115 axes = list(range(ndim))
|
Chris@87
|
116 elif isinstance(axes, integer_types):
|
Chris@87
|
117 axes = (axes,)
|
Chris@87
|
118 y = tmp
|
Chris@87
|
119 for k in axes:
|
Chris@87
|
120 n = tmp.shape[k]
|
Chris@87
|
121 p2 = n-(n+1)//2
|
Chris@87
|
122 mylist = concatenate((arange(p2, n), arange(p2)))
|
Chris@87
|
123 y = take(y, mylist, k)
|
Chris@87
|
124 return y
|
Chris@87
|
125
|
Chris@87
|
126
|
Chris@87
|
127 def fftfreq(n, d=1.0):
|
Chris@87
|
128 """
|
Chris@87
|
129 Return the Discrete Fourier Transform sample frequencies.
|
Chris@87
|
130
|
Chris@87
|
131 The returned float array `f` contains the frequency bin centers in cycles
|
Chris@87
|
132 per unit of the sample spacing (with zero at the start). For instance, if
|
Chris@87
|
133 the sample spacing is in seconds, then the frequency unit is cycles/second.
|
Chris@87
|
134
|
Chris@87
|
135 Given a window length `n` and a sample spacing `d`::
|
Chris@87
|
136
|
Chris@87
|
137 f = [0, 1, ..., n/2-1, -n/2, ..., -1] / (d*n) if n is even
|
Chris@87
|
138 f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (d*n) if n is odd
|
Chris@87
|
139
|
Chris@87
|
140 Parameters
|
Chris@87
|
141 ----------
|
Chris@87
|
142 n : int
|
Chris@87
|
143 Window length.
|
Chris@87
|
144 d : scalar, optional
|
Chris@87
|
145 Sample spacing (inverse of the sampling rate). Defaults to 1.
|
Chris@87
|
146
|
Chris@87
|
147 Returns
|
Chris@87
|
148 -------
|
Chris@87
|
149 f : ndarray
|
Chris@87
|
150 Array of length `n` containing the sample frequencies.
|
Chris@87
|
151
|
Chris@87
|
152 Examples
|
Chris@87
|
153 --------
|
Chris@87
|
154 >>> signal = np.array([-2, 8, 6, 4, 1, 0, 3, 5], dtype=float)
|
Chris@87
|
155 >>> fourier = np.fft.fft(signal)
|
Chris@87
|
156 >>> n = signal.size
|
Chris@87
|
157 >>> timestep = 0.1
|
Chris@87
|
158 >>> freq = np.fft.fftfreq(n, d=timestep)
|
Chris@87
|
159 >>> freq
|
Chris@87
|
160 array([ 0. , 1.25, 2.5 , 3.75, -5. , -3.75, -2.5 , -1.25])
|
Chris@87
|
161
|
Chris@87
|
162 """
|
Chris@87
|
163 if not isinstance(n, integer_types):
|
Chris@87
|
164 raise ValueError("n should be an integer")
|
Chris@87
|
165 val = 1.0 / (n * d)
|
Chris@87
|
166 results = empty(n, int)
|
Chris@87
|
167 N = (n-1)//2 + 1
|
Chris@87
|
168 p1 = arange(0, N, dtype=int)
|
Chris@87
|
169 results[:N] = p1
|
Chris@87
|
170 p2 = arange(-(n//2), 0, dtype=int)
|
Chris@87
|
171 results[N:] = p2
|
Chris@87
|
172 return results * val
|
Chris@87
|
173 #return hstack((arange(0,(n-1)/2 + 1), arange(-(n/2),0))) / (n*d)
|
Chris@87
|
174
|
Chris@87
|
175
|
Chris@87
|
176 def rfftfreq(n, d=1.0):
|
Chris@87
|
177 """
|
Chris@87
|
178 Return the Discrete Fourier Transform sample frequencies
|
Chris@87
|
179 (for usage with rfft, irfft).
|
Chris@87
|
180
|
Chris@87
|
181 The returned float array `f` contains the frequency bin centers in cycles
|
Chris@87
|
182 per unit of the sample spacing (with zero at the start). For instance, if
|
Chris@87
|
183 the sample spacing is in seconds, then the frequency unit is cycles/second.
|
Chris@87
|
184
|
Chris@87
|
185 Given a window length `n` and a sample spacing `d`::
|
Chris@87
|
186
|
Chris@87
|
187 f = [0, 1, ..., n/2-1, n/2] / (d*n) if n is even
|
Chris@87
|
188 f = [0, 1, ..., (n-1)/2-1, (n-1)/2] / (d*n) if n is odd
|
Chris@87
|
189
|
Chris@87
|
190 Unlike `fftfreq` (but like `scipy.fftpack.rfftfreq`)
|
Chris@87
|
191 the Nyquist frequency component is considered to be positive.
|
Chris@87
|
192
|
Chris@87
|
193 Parameters
|
Chris@87
|
194 ----------
|
Chris@87
|
195 n : int
|
Chris@87
|
196 Window length.
|
Chris@87
|
197 d : scalar, optional
|
Chris@87
|
198 Sample spacing (inverse of the sampling rate). Defaults to 1.
|
Chris@87
|
199
|
Chris@87
|
200 Returns
|
Chris@87
|
201 -------
|
Chris@87
|
202 f : ndarray
|
Chris@87
|
203 Array of length ``n//2 + 1`` containing the sample frequencies.
|
Chris@87
|
204
|
Chris@87
|
205 Examples
|
Chris@87
|
206 --------
|
Chris@87
|
207 >>> signal = np.array([-2, 8, 6, 4, 1, 0, 3, 5, -3, 4], dtype=float)
|
Chris@87
|
208 >>> fourier = np.fft.rfft(signal)
|
Chris@87
|
209 >>> n = signal.size
|
Chris@87
|
210 >>> sample_rate = 100
|
Chris@87
|
211 >>> freq = np.fft.fftfreq(n, d=1./sample_rate)
|
Chris@87
|
212 >>> freq
|
Chris@87
|
213 array([ 0., 10., 20., 30., 40., -50., -40., -30., -20., -10.])
|
Chris@87
|
214 >>> freq = np.fft.rfftfreq(n, d=1./sample_rate)
|
Chris@87
|
215 >>> freq
|
Chris@87
|
216 array([ 0., 10., 20., 30., 40., 50.])
|
Chris@87
|
217
|
Chris@87
|
218 """
|
Chris@87
|
219 if not isinstance(n, integer_types):
|
Chris@87
|
220 raise ValueError("n should be an integer")
|
Chris@87
|
221 val = 1.0/(n*d)
|
Chris@87
|
222 N = n//2 + 1
|
Chris@87
|
223 results = arange(0, N, dtype=int)
|
Chris@87
|
224 return results * val
|