by Ian Czekala (iczekala@cfa.harvard.edu)
Questions to address:
A bunch of information here: http://docs.scipy.org/doc/numpy/reference/routines.fft.html
In [33]:
import numpy as np
from numpy.fft import fft,ifft,fftshift,ifftshift,fftfreq
import matplotlib.pyplot as plt
from IPython.display import Image
from scipy.signal import hann,kaiser,boxcar
%matplotlib inline
The numpy.fft convention for the fft and ifft is as follows, respectively
In [34]:
Image(url="http://docs.scipy.org/doc/numpy/_images/math/dd1819fee3d9660e590d152f0c7cb7aa0c441ff1.png")
Out[34]:
In [35]:
Image(url="http://docs.scipy.org/doc/numpy/_images/math/20aa8eba345ae91d3fff514f719d1424bccf0850.png")
Out[35]:
Bracewell's convention has the 1/n term on the fft, instead, and no 1/n on the ifft, but otherwise the conventions are the same.
In [36]:
x_odd = np.arange(-5,5.1,1.)
In [37]:
len(x_odd)
x_odd
Out[37]:
In [38]:
x_even = np.arange(-5,5,1.)
In [39]:
len(x_even)
x_even
Out[39]:
In [40]:
fftshift(x_odd)
Out[40]:
In [41]:
ifftshift(x_odd)
Out[41]:
In [42]:
fftshift(fftshift(x_odd))
Out[42]:
In [43]:
fftshift(x_even)
Out[43]:
In [44]:
ifftshift(x_even)
Out[44]:
In [45]:
fftshift(fftshift(x_even))
Out[45]:
Conclusion: for an even length array, fftshift and ifftshift are identical operations, and simply swap halves of the array. The difference comes in when you have an odd-length array. From the numpy.fft documentation, fftshift is designed to "Shift the zero-frequency component to the center of the spectrum." Meaning, that this takes the output of the fft routine which always has the zeroth element as the zero-frequency component and shifts it from the beginning to the middle of the array. Therefore, ifftshift must be designed to do the opposite of this, and is what you want to do when you have an array that stretches from -x_lim to +x_lim and you wish to perform an fft on it.
In [46]:
xs_odd = np.linspace(-2.,2.,num=81)
xs_odd
Out[46]:
In [47]:
@np.vectorize
def boxcar(x):
if np.abs(x) < 0.5:
return 1.
if np.abs(x) == 0.5:
return 0.5
else:
return 0.
In [48]:
ys = boxcar(xs_odd)
In [49]:
ys
Out[49]:
Just to see what this looks like
In [50]:
plt.clf()
plt.plot(xs_odd,ys,drawstyle='steps-mid')
Out[50]:
In [51]:
out_noshift = fft(ys)
#out_noshift
In [52]:
plt.clf()
plt.plot(np.real(out_noshift))
plt.plot(np.imag(out_noshift))
Out[52]:
Following the lessons from the previous section, lets fftshift back, since we are in the frequency domain and we wish to put the zeroth-frequency component in the center!
In [53]:
shift_out_noshift = fftshift(out_noshift)
plt.clf()
plt.plot(np.real(shift_out_noshift))
plt.plot(np.imag(shift_out_noshift))
Out[53]:
Here we see that there is a serious imaginary component to our spectrum, something that we should not have if we indeed took the Fourier transfrom of a boxcar, which is a real and even function, therefore its Fourier transform should also be real and even. If we plot the absolute magnitude of this, however, we should be able to recover what we want
In [54]:
plt.clf()
plt.plot(np.abs(shift_out_noshift))
Out[54]:
However, we see that we don't have any negative values, we've just recovered np.abs(sinc)
In [55]:
out_shift = fft(ifftshift(ys))
#out_shift
In [56]:
plt.clf()
plt.plot(np.real(fftshift(out_shift)))
plt.plot(np.imag(fftshift(out_shift)))
Out[56]:
Here we recover exactly what we wanted, a sinc function. We see that the output produces no phase, this is because we packed the fft correctly with the zeroth term in the 0 index, and the routine assumed the function was periodic. What if we tried this with ifft?
In [57]:
iout_noshift = ifft(ys)
#iout_noshift
In [58]:
plt.clf()
plt.plot(np.real(iout_noshift))
plt.plot(np.imag(iout_noshift))
Out[58]:
In [59]:
shift_iout_noshift = fftshift(iout_noshift)
plt.clf()
plt.plot(np.real(shift_iout_noshift))
plt.plot(np.imag(shift_iout_noshift))
Out[59]:
It seems as though we have recovered the same thing as before, but now the amplitude is scaled by 1/n
In [60]:
iout_shift = ifft(ifftshift(ys))
#iout_shift
In [60]:
In [61]:
shift_iout_shift = fftshift(iout_shift)
plt.clf()
plt.plot(np.real(shift_iout_shift))
plt.plot(np.imag(shift_iout_shift))
Out[61]:
And we see that we recover our sinc function
We should be able to recover the same function by doing ifft on the fft, even if we didn't pre-process correctly at the beginning by doing ifftshift, as long as we preserve the imaginary components, because really by not doing ifftshift all we did was introduce a phase offset. But have we introduced high-frequency components that shouldn't be there? Let's try with a Gaussian that is windowed and not windowed.
In [62]:
xs_odd
Out[62]:
In [63]:
@np.vectorize
def gauss(x):
return np.exp(-0.5 * x**2)
In [64]:
ygs = np.exp(-0.5 * xs_odd**2)
wygs = hann(len(ygs))*ygs
In [65]:
plt.clf()
plt.plot(ygs)
plt.plot(wygs)
plt.show()
In [66]:
gout = fftshift(fft(ifftshift(ygs)))
plt.clf()
plt.plot(np.real(gout))
plt.plot(np.imag(gout))
plt.show()
There is no imaginary component, however there is a bit of ringing thanks to the fact that the series didn't go smoothly to zero at the edges. If we show the one with a window, then there should be no ringing. I think the padding is only important if you really care about the frequency spectrum. If you're going to do something like a Gaussian taper/cutoff in the frequency domain, then transform back to the real domain, you shouldn't need to window.
In [67]:
wgout = fft(ifftshift(wygs))
wgout_shift = fftshift(wgout)
plt.clf()
plt.plot(np.real(wgout_shift))
plt.plot(np.imag(wgout_shift))
plt.show()
In [68]:
wgout_raw = fft(wygs)
plt.clf()
plt.plot(np.real(wgout_raw))
plt.plot(np.imag(wgout_raw))
plt.show()
In [69]:
wback = ifft(wgout_raw)
plt.clf()
plt.plot(np.real(wback))
plt.plot(np.imag(wback))
plt.show()
Conclusion: If your end goal is to have what you started with back at the end of the day, there is no need to do ifftshift/fftshift, everything will be ok. This is important if you want the Fourier domain representation to be nice.
I think this should work that you always want to pad nearest the high frequencies. This means that you basically insert zeros in the middle.
In [70]:
xs_odd
Out[70]:
In [71]:
N = len(xs_odd)
print(N)
freqs = fftfreq(N,d=0.05)
print(freqs)
freqs[40], freqs[41]
print(1./(0.05 * 2))
print(80./81 * 10.0)
Can we recover xs_odd from fftfreqs? The answer is yes!
In [72]:
D = freqs[1]
print(D)
fftfreq(N,d=D)
Out[72]:
In [73]:
wpad = np.concatenate((wgout[:41],np.zeros((80,)),wgout[41:])) #want to keep the FFT odd
rest = fftshift(ifft(wpad))
plt.clf()
plt.plot(np.real(rest))
plt.plot(np.imag(rest))
plt.show()
In [74]:
print(rest[0],rest[-1])
print(wygs[0],wygs[-1])
Inserting zeros in the middle and transforming back seems to have worked correctly.
In [75]:
Npad = len(wpad)
print(Npad)
xs_pad_freq = fftshift(fftfreq(Npad,d=D))
print(xs_pad_freq)
In [76]:
xs_odd_pad = np.linspace(-2,2.,num=161)
xs_high = np.linspace(-2,2,num=2000)
wy_high = hann(2000) * gauss(xs_high)
In [77]:
plt.clf()
plt.plot(xs_high, wy_high)
#plt.plot(xs_odd, wygs,"bo")
plt.plot(xs_pad_freq, 161/81. * rest,"o")
plt.xlim(-1.5,-1.)
plt.show()
The returned float array contains the frequency bins in cycles/unit (with zero at the start) given a window length n and a sample spacing d:
f = [0, 1, ..., n/2-1, -n/2, ..., -1] / (d*n) if n is even
f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (d*n) if n is odd
How should you change the fftfreqs if you've padded in the Fourier domain and wish to get the correct real spacings back? Turns out that the correct thing is to use fftfreq with the original spacing of the unpadded array but the length of the padded array.
In [78]:
wave_grid = np.load('wave_grid_trunc.npy')
f_grid = np.load('f_grid.npy')
len(wave_grid)
Out[78]:
In [79]:
plt.clf()
plt.plot(wave_grid,f_grid)
plt.xlabel(r"$\lambda$")
Out[79]:
In [80]:
@np.vectorize
def gauss_taper(s,sigma=2.89):
'''This is the FT of a gaussian w/ this sigma. Sigma in km/s'''
return np.exp(-2 * np.pi**2 * sigma*2 * s**2)
def convolve_gauss(wl,fl,sigma=2.89,spacing=0.35):
##Take FFT of f_grid
out = fft(fl)
N = len(fl)
freqs = fftfreq(N,d=spacing)
taper = gauss_taper(freqs,sigma)
tout = out * taper
blended = np.fft.ifft(tout)
return blended #remove tiny complex component
In [81]:
f_blend = convolve_gauss(wave_grid,f_grid)
In [82]:
plt.clf()
plt.plot(wave_grid,np.real(f_blend))
plt.xlabel(r"$\lambda$")
Out[82]:
In [83]:
plt.clf()
plt.plot(wave_grid,np.imag(f_blend),'g')
plt.xlabel(r"$\lambda$")
Out[83]:
In [83]: