Skip to content

Commit efc1c30

Browse files
roxyboyroxyboy
authored
Merge pull request #222 from miniufo/master
Change isotropize to satisfy Parseval's thereom. This PR addresses issue #219.
2 parents 9ac46f9 + aef1484 commit efc1c30

8 files changed

Lines changed: 1788 additions & 241 deletions

File tree

‎doc/DFT-iDFT_example.ipynb‎

Lines changed: 115 additions & 73 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@
2121
"import scipy.signal as signal\n",
2222
"import dask.array as dsar\n",
2323
"import matplotlib.pyplot as plt\n",
24+
"\n",
2425
"%matplotlib inline"
2526
]
2627
},
@@ -46,13 +47,13 @@
4647
"metadata": {},
4748
"outputs": [],
4849
"source": [
49-
"k0 = 1/0.52\n",
50-
"T = 4.\n",
50+
"k0 = 1 / 0.52\n",
51+
"T = 4.0\n",
5152
"dx = 0.02\n",
52-
"x = np.arange(-2*T,2*T,dx) \n",
53-
"y = np.cos(2*np.pi*k0*x) \n",
54-
"y[np.abs(x)>T/2]=0.\n",
55-
"da = xr.DataArray(y, dims=('x',), coords={'x':x})"
53+
"x = np.arange(-2 * T, 2 * T, dx)\n",
54+
"y = np.cos(2 * np.pi * k0 * x)\n",
55+
"y[np.abs(x) > T / 2] = 0.0\n",
56+
"da = xr.DataArray(y, dims=(\"x\",), coords={\"x\": x})"
5657
]
5758
},
5859
{
@@ -74,10 +75,10 @@
7475
}
7576
],
7677
"source": [
77-
"fig, ax = plt.subplots(figsize=(12,4)) \n",
78+
"fig, ax = plt.subplots(figsize=(12, 4))\n",
7879
"fig.set_tight_layout(True)\n",
79-
"da.plot(ax=ax, marker='+', label='original signal') \n",
80-
"ax.set_xlim([-8,8]);"
80+
"da.plot(ax=ax, marker=\"+\", label=\"original signal\")\n",
81+
"ax.set_xlim([-8, 8]);"
8182
]
8283
},
8384
{
@@ -95,8 +96,10 @@
9596
"metadata": {},
9697
"outputs": [],
9798
"source": [
98-
"da_dft = xrft.dft(da, true_phase=True, true_amplitude=True) # Fourier Transform w/ consideration of phase\n",
99-
"da_fft = xrft.fft(da) # Fourier Transform w/ numpy.fft-like behavior\n",
99+
"da_dft = xrft.dft(\n",
100+
" da, true_phase=True, true_amplitude=True\n",
101+
") # Fourier Transform w/ consideration of phase\n",
102+
"da_fft = xrft.fft(da) # Fourier Transform w/ numpy.fft-like behavior\n",
100103
"da_npft = npft.fft(da)"
101104
]
102105
},
@@ -106,8 +109,10 @@
106109
"metadata": {},
107110
"outputs": [],
108111
"source": [
109-
"k = da_dft.freq_x # wavenumber axis\n",
110-
"TF_s = T/2*(np.sinc(T*(k-k0)) + np.sinc(T*(k+k0))) # Theoretical result of the Fourier transform"
112+
"k = da_dft.freq_x # wavenumber axis\n",
113+
"TF_s = (\n",
114+
" T / 2 * (np.sinc(T * (k - k0)) + np.sinc(T * (k + k0)))\n",
115+
") # Theoretical result of the Fourier transform"
111116
]
112117
},
113118
{
@@ -129,26 +134,34 @@
129134
}
130135
],
131136
"source": [
132-
"fig, (ax1,ax2) = plt.subplots(figsize=(12,8), nrows=2, ncols=1)\n",
137+
"fig, (ax1, ax2) = plt.subplots(figsize=(12, 8), nrows=2, ncols=1)\n",
133138
"fig.set_tight_layout(True)\n",
134139
"\n",
135-
"(da_dft.real).plot(ax=ax1, linestyle='-', lw=3, c='k', label='phase preservation') \n",
136-
"((da_fft*dx).real).plot(ax=ax1, linestyle='', marker='+',label='no phase preservation') \n",
137-
"ax1.plot(k, (npft.fftshift(da_npft)*dx).real, linestyle='', marker='x',label='numpy fft')\n",
138-
"ax1.plot(k, TF_s.real, linestyle='--', label='Theory')\n",
139-
"ax1.set_xlim([-10,10])\n",
140-
"ax1.set_ylim([-2,2])\n",
140+
"(da_dft.real).plot(ax=ax1, linestyle=\"-\", lw=3, c=\"k\", label=\"phase preservation\")\n",
141+
"((da_fft * dx).real).plot(\n",
142+
" ax=ax1, linestyle=\"\", marker=\"+\", label=\"no phase preservation\"\n",
143+
")\n",
144+
"ax1.plot(\n",
145+
" k, (npft.fftshift(da_npft) * dx).real, linestyle=\"\", marker=\"x\", label=\"numpy fft\"\n",
146+
")\n",
147+
"ax1.plot(k, TF_s.real, linestyle=\"--\", label=\"Theory\")\n",
148+
"ax1.set_xlim([-10, 10])\n",
149+
"ax1.set_ylim([-2, 2])\n",
141150
"ax1.legend()\n",
142-
"ax1.set_title('REAL PART')\n",
151+
"ax1.set_title(\"REAL PART\")\n",
143152
"\n",
144-
"(da_dft.imag).plot(ax=ax2, linestyle='-', lw=3, c='k', label='phase preservation') \n",
145-
"((da_fft*dx).imag).plot(ax=ax2, linestyle='', marker='+', label='no phase preservation') \n",
146-
"ax2.plot(k, (npft.fftshift(da_npft)*dx).imag, linestyle='', marker='x',label='numpy fft')\n",
147-
"ax2.plot(k, TF_s.imag, linestyle='--', label='Theory')\n",
148-
"ax2.set_xlim([-10,10])\n",
149-
"ax2.set_ylim([-2,2])\n",
153+
"(da_dft.imag).plot(ax=ax2, linestyle=\"-\", lw=3, c=\"k\", label=\"phase preservation\")\n",
154+
"((da_fft * dx).imag).plot(\n",
155+
" ax=ax2, linestyle=\"\", marker=\"+\", label=\"no phase preservation\"\n",
156+
")\n",
157+
"ax2.plot(\n",
158+
" k, (npft.fftshift(da_npft) * dx).imag, linestyle=\"\", marker=\"x\", label=\"numpy fft\"\n",
159+
")\n",
160+
"ax2.plot(k, TF_s.imag, linestyle=\"--\", label=\"Theory\")\n",
161+
"ax2.set_xlim([-10, 10])\n",
162+
"ax2.set_ylim([-2, 2])\n",
150163
"ax2.legend()\n",
151-
"ax2.set_title('IMAGINARY PART');"
164+
"ax2.set_title(\"IMAGINARY PART\");"
152165
]
153166
},
154167
{
@@ -169,7 +182,9 @@
169182
"metadata": {},
170183
"outputs": [],
171184
"source": [
172-
"ida_dft = xrft.idft(da_dft, true_phase=True, true_amplitude=True) # Signal in direct space \n",
185+
"ida_dft = xrft.idft(\n",
186+
" da_dft, true_phase=True, true_amplitude=True\n",
187+
") # Signal in direct space\n",
173188
"ida_fft = xrft.ifft(da_fft)"
174189
]
175190
},
@@ -192,14 +207,16 @@
192207
}
193208
],
194209
"source": [
195-
"fig, ax = plt.subplots(figsize=(12,4))\n",
210+
"fig, ax = plt.subplots(figsize=(12, 4))\n",
196211
"fig.set_tight_layout(True)\n",
197-
"ida_dft.real.plot(ax=ax, linestyle='-', c='k', lw=4, label='phase preservation')\n",
198-
"ax.plot(x, ida_fft.real, linestyle='', marker='+', label='no phase preservation', alpha=.6) # w/out the phase information, the coordinates are lost\n",
199-
"da.plot(ax=ax, ls='--', lw=3, label='original signal')\n",
200-
"ax.plot(x, npft.ifft(da_npft).real, ls=':', label='inverse of numpy fft')\n",
201-
"ax.set_xlim([-8,8])\n",
202-
"ax.legend(loc='upper left');"
212+
"ida_dft.real.plot(ax=ax, linestyle=\"-\", c=\"k\", lw=4, label=\"phase preservation\")\n",
213+
"ax.plot(\n",
214+
" x, ida_fft.real, linestyle=\"\", marker=\"+\", label=\"no phase preservation\", alpha=0.6\n",
215+
") # w/out the phase information, the coordinates are lost\n",
216+
"da.plot(ax=ax, ls=\"--\", lw=3, label=\"original signal\")\n",
217+
"ax.plot(x, npft.ifft(da_npft).real, ls=\":\", label=\"inverse of numpy fft\")\n",
218+
"ax.set_xlim([-8, 8])\n",
219+
"ax.legend(loc=\"upper left\");"
203220
]
204221
},
205222
{
@@ -235,9 +252,9 @@
235252
"metadata": {},
236253
"outputs": [],
237254
"source": [
238-
"nshift = 70 # defining a shift\n",
239-
"x0 = dx*nshift \n",
240-
"nda = da.shift(x=nshift).dropna('x')"
255+
"nshift = 70 # defining a shift\n",
256+
"x0 = dx * nshift\n",
257+
"nda = da.shift(x=nshift).dropna(\"x\")"
241258
]
242259
},
243260
{
@@ -259,11 +276,11 @@
259276
}
260277
],
261278
"source": [
262-
"fig, ax = plt.subplots(figsize=(12,4))\n",
279+
"fig, ax = plt.subplots(figsize=(12, 4))\n",
263280
"fig.set_tight_layout(True)\n",
264-
"da.plot(ax=ax, label='original (centered) signal') \n",
265-
"nda.plot(ax=ax, marker='+', label='shifted signal', alpha=.6) \n",
266-
"ax.set_xlim([-8,nda.x.max()])\n",
281+
"da.plot(ax=ax, label=\"original (centered) signal\")\n",
282+
"nda.plot(ax=ax, marker=\"+\", label=\"shifted signal\", alpha=0.6)\n",
283+
"ax.set_xlim([-8, nda.x.max()])\n",
267284
"ax.legend();"
268285
]
269286
},
@@ -280,8 +297,10 @@
280297
"metadata": {},
281298
"outputs": [],
282299
"source": [
283-
"nda_dft = xrft.dft(nda, true_phase=True, true_amplitude=True) # Fourier Transform w/ phase preservation \n",
284-
"nda_fft = xrft.fft(nda) # Fourier Transform w/out phase preservation\n",
300+
"nda_dft = xrft.dft(\n",
301+
" nda, true_phase=True, true_amplitude=True\n",
302+
") # Fourier Transform w/ phase preservation\n",
303+
"nda_fft = xrft.fft(nda) # Fourier Transform w/out phase preservation\n",
285304
"nda_npft = npft.fft(nda)"
286305
]
287306
},
@@ -291,8 +310,13 @@
291310
"metadata": {},
292311
"outputs": [],
293312
"source": [
294-
"nk = nda_dft.freq_x # wavenumber axis\n",
295-
"TF_ns = T/2*(np.sinc(T*(nk-k0)) + np.sinc(T*(nk+k0)))*np.exp(-2j*np.pi*nk*x0) # Theoretical FT (Note the additional phase)"
313+
"nk = nda_dft.freq_x # wavenumber axis\n",
314+
"TF_ns = (\n",
315+
" T\n",
316+
" / 2\n",
317+
" * (np.sinc(T * (nk - k0)) + np.sinc(T * (nk + k0)))\n",
318+
" * np.exp(-2j * np.pi * nk * x0)\n",
319+
") # Theoretical FT (Note the additional phase)"
296320
]
297321
},
298322
{
@@ -314,26 +338,34 @@
314338
}
315339
],
316340
"source": [
317-
"fig, (ax1,ax2) = plt.subplots(figsize=(12,8), nrows=2, ncols=1)\n",
341+
"fig, (ax1, ax2) = plt.subplots(figsize=(12, 8), nrows=2, ncols=1)\n",
318342
"fig.set_tight_layout(True)\n",
319343
"\n",
320-
"(nda_dft.real).plot(ax=ax1, linestyle='-', lw=3, c='k', label='phase preservation') \n",
321-
"((nda_fft*dx).real).plot(ax=ax1, linestyle='', marker='+',label='no phase preservation')\n",
322-
"ax1.plot(nk, (npft.fftshift(nda_npft)*dx).real, linestyle='', marker='x',label='numpy fft')\n",
323-
"ax1.plot(nk, TF_ns.real, linestyle='--', label='Theory')\n",
324-
"ax1.set_xlim([-10,10])\n",
325-
"ax1.set_ylim([-2.,2])\n",
344+
"(nda_dft.real).plot(ax=ax1, linestyle=\"-\", lw=3, c=\"k\", label=\"phase preservation\")\n",
345+
"((nda_fft * dx).real).plot(\n",
346+
" ax=ax1, linestyle=\"\", marker=\"+\", label=\"no phase preservation\"\n",
347+
")\n",
348+
"ax1.plot(\n",
349+
" nk, (npft.fftshift(nda_npft) * dx).real, linestyle=\"\", marker=\"x\", label=\"numpy fft\"\n",
350+
")\n",
351+
"ax1.plot(nk, TF_ns.real, linestyle=\"--\", label=\"Theory\")\n",
352+
"ax1.set_xlim([-10, 10])\n",
353+
"ax1.set_ylim([-2.0, 2])\n",
326354
"ax1.legend()\n",
327-
"ax1.set_title('REAL PART')\n",
355+
"ax1.set_title(\"REAL PART\")\n",
328356
"\n",
329-
"(nda_dft.imag).plot(ax=ax2, linestyle='-', lw=3, c='k', label='phase preservation') \n",
330-
"((nda_fft*dx).imag).plot(ax=ax2, linestyle='', marker='+', label='no phase preservation') \n",
331-
"ax2.plot(nk, (npft.fftshift(nda_npft)*dx).imag, linestyle='', marker='x',label='numpy fft')\n",
332-
"ax2.plot(nk, TF_ns.imag, linestyle='--', label='Theory')\n",
333-
"ax2.set_xlim([-10,10])\n",
334-
"ax2.set_ylim([-2.,2.])\n",
357+
"(nda_dft.imag).plot(ax=ax2, linestyle=\"-\", lw=3, c=\"k\", label=\"phase preservation\")\n",
358+
"((nda_fft * dx).imag).plot(\n",
359+
" ax=ax2, linestyle=\"\", marker=\"+\", label=\"no phase preservation\"\n",
360+
")\n",
361+
"ax2.plot(\n",
362+
" nk, (npft.fftshift(nda_npft) * dx).imag, linestyle=\"\", marker=\"x\", label=\"numpy fft\"\n",
363+
")\n",
364+
"ax2.plot(nk, TF_ns.imag, linestyle=\"--\", label=\"Theory\")\n",
365+
"ax2.set_xlim([-10, 10])\n",
366+
"ax2.set_ylim([-2.0, 2.0])\n",
335367
"ax2.legend()\n",
336-
"ax2.set_title('IMAGINARY PART');"
368+
"ax2.set_title(\"IMAGINARY PART\");"
337369
]
338370
},
339371
{
@@ -356,7 +388,9 @@
356388
"metadata": {},
357389
"outputs": [],
358390
"source": [
359-
"inda_dft = xrft.idft(nda_dft, true_phase=True, true_amplitude=True) # Signal in direct space \n",
391+
"inda_dft = xrft.idft(\n",
392+
" nda_dft, true_phase=True, true_amplitude=True\n",
393+
") # Signal in direct space\n",
360394
"inda_fft = xrft.ifft(nda_fft)"
361395
]
362396
},
@@ -379,16 +413,24 @@
379413
}
380414
],
381415
"source": [
382-
"fig, ax = plt.subplots(figsize=(12,4))\n",
416+
"fig, ax = plt.subplots(figsize=(12, 4))\n",
383417
"fig.set_tight_layout(True)\n",
384-
"inda_dft.real.plot(ax=ax, linestyle='-', c='k', lw=4, label='phase preservation')\n",
385-
"ax.plot(x[:len(inda_fft.real)], inda_fft.real, linestyle='', marker='o', alpha=.7, \n",
386-
" label='no phase preservation (w/out shifting)')\n",
387-
"ax.plot(x[nshift:], inda_fft.real, linestyle='', marker='+', label='no phase preservation')\n",
388-
"nda.plot(ax=ax, ls='--', lw=3, label='original signal')\n",
389-
"ax.plot(x[nshift:], npft.ifft(nda_npft).real, ls=':', label='inverse of numpy fft')\n",
390-
"ax.set_xlim([nda.x.min(),nda.x.max()])\n",
391-
"ax.legend(loc='upper left');"
418+
"inda_dft.real.plot(ax=ax, linestyle=\"-\", c=\"k\", lw=4, label=\"phase preservation\")\n",
419+
"ax.plot(\n",
420+
" x[: len(inda_fft.real)],\n",
421+
" inda_fft.real,\n",
422+
" linestyle=\"\",\n",
423+
" marker=\"o\",\n",
424+
" alpha=0.7,\n",
425+
" label=\"no phase preservation (w/out shifting)\",\n",
426+
")\n",
427+
"ax.plot(\n",
428+
" x[nshift:], inda_fft.real, linestyle=\"\", marker=\"+\", label=\"no phase preservation\"\n",
429+
")\n",
430+
"nda.plot(ax=ax, ls=\"--\", lw=3, label=\"original signal\")\n",
431+
"ax.plot(x[nshift:], npft.ifft(nda_npft).real, ls=\":\", label=\"inverse of numpy fft\")\n",
432+
"ax.set_xlim([nda.x.min(), nda.x.max()])\n",
433+
"ax.legend(loc=\"upper left\");"
392434
]
393435
},
394436
{
@@ -1645,7 +1687,7 @@
16451687
}
16461688
],
16471689
"source": [
1648-
"fig, (ax1,ax2) = plt.subplots(figsize=(12,4), nrows=1, ncols=2)\n",
1690+
"fig, (ax1, ax2) = plt.subplots(figsize=(12, 4), nrows=1, ncols=2)\n",
16491691
"da.isel(time=0).plot(ax=ax1)\n",
16501692
"Fda_1.real.plot(ax=ax2)"
16511693
]

0 commit comments

Comments
 (0)