Playing around with FFTs in Python



In [1]:

    
print arange(5)









    



[0 1 2 3 4]



In [2]:

    
spectral_size = 50
size = 1*spectral_size



In [3]:

    
f = []
g = []
for i in arange(size):
    f.append(math.sin(2*pi*i/size))
    g.append(math.cos(2*pi*i/size))
    for j in arange(spectral_size):
        f[len(f)-1] += random.random()*math.sin(2*pi*j*i/size)
        g[len(g)-1] += random.random()*math.cos(2*pi*j*i/size)

We want to transform these functions to space, multiply them together and then transform them back



In [4]:

    
print f









    



[0.0, 0.5105310830714515, 0.8291321137694428, -0.4350561486893793, 1.3505277626315526, -0.4651721473368178, 2.4865390458216163, 1.1419317787015832, 1.8592328782027467, 1.0815648513198204, -1.6624365969652886, 2.095700329450023, 0.8160044696871677, -1.2209174844235904, -0.33412522990680316, 0.6401300385133023, -1.028171706476379, 3.780464244648382, 1.2844964678809627, -1.75307224156116, 2.332595400200187, -1.8912082246491622, 0.5057445252576036, 1.030290849058479, -0.631576971499558, -2.5539220183814744e-14, -0.14281389799408126, 3.0529748953602174, -0.0004680652924593104, 1.6200236558385024, -0.5633173741304043, -0.19006665172024007, -0.8963779315762418, -1.3444073905177973, -1.726909436035283, 2.1707726783616175, -0.9437721073837821, -1.109121461554114, -0.7982622644779716, -0.925229060388547, -0.8900923120283409, -2.365462034668658, -0.9746134280599484, -2.174914219103794, -0.25952128214209913, -0.7428297574092118, 0.040932356221215216, -1.1426751925422438, 1.4050789888948505, -0.631508190899547]



In [5]:

    
pyplot.plot(arange(size)*1.0/size, f)
pyplot.plot(arange(size)*1.0/size, g)
pyplot.show()



In [6]:

    
npf = np.array(f)



In [7]:

    
print npf









    



[  0.00000000e+00   5.10531083e-01   8.29132114e-01  -4.35056149e-01
   1.35052776e+00  -4.65172147e-01   2.48653905e+00   1.14193178e+00
   1.85923288e+00   1.08156485e+00  -1.66243660e+00   2.09570033e+00
   8.16004470e-01  -1.22091748e+00  -3.34125230e-01   6.40130039e-01
  -1.02817171e+00   3.78046424e+00   1.28449647e+00  -1.75307224e+00
   2.33259540e+00  -1.89120822e+00   5.05744525e-01   1.03029085e+00
  -6.31576971e-01  -2.55392202e-14  -1.42813898e-01   3.05297490e+00
  -4.68065292e-04   1.62002366e+00  -5.63317374e-01  -1.90066652e-01
  -8.96377932e-01  -1.34440739e+00  -1.72690944e+00   2.17077268e+00
  -9.43772107e-01  -1.10912146e+00  -7.98262264e-01  -9.25229060e-01
  -8.90092312e-01  -2.36546203e+00  -9.74613428e-01  -2.17491422e+00
  -2.59521282e-01  -7.42829757e-01   4.09323562e-02  -1.14267519e+00
   1.40507899e+00  -6.31508191e-01]



In [9]:

    
npf_fft = numpy.fft.fft(npf)



In [10]:

    
print npf_fft









    



[  2.79056960 +0.00000000e+00j  -3.15314402 -1.90795915e+01j
   8.14498421 -9.83138446e+00j   0.68583050 -7.19857641e-01j
  -1.67798642 -1.41768849e+00j  -2.64785603 +1.45942268e+01j
   0.41580653 -5.09887587e+00j   2.53175875 -6.09497942e-01j
  -1.95032747 +3.72968618e+00j   6.89747540 -7.73667438e+00j
  -4.84183809 -4.50463837e+00j  -2.57526731 +8.32556886e+00j
  -5.93418119 +6.89127542e-01j   4.91456041 -3.76536610e+00j
  -0.98837261 +2.38161133e+00j  -3.98006381 -4.75262719e+00j
 -12.22451581 -3.56666387e+00j  11.20722388 +8.14679072e+00j
   1.20938688 -1.25207257e+01j -11.58807319 -7.80546386e+00j
   5.49567811 +5.62556264e+00j  -9.91129879 -3.16365137e+00j
  -1.64778715 +1.05731776e+01j   6.95631360 -9.83701621e+00j
  12.60386819 +1.67131923e+00j   1.32508121 -1.86517468e-14j
  12.60386819 -1.67131923e+00j   6.95631360 +9.83701621e+00j
  -1.64778715 -1.05731776e+01j  -9.91129879 +3.16365137e+00j
   5.49567811 -5.62556264e+00j -11.58807319 +7.80546386e+00j
   1.20938688 +1.25207257e+01j  11.20722388 -8.14679072e+00j
 -12.22451581 +3.56666387e+00j  -3.98006381 +4.75262719e+00j
  -0.98837261 -2.38161133e+00j   4.91456041 +3.76536610e+00j
  -5.93418119 -6.89127542e-01j  -2.57526731 -8.32556886e+00j
  -4.84183809 +4.50463837e+00j   6.89747540 +7.73667438e+00j
  -1.95032747 -3.72968618e+00j   2.53175875 +6.09497942e-01j
   0.41580653 +5.09887587e+00j  -2.64785603 -1.45942268e+01j
  -1.67798642 +1.41768849e+00j   0.68583050 +7.19857641e-01j
   8.14498421 +9.83138446e+00j  -3.15314402 +1.90795915e+01j]



In [73]:

    
pyplot.plot(arange(size)*1.0/size, imag(npf_fft), 'b')
pyplot.plot(arange(size)*1.0/size, real(npf_fft), 'r')
pyplot.plot(arange(size)*1.0/size, abs(npf_fft), 'k')
pyplot.show()



In [75]:

    
npf_fft_ifft = numpy.fft.ifft(npf_fft)



In [77]:

    
pyplot.plot(arange(size)*1.0/size, npf, 'b')
pyplot.plot(arange(size)*1.0/size, real(npf_fft_ifft), 'r')
pyplot.show()



In [78]:

    
print npf-npf_fft_ifft









    



[  0.00000000e+00 -2.04281037e-16j  -2.22044605e-15 -1.32338585e-15j
  -4.99600361e-16 -4.12710509e-16j  -8.88178420e-16 +9.02549498e-16j
   1.38083989e-15 -2.22933760e-16j  -1.44328993e-15 +1.46879540e-15j
   7.77156117e-16 +2.15991895e-15j  -1.11022302e-15 +2.49201297e-15j
  -4.44089210e-16 +8.20516343e-16j  -4.44089210e-16 +2.81050368e-15j
   2.22044605e-15 +9.29320162e-16j   5.55111512e-16 -4.38646016e-17j
   8.88178420e-16 -1.34977040e-15j  -1.11022302e-15 +5.20228387e-16j
  -2.57432964e-15 +1.08234970e-16j   2.66453526e-15 -2.59016536e-15j
  -2.22044605e-16 +2.05734306e-16j   1.24900090e-15 -4.06898339e-16j
   1.83186799e-15 +2.77858616e-15j  -2.49800181e-16 +4.35002785e-15j
  -2.27595720e-15 -6.39160943e-15j  -1.22124533e-15 +2.69733671e-16j
  -1.99840144e-15 +3.48149960e-16j   4.44089210e-16 -1.53445015e-15j
   5.88418203e-15 -3.00444760e-15j  -1.98427968e-15 +9.27859751e-16j
  -4.44089210e-16 -2.61597452e-15j   2.10942375e-15 +5.78447400e-16j
  -1.32185929e-15 -1.31156385e-16j  -3.10862447e-15 +2.28237032e-15j
  -1.11022302e-15 -7.65018443e-16j   8.88178420e-16 +2.05922957e-16j
  -2.33146835e-15 +1.46583803e-16j  -1.47104551e-15 +1.79613837e-16j
   2.66453526e-15 +4.43651352e-15j   6.68909372e-15 +1.21322994e-15j
   1.38777878e-15 +3.94361350e-15j   8.88178420e-16 -5.55972817e-15j
   1.36002321e-15 +6.98043316e-16j  -2.44249065e-15 -2.67601043e-16j
  -2.66453526e-15 +3.98909809e-15j   1.70870262e-15 -7.51536257e-16j
   3.75394160e-15 -6.18169722e-16j   8.88178420e-16 -1.23953596e-15j
  -4.44089210e-15 -5.98474960e-15j   3.55271368e-15 -5.75496922e-15j
  -2.49800181e-15 -1.42021030e-15j  -3.10862447e-15 +1.14537626e-15j
  -1.05471187e-15 +3.66694165e-15j  -3.10862447e-15 -9.84759993e-16j]

Figuring out how aliasing works

Now we're going to try transforming two functions to space, taking their product, and then transforming them back to spatial coordinates. We'll do this with two functions that are band-limited (eg. only N Fourier components) - one case where their product should still be within the N Fourier components, and one case where some of the product should have frequencies above our cutoff. We'll try to reconstruct the Fourier components of the product from the analytical results obtained by hand.



In [1]:

    
spectral_n = 10
spectral_size = spectral_n + 1
spatial_size = 2*spectral_size + 1



In [106]:

    
a = zeros(spectral_size)
b = zeros(spectral_size)
print a









    



[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]



In [107]:

    
a[0] = 1.0
a[1] = 2.0
a[6] = 1.0
b[1] = 0.5
b[2] = 3.0
print a, b









    



[ 1.  2.  0.  0.  0.  0.  1.  0.  0.  0.  0.] [ 0.   0.5  3.   0.   0.   0.   0.   0.   0.   0.   0. ]



In [4]:

    
a_spatial = []
b_spatial = []
for i in arange(spatial_size):
    a_spatial.append(0.0)
    b_spatial.append(0.0)
    for n in arange(spectral_size):
        a_spatial[i] += a[n]*math.cos(2*pi*n*i/(spatial_size))
        b_spatial[i] += b[n]*math.sin(2*pi*n*i/(spatial_size))



In [5]:

    
pyplot.plot(arange(spatial_size), a_spatial, 'b')
pyplot.plot(arange(spatial_size), b_spatial, 'r')
pyplot.show()

There are several ways of FFT'ing the data here. We can do an FFT of size spectral_size on the individual arrays a & b (of course producing output array also of length spectral_size). The issue here is that although this FFT is completely reversible (we get the same frequencies out of the inverse transform), the spatial information can have higher frequencies hidden in it (aliased) that will not show up.



In [6]:

    
a_ifft = numpy.fft.ifft(np.array(a))
b_ifft = numpy.fft.ifft(np.array(b))



In [7]:

    
a_spatial_undersampled = []
b_spatial_undersampled = []
for i in arange(spectral_size):
    a_spatial_undersampled.append(0.0)
    b_spatial_undersampled.append(0.0)
    for n in arange(spectral_size):
        a_spatial_undersampled[i] += a[n]*math.cos(2*pi*n*i/(spectral_size))
        b_spatial_undersampled[i] += b[n]*math.sin(2*pi*n*i/(spectral_size))



In [9]:

    
pyplot.plot(arange(spectral_size)*1.0/(spectral_size), (spectral_size)*real(a_ifft), 'b')
pyplot.plot(arange(spectral_size)*1.0/(spectral_size), (spectral_size)*imag(a_ifft), 'r')
pyplot.plot(arange(spectral_size)*1.0/(spectral_size), a_spatial_undersampled, 'k')
pyplot.show()

Here we see that that real part of the iFFT (of our real array representing cosine terms) is the same as the Fourier sum above (sampled at spectral_size evenly spaced points rather than spatial_size evenly spaced points). Not sure how to interpret the imaginary points yet... Additionally, the iFFT routine in Python appears to have the entire normalization term, so we need to multiply by spectral_size to get the correct spatial terms. We can do a similar exercise for the sine terms:



In [9]:

    
pyplot.plot(arange(spectral_size)*1.0/(spectral_size), (spectral_size)*real(b_ifft), 'b')
pyplot.plot(arange(spectral_size)*1.0/(spectral_size), (spectral_size)*imag(b_ifft), 'r')
pyplot.plot(arange(spectral_size)*1.0/(spectral_size), b_spatial_undersampled, 'k')
pyplot.show()

In this case, it's the imaginary terms that represent the actual spatial function we're representing with Fourier terms. Again, not entirely sure about the real terms in this case. Perhaps it's an artifact of the FFT making the signal periodic over all space (it looks a bit like some kind of odd/even extention).

If we transform these 'undersampled' results back to frequency space with the standard FFT then we should get our original frequency arrays back:



In [6]:

    
a_ifft_fft = numpy.fft.fft(a_ifft)
b_ifft_fft = numpy.fft.fft(b_ifft)



In [7]:

    
pyplot.plot(arange(spectral_size), real(a_ifft_fft), 'b')
pyplot.plot(arange(spectral_size), imag(a_ifft_fft), 'r')
pyplot.plot(arange(spectral_size), a, 'k')
pyplot.show()



In [8]:

    
pyplot.plot(arange(spectral_size), real(b_ifft_fft), 'b')
pyplot.plot(arange(spectral_size), imag(b_ifft_fft), 'r')
pyplot.plot(arange(spectral_size), b, 'k')
pyplot.show()

In both cases we indeed get the frequencies we started with, as well as having only nonzero real parts)

Now, what happens if we inject a signal above the highest frequency in our Fourier sum?



In [13]:

    
a_spatial_undersampled_inject = []
b_spatial_undersampled_inject = []
for i in arange(spectral_size):
    a_spatial_undersampled_inject.append(0.0)
    b_spatial_undersampled_inject.append(0.0)
    for n in arange(spectral_size):
        a_spatial_undersampled_inject[i] += a[n]*math.cos(2*pi*n*i/spectral_size)
        b_spatial_undersampled_inject[i] += b[n]*math.sin(2*pi*n*i/spectral_size)
    a_spatial_undersampled_inject[i] += 1.0*math.cos(2*pi*(spectral_size)*i/spectral_size)
    b_spatial_undersampled_inject[i] += 1.0*math.sin(2*pi*(spectral_size)*i/spectral_size)



In [14]:

    
pyplot.plot(arange(spectral_size)*1.0/spectral_size, spectral_size*real(a_ifft), 'b')
pyplot.plot(arange(spectral_size)*1.0/spectral_size, spectral_size*imag(a_ifft), 'r')
pyplot.plot(arange(spectral_size)*1.0/spectral_size, a_spatial_undersampled_inject, 'k')
pyplot.show()

We see the signal is translated up by whatever the coefficients of the injected frequency were and therfore when we try and extract the Fourier components from this spatial signal, we won't be able to tell that this frequency was added in because it just appears as a constant offset here. Similar results hold for other higher frequencies (eg. n = spectral_size + 4 would show up as n = 4).

Now we can try to multiply our two functions in space and transform them back, checking whether we get our expected results, and seeing how the outcome depends on the number of spatial points we use.



In [17]:

    
product_size = len(a_spatial)
ab_spatial = zeros(product_size)
for i in arange(len(a_spatial)):
    ab_spatial[i] = a_spatial[i]*b_spatial[i]
pyplot.plot(arange(spatial_size), a_spatial, 'b')
pyplot.plot(arange(spatial_size), b_spatial, 'r')
pyplot.plot(arange(product_size), ab_spatial, 'k')
pyplot.show()



In [10]:

    
ab_spatial_fft = numpy.fft.fft(np.array(ab_spatial))



In [11]:

    
pyplot.plot(arange(product_size), imag(ab_spatial_fft)/product_size, 'b')
pyplot.plot(arange(product_size), real(ab_spatial_fft)/product_size, 'r')
pyplot.show()



In [108]:

    
#This if for the cos*sin products
product_spectral_full = zeros(product_size)
product_spectral_full[1] = (a[0]*b[1] + 0.5*a[1]*b[2])
product_spectral_full[2] = (a[0]*b[2] + 0.5*a[1]*b[1])
product_spectral_full[3] = 0.5*a[1]*b[2]
product_spectral_full[4] = -0.5*a[6]*b[2]
product_spectral_full[5] = -0.5*a[6]*b[1]
product_spectral_full[7] = 0.5*a[6]*b[1]
product_spectral_full[8] = 0.5*a[6]*b[2]



In [14]:

    
print spectral_size, product_size
pyplot.plot(arange(product_size), -2*imag(ab_spatial_fft)/(product_size), 'b')
pyplot.plot(arange(product_size), real(ab_spatial_fft)/(product_size), 'r')
pyplot.plot(arange(product_size), product_spectral_full, 'k')
pyplot.show()

So this doesn't quite match up, but things are only off by a factor of 2. The terms agree when we multiply the fft terms by 2. I need to check whether this is equivalent to using the frequencies in the upper half of the array. Here we only had frequencies in our product that were less than the original max frequency in the a and b arrays alone. The FFTs have extra frequencies (ones that can be hidden by aliasing?)

We might be able to see what's going on here by looking at what happens when we take the iFFT of the a & b functions with our full spatial_size and see what frequencies we get out. Is it just the ones from 0 to spectral_size?



In [15]:

    
a_spatial_fft = numpy.fft.fft(np.array(a_spatial))
b_spatial_fft = numpy.fft.fft(np.array(b_spatial))



In [27]:

    
pyplot.plot(arange(spatial_size), imag(a_spatial_fft), 'b')
pyplot.plot(arange(spatial_size), real(b_spatial_fft), 'r')
pyplot.plot(arange(spatial_size), zeros(spatial_size), 'k')
pyplot.show()

(IGNORE-INCORRECT): This does not seem to work (i.e. we don't get the original frequencies back). I think this is because they're not part of the possible frequencies in this larger transform. I think we'd want something like initial frequencies from 0-10 mean you have frequency spacings of deltaf = 1/11. The only way to get these intial frequecies back is to use some multiple of the intial points. Eg. frequecies 1/11 and 6/11 can also be 2/22 and 12/22 if we use 22 points. Let's try this again with spatial_size = 2*spectral_size instead of 2*spectral_size - 1. (FIX:) I think we weren't seeing anything since I was looking at the wrong pices of the FFTs (should have plotted real(a_spatial_fft) and imag(b_spatial_fft).



In [28]:

    
a_spatial = []
b_spatial = []
for i in arange(spatial_size):
    a_spatial.append(0.0)
    b_spatial.append(0.0)
    for n in arange(spectral_size):
        a_spatial[i] += a[n]*math.cos(2*pi*n*i/(spatial_size))
        b_spatial[i] += b[n]*math.sin(2*pi*n*i/(spatial_size))

a_spatial_fft = numpy.fft.fft(np.array(a_spatial))
b_spatial_fft = numpy.fft.fft(np.array(b_spatial))
pyplot.plot(arange(spatial_size), real(a_spatial_fft)/(spatial_size), 'b')
pyplot.plot(arange(spatial_size), imag(a_spatial_fft)/(spatial_size), 'r')
pyplot.show()
pyplot.plot(arange(spatial_size), real(b_spatial_fft)/(spatial_size), 'b')
pyplot.plot(arange(spatial_size), imag(b_spatial_fft)/(spatial_size), 'r')
pyplot.show()

Aha, great! This looks like what I expected (that is - the original frequecies appear, but so do other higher ones (even in the original functions - not just the ab product). It looks like the a terms worked out to have a[0] reproduced correctly, with a[1] and a[6] being off by a factor of 2. Is this because of doubling spectral_size to spatial_size or because of a sin/cos split?



In [25]:

    
a_spatial = []
b_spatial = []
for i in arange(2*spatial_size):
    a_spatial.append(0.0)
    b_spatial.append(0.0)
    for n in arange(spectral_size):
        a_spatial[i] += a[n]*math.cos(2*pi*n*i/(2*spatial_size))
        b_spatial[i] += b[n]*math.sin(2*pi*n*i/(2*spatial_size))

a_spatial_fft = numpy.fft.fft(np.array(a_spatial))
b_spatial_fft = numpy.fft.fft(np.array(b_spatial))
pyplot.plot(arange(2*spatial_size), real(a_spatial_fft)/(2*spatial_size), 'b')
pyplot.plot(arange(2*spatial_size), imag(a_spatial_fft)/(2*spatial_size), 'r')
pyplot.show()
pyplot.plot(arange(2*spatial_size), real(b_spatial_fft)/(2*spatial_size), 'b')
pyplot.plot(arange(2*spatial_size), imag(b_spatial_fft)/(2*spatial_size), 'r')
pyplot.show()

Yep, it looks like it doesn't have anything to do with the number of points, but instead has to do with the sin/cos reconstruction. I'll have to look into that on paper later. The remaining question here is where the higher frequency components in the FFT are coming from. They're explicity not in the original functions. Do they mean something here? Can I throw them away?

I'm assuming that these will stick around even if I just use a single sin/cos term:



In [20]:

    
a_single = []
b_single = []
temp_spatial_size = spectral_size
for i in arange(temp_spatial_size):
    a_single.append(0.0)
    b_single.append(0.0)
    a_single[i] = 3*math.cos(4*pi*i/(temp_spatial_size))
    b_single[i] = 2*math.sin(6*pi*i/(temp_spatial_size))
    
a_single_fft = numpy.fft.fft(np.array(a_single))
b_single_fft = numpy.fft.fft(np.array(b_single))
pyplot.plot(arange(temp_spatial_size), real(a_single_fft)/(temp_spatial_size), 'b')
pyplot.plot(arange(temp_spatial_size), imag(a_single_fft)/(temp_spatial_size), 'r')
pyplot.show()
pyplot.plot(arange(temp_spatial_size), real(b_single_fft)/(temp_spatial_size), 'b')
pyplot.plot(arange(temp_spatial_size), imag(b_single_fft)/(temp_spatial_size), 'r')
pyplot.show()

It actually looks like this is a difference between iFFT'ing and not. When we start with frequency components, make a spatial sum out of them with the iFFT and then FFT back to get the frequency components back then everything looks good. This is what we did at the beginning of this whole section. If, instead, we take the spatial sum directly and then try to FFT back, then we get these weird higher-frequency terms. Perhaps this is what the unused real/imaginary parts of the iFFT are there to account for? Anyway, let's try it again here. One thing we need to do is explicity extend the frequency range of the initial a and b arrays - pad with zeros.



In [21]:

    
print a
a_pad = numpy.pad(a,(0,spatial_size-len(a)),'constant', constant_values=(0))
b_pad = numpy.pad(b,(0,spatial_size-len(b)),'constant', constant_values=(0))
print a_pad
print b_pad
print len(a), len(a_pad), len(b_pad)









    



[ 1.  2.  0.  0.  0.  0.  1.  0.  0.  0.  0.]
[ 1.  2.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.]
[ 0.   0.5  3.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.
  0.   0.   0.   0.   0.   0.   0.   0. ]
11 23 23



In [22]:

    
a_pad_ifft = numpy.fft.ifft(numpy.array(a_pad))
b_pad_ifft = numpy.fft.ifft(numpy.array(b_pad))



In [23]:

    
pyplot.plot(arange(spatial_size), real(a_pad_ifft)/(spatial_size), 'b')
pyplot.plot(arange(spatial_size), imag(a_pad_ifft)/(spatial_size), 'r')
pyplot.show()



In [24]:

    
pyplot.plot(arange(spatial_size), real(b_pad_ifft)/(spatial_size), 'b')
pyplot.plot(arange(spatial_size), imag(b_pad_ifft)/(spatial_size), 'r')
pyplot.show()

We should get matching with the real part of a_pad_ifft with a_spatial and the imaginary part of b_pad_ifft with b_spatial



In [34]:

    
pyplot.plot(arange(spatial_size), spatial_size*real(a_pad_ifft), 'b')
pyplot.plot(arange(spatial_size), spatial_size*imag(a_pad_ifft), 'r')
pyplot.plot(arange(spatial_size), a_spatial, 'k')
pyplot.show()



In [35]:

    
pyplot.plot(arange(spatial_size), spatial_size*real(b_pad_ifft), 'b')
pyplot.plot(arange(spatial_size), spatial_size*imag(b_pad_ifft), 'r')
pyplot.plot(arange(spatial_size), numpy.array(b_spatial), 'k')
pyplot.show()

Great, just as expected! Our iFFTs have weird imaginary/real components too but they should transform back to their frequency components without issue:



In [36]:

    
a_pad_ifft_fft = numpy.fft.fft(numpy.array(a_pad_ifft))
b_pad_ifft_fft = numpy.fft.fft(numpy.array(b_pad_ifft))



In [37]:

    
pyplot.plot(arange(spatial_size), real(a_pad_ifft_fft), 'b')
pyplot.plot(arange(spatial_size), imag(a_pad_ifft_fft), 'r')
pyplot.plot(arange(spatial_size), a_pad, 'k')
pyplot.show()



In [38]:

    
pyplot.plot(arange(spatial_size), real(b_pad_ifft_fft), 'b')
pyplot.plot(arange(spatial_size), imag(b_pad_ifft_fft), 'r')
pyplot.plot(arange(spatial_size), b_pad, 'k')
pyplot.show()

Perfect - just as expected. Now we can finally test what's going on with the products using iFFTs instead of explicit sin/cos summation.



In [39]:

    
ab_ifft = spatial_size**2*a_pad_ifft*b_pad_ifft
pyplot.plot(arange(spatial_size), real(ab_ifft), 'b')
pyplot.plot(arange(spatial_size), imag(ab_ifft), 'r')
pyplot.plot(arange(spatial_size), ab_spatial, 'k')
pyplot.show()

That was not the right function. Instead we should do ab_ifft = spatial_size**2*real(a_pad_ifft)*imag(b_pad_ifft). What to do about the imaginary part of ab_ifft?



In [40]:

    
#ab_ifft = spatial_size**2*real(a_pad_ifft)*imag(b_pad_ifft) #+ (spatial_size**2*imag(a_pad_ifft)*real(b_pad_ifft))*1j
ab_ifft = spatial_size**2*a_pad_ifft*b_pad_ifft #Product of complex arrays
pyplot.plot(arange(spatial_size), real(ab_ifft), 'b')
pyplot.plot(arange(spatial_size), imag(ab_ifft), 'r')
pyplot.plot(arange(spatial_size), ab_spatial, 'k')
pyplot.show()



In [41]:

    
ab_ifft_fft = numpy.fft.fft(numpy.array(ab_ifft))
pyplot.plot(arange(spatial_size), real(ab_ifft_fft)/(spatial_size), 'b')
pyplot.plot(arange(spatial_size), imag(ab_ifft_fft)/(spatial_size), 'r')
pyplot.plot(arange(spatial_size), product_spectral_full, 'k')
pyplot.show()



In [42]:

    
print a_pad
print b_pad
print real(ab_ifft_fft)/spatial_size
print imag(ab_ifft_fft)/spatial_size
print product_spectral_full









    



[ 1.  2.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.]
[ 0.   0.5  3.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.
  0.   0.   0.   0.   0.   0.   0.   0. ]
[ -8.68870193e-16   5.00000000e-01   4.00000000e+00   6.00000000e+00
   1.68946982e-16  -6.06278313e-15   1.15849359e-15   5.00000000e-01
   3.00000000e+00   1.03299012e-15  -4.80292135e-16   1.09936215e-15
  -4.06679521e-16  -1.85117622e-15  -1.12953125e-15   6.17863248e-16
  -1.15849359e-15  -3.28239851e-16  -1.19711004e-15  -1.11022302e-16
  -1.85358975e-15  -1.54465812e-15   2.31698718e-16]
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.]
[ 0.    3.5   3.5   3.   -1.5  -0.25  0.    0.25  1.5   0.    0.    0.    0.
  0.    0.    0.    0.    0.    0.    0.    0.    0.    0.  ]

What if we make the b frequency components imaginary so that the transform is a real function? Then we don't have to worry about multiplying the real part of one function with the imaginary part of another.



In [49]:

    
a_pad_ifft = numpy.fft.ifft(numpy.array(a_pad))
b_pad_ifft = numpy.fft.ifft(numpy.array(b_pad*-1.0j))



In [63]:

    
pyplot.plot(arange(spatial_size), spatial_size*real(a_pad_ifft), 'b')
pyplot.plot(arange(spatial_size), spatial_size*imag(a_pad_ifft), 'r')
pyplot.plot(arange(spatial_size), a_spatial, 'k')
pyplot.show()



In [64]:

    
pyplot.plot(arange(spatial_size), spatial_size*real(b_pad_ifft), 'b')
pyplot.plot(arange(spatial_size), spatial_size*imag(b_pad_ifft), 'r')
pyplot.plot(arange(spatial_size), numpy.array(b_spatial), 'k')
pyplot.show()



In [57]:

    
a_pad_ifft_fft = numpy.fft.fft(numpy.array(a_pad_ifft))
b_pad_ifft_fft = numpy.fft.fft(numpy.array(b_pad_ifft*-1.0))



In [58]:

    
pyplot.plot(arange(spatial_size), real(a_pad_ifft_fft), 'b')
pyplot.plot(arange(spatial_size), imag(a_pad_ifft_fft), 'r')
pyplot.plot(arange(spatial_size), a_pad, 'k')
pyplot.show()



In [59]:

    
pyplot.plot(arange(spatial_size), real(b_pad_ifft_fft), 'b')
pyplot.plot(arange(spatial_size), imag(b_pad_ifft_fft), 'r')
pyplot.plot(arange(spatial_size), b_pad, 'k')
pyplot.show()



In [67]:

    
#ab_ifft = spatial_size**2*a_pad_ifft*b_pad_ifft
ab_ifft = spatial_size**2*real(a_pad_ifft)*real(b_pad_ifft)
pyplot.plot(arange(spatial_size), real(ab_ifft), 'b')
pyplot.plot(arange(spatial_size), imag(ab_ifft), 'r')
pyplot.plot(arange(spatial_size), ab_spatial, 'k')
pyplot.show()



In [69]:

    
ab_ifft_fft = numpy.fft.fft(numpy.array(ab_ifft))
pyplot.plot(arange(spatial_size), real(ab_ifft_fft)/(spatial_size), 'b')
pyplot.plot(arange(spatial_size), imag(ab_ifft_fft)/(spatial_size), 'r')
pyplot.plot(arange(spatial_size), product_spectral_full, 'k')
pyplot.show()

Ok, this doesn't seem to work either because a_pad_ifft and b_pad_ifft both have real and imag parts so real(ab_ifft) contains contributions from both real(a_ifft)real(b_ifft) and -imag(a_ifft)imag(b_ifft). What can we do about this? Let's see what happens when we iFFT the expected spectral terms in our product function.



In [72]:

    
product_spectral_full_ifft = numpy.fft.ifft(numpy.array(product_spectral_full))
pyplot.plot(arange(spatial_size), spatial_size*real(product_spectral_full_ifft), 'b')
pyplot.plot(arange(spatial_size), spatial_size*imag(product_spectral_full_ifft), 'r')
pyplot.plot(arange(spatial_size), ab_spatial, 'k')
pyplot.show()

So this works - the imagiary part of the iFFT of our expected spectral terms give the correct product function in space. Now, how do we get the real part from the pieces of the individual iFFTs?

Perhaps we should be working with the discrete cosine/sine transforms?



In [5]:

    
import scipy.fftpack



In [6]:

    
a_spatial = []
b_spatial = []
for i in arange(spatial_size):
    a_spatial.append(0.0)
    b_spatial.append(0.0)
    for n in arange(spectral_size):
        a_spatial[i] += a[n]*math.cos(pi*n*i/(spatial_size-1))
        b_spatial[i] += b[n]*math.sin(pi*n*i/(spatial_size-1))
        
pyplot.plot(arange(spatial_size), a_spatial, 'b')
pyplot.plot(arange(spatial_size), b_spatial, 'r')
pyplot.show()



In [12]:

    
a_pad_idct = scipy.fftpack.idct(numpy.array(a_pad), type=1)
b_pad_idst = scipy.fftpack.idst(numpy.array(b_pad), type=1)



In [13]:

    
pyplot.plot(arange(spatial_size), a_pad_idct, 'b')
pyplot.plot(arange(spatial_size), a_spatial, 'r')
pyplot.show()



In [11]:

    
pyplot.plot(arange(spatial_size), b_pad_idst, 'b')
pyplot.plot(arange(spatial_size), b_spatial, 'r')
pyplot.show()

Realization that we have negative frequencies up in here

Ok - so as we thought, we need to convert our initial real frequency components into complex values



In [10]:

    
a_complex = zeros(spatial_size)*1j
b_complex = zeros(spatial_size)*1j
a_complex[0] = a[0]
for index, val in enumerate(a[1:]):
    a_complex[index+1] = val/2
    a_complex[spatial_size-(index+1)] = val/2
for index, val in enumerate(b[1:]):
    b_complex[index+1] = val/(2j)
    b_complex[spatial_size-(index+1)] = -val/(2j)

print a
print a_complex
print b
print b_complex









    



[ 1.  2.  0.  0.  0.  0.  1.  0.  0.  0.  0.]
[ 1.0+0.j  1.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.5+0.j  0.0+0.j
  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j
  0.0+0.j  0.5+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  1.0+0.j]
[ 0.   0.5  3.   0.   0.   0.   0.   0.   0.   0.   0. ]
[ 0.+0.j    0.-0.25j  0.-1.5j   0.+0.j    0.+0.j    0.+0.j    0.+0.j
  0.+0.j    0.+0.j    0.+0.j    0.+0.j    0.+0.j    0.+0.j    0.+0.j
  0.+0.j    0.+0.j    0.+0.j    0.+0.j    0.+0.j    0.+0.j    0.+0.j
  0.+1.5j   0.+0.25j]



In [11]:

    
a_complex_ifft = numpy.fft.ifft(numpy.array(a_complex))
b_complex_ifft = numpy.fft.ifft(numpy.array(b_complex))



In [129]:

    
pyplot.plot(arange(spatial_size), spatial_size*real(a_complex_ifft), 'b')
pyplot.plot(arange(spatial_size), spatial_size*imag(a_complex_ifft), 'r')
pyplot.plot(arange(spatial_size), a_spatial, 'k')
pyplot.show()



In [128]:

    
pyplot.plot(arange(spatial_size), spatial_size*real(b_complex_ifft), 'b')
pyplot.plot(arange(spatial_size), spatial_size*imag(b_complex_ifft), 'r')
pyplot.plot(arange(spatial_size), b_spatial, 'k')
pyplot.show()

Fuckin' A - finally! Now let's multiply this ish.



In [14]:

    
ab_complex_ifft_fft = numpy.fft.fft(numpy.array(spatial_size**2*a_complex_ifft*b_complex_ifft))



In [15]:

    
pyplot.plot(arange(spatial_size), real(ab_complex_ifft_fft)/spatial_size, 'b')
pyplot.plot(arange(spatial_size), imag(ab_complex_ifft_fft)/spatial_size, 'r')
pyplot.plot(arange(spatial_size), product_spectral_full, 'k')
pyplot.show()









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-15-d86649c33a76> in <module>()
      1 pyplot.plot(arange(spatial_size), real(ab_complex_ifft_fft)/spatial_size, 'b')
      2 pyplot.plot(arange(spatial_size), imag(ab_complex_ifft_fft)/spatial_size, 'r')
----> 3 pyplot.plot(arange(spatial_size), product_spectral_full, 'k')
      4 pyplot.show()

NameError: name 'product_spectral_full' is not defined

Finally - negative frequencies are still there, but now I understand why. They must be for real spatial functions!

Let's look at what happens when there are now frequencies in the product that are above the original frequency maximum of our a and b functions. We'll do this by replacing (eg. b[2] with b[6], so that we get a non-zero product_spectral_full[12] term)



In [99]:

    
a = zeros(spectral_size)
b = zeros(spectral_size)
a[0] = 1.0
a[1] = 2.0
a[6] = 1.0
b[1] = 0.5
b[6] = 3.0
print a, b









    



[ 1.  2.  0.  0.  0.  0.  1.  0.  0.  0.  0.] [ 0.   0.5  0.   0.   0.   0.   3.   0.   0.   0.   0. ]



In [69]:

    
a_spatial = []
b_spatial = []
ab_spatial = []
for i in arange(spatial_size):
    a_spatial.append(0.0)
    b_spatial.append(0.0)
    ab_spatial.append(0.0)
    for n in arange(spectral_size):
        a_spatial[i] += a[n]*math.cos(2*pi*n*i/(spatial_size))
        b_spatial[i] += b[n]*math.sin(2*pi*n*i/(spatial_size))
        ab_spatial[i] += a_spatial[i]*b_spatial[i]



In [100]:

    
a_complex = zeros(spatial_size)*1j
b_complex = zeros(spatial_size)*1j
a_complex[0] = a[0]
for index, val in enumerate(a[1:]):
    a_complex[index+1] = val/2
    a_complex[spatial_size-(index+1)] = val/2
for index, val in enumerate(b[1:]):
    b_complex[index+1] = val/(2j)
    b_complex[spatial_size-(index+1)] = -val/(2j)

print a
print a_complex
print b
print b_complex









    



[ 1.  2.  0.  0.  0.  0.  1.  0.  0.  0.  0.]
[ 1.0+0.j  1.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.5+0.j  0.0+0.j
  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j
  0.0+0.j  0.5+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  1.0+0.j]
[ 0.   0.5  0.   0.   0.   0.   3.   0.   0.   0.   0. ]
[ 0.+0.j    0.-0.25j  0.+0.j    0.+0.j    0.+0.j    0.+0.j    0.-1.5j
  0.+0.j    0.+0.j    0.+0.j    0.+0.j    0.+0.j    0.+0.j    0.+0.j
  0.+0.j    0.+0.j    0.+0.j    0.+1.5j   0.+0.j    0.+0.j    0.+0.j
  0.+0.j    0.+0.25j]



In [63]:

    
a_complex_ifft = numpy.fft.ifft(numpy.array(a_complex))
b_complex_ifft = numpy.fft.ifft(numpy.array(b_complex))



In [140]:

    
pyplot.plot(arange(spatial_size), spatial_size*real(a_complex_ifft), 'b')
pyplot.plot(arange(spatial_size), spatial_size*imag(a_complex_ifft), 'r')
pyplot.plot(arange(spatial_size), a_spatial, 'k')
pyplot.show()



In [139]:

    
pyplot.plot(arange(spatial_size), spatial_size*real(b_complex_ifft), 'b')
pyplot.plot(arange(spatial_size), spatial_size*imag(b_complex_ifft), 'r')
pyplot.plot(arange(spatial_size), b_spatial, 'k')
pyplot.show()



In [70]:

    
ab_complex_ifft_fft = numpy.fft.fft(numpy.array(spatial_size**2*a_complex_ifft*b_complex_ifft))



In [19]:

    
product_spectral_full = zeros(product_size)
product_spectral_full[1] = (a[0]*b[1] + 0.5*a[1]*b[2])
product_spectral_full[2] = (a[0]*b[2] + 0.5*a[1]*b[1])
product_spectral_full[3] = 0.5*a[1]*b[2]
product_spectral_full[4] = -0.5*a[6]*b[2]
product_spectral_full[5] = -0.5*a[6]*b[1] + 0.5*a[1]*b[6]
product_spectral_full[6] = a[0]*b[6]
product_spectral_full[7] = 0.5*a[6]*b[1] + 0.5*a[1]*b[6]
product_spectral_full[8] = 0.5*a[6]*b[2]
product_spectral_full[12] = 0.5*a[6]*b[6]



In [61]:

    
pyplot.plot(arange(spatial_size), real(ab_complex_ifft_fft)/spatial_size, 'b')
pyplot.plot(arange(spatial_size), imag(ab_complex_ifft_fft)/spatial_size, 'r')
pyplot.plot(arange(spatial_size), product_spectral_full, 'k')
pyplot.show()



In [84]:

    
print spectral_size, spatial_size

So here's how you multiply terms in the spectral transform method. You start out with Fourier terms representing sines/cosines. These need to be converted to coefficients in a complex exponetial expansion for use with the FFT. This is pretty straightforward. Here's how we do it for our a and b arrays defined at the beginning:



In [109]:

    
a_complex = zeros(spatial_size)*1j
b_complex = zeros(spatial_size)*1j
a_complex[0] = a[0]
for index, val in enumerate(a[1:]):
    a_complex[index+1] = val/2
    a_complex[spatial_size-(index+1)] = val/2
for index, val in enumerate(b[1:]):
    b_complex[index+1] = val/(2j)
    b_complex[spatial_size-(index+1)] = -val/(2j)

print a
print a_complex
print b
print b_complex









    



[ 1.  2.  0.  0.  0.  0.  1.  0.  0.  0.  0.]
[ 1.0+0.j  1.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.5+0.j  0.0+0.j
  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j
  0.0+0.j  0.5+0.j  0.0+0.j  0.0+0.j  0.0+0.j  0.0+0.j  1.0+0.j]
[ 0.   0.5  3.   0.   0.   0.   0.   0.   0.   0.   0. ]
[ 0.+0.j    0.-0.25j  0.-1.5j   0.+0.j    0.+0.j    0.+0.j    0.+0.j
  0.+0.j    0.+0.j    0.+0.j    0.+0.j    0.+0.j    0.+0.j    0.+0.j
  0.+0.j    0.+0.j    0.+0.j    0.+0.j    0.+0.j    0.+0.j    0.+0.j
  0.+1.5j   0.+0.25j]

The next step is to multiply the functions together in spatial space via the iFFT (noting the normalization factor in the FFT definition used by SciPy - this is not in FFTW):



In [110]:

    
a_complex_ifft = numpy.fft.ifft(numpy.array(a_complex))
b_complex_ifft = numpy.fft.ifft(numpy.array(b_complex))

Finally, we transform the product back to coefficients in a complex exponential expansion, and then back to sine/cosine coefficients if that's what we started with:



In [111]:

    
ab_complex_ifft_fft = numpy.fft.fft(numpy.array(spatial_size**2*a_complex_ifft*b_complex_ifft))
ab = zeros(spectral_size)
ab = -2*imag(ab_complex_ifft_fft[0:spectral_size])/(spatial_size)
print ab









    



[ -0.00000000e+00   3.50000000e+00   3.50000000e+00   3.00000000e+00
  -1.50000000e+00  -2.50000000e-01  -1.70877805e-15   2.50000000e-01
   1.50000000e+00   2.22044605e-15   1.43363582e-15]

Check that everything worked as planned:



In [112]:

    
pyplot.plot(arange(spatial_size), real(ab_complex_ifft_fft)/spatial_size, 'b')
pyplot.plot(arange(spatial_size), imag(ab_complex_ifft_fft)/spatial_size, 'r')
pyplot.plot(arange(spatial_size), product_spectral_full, 'k')
pyplot.show()



In [152]:

    
a_spatial = []
b_spatial = []
ab_spatial = []
ab_analytic = []
ab_transform = []
for i in arange(spatial_size):
    a_spatial.append(0.0)
    b_spatial.append(0.0)
    ab_spatial.append(0.0)
    ab_analytic.append(0.0)
    ab_transform.append(0.0)
    for n in arange(spectral_size):
        a_spatial[i] += a[n]*math.cos(2*pi*n*i/(spatial_size))
        b_spatial[i] += b[n]*math.sin(2*pi*n*i/(spatial_size))
        ab_analytic[i] += product_spectral_full[n]*math.sin(2*pi*n*i/(spatial_size))
        ab_transform[i] += ab[n]*math.sin(2*pi*n*i/(spatial_size))
    ab_spatial[i] += a_spatial[i]*b_spatial[i]

pyplot.plot(arange(spatial_size), ab_analytic, 'b')
pyplot.plot(arange(spatial_size), ab_transform, 'r')
pyplot.plot(arange(spatial_size), ab_spatial, 'k')
pyplot.plot(arange(spatial_size), numpy.array(spatial_size**2*real(a_complex_ifft)*real(b_complex_ifft)), 'o')
pyplot.show()



In [ ]: