Shift theorem in Discrete Fourier Transform

问题:

I'm trying to solve a problem with python+numpy in which I've some functions of type that I need to convolve with another function . In order to optimize code, I performed the fft of f and g, I multiplied them and then I performed the inverse transformation to obtain the result.

As a further optimization I realized that, thanks to the shift theorem, I could simply compute once the fft of f(x,y,z) and then multiply it by a phase factor that depends on to obtain the fft of . In particular, , where N is the length of both x and y.

I tried to implement this simple formula with python+numpy, but it fails for some reason that is obscure for me at the moment, so I'm asking the help of SO community in order to figure out what I'm missing.

I'm providing also a simple example.

In [1]: import numpy as np
In [2]: x = np.arange(-10, 11)
In [3]: base = np.fft.fft(np.cos(x))
In [4]: shifted = np.fft.fft(np.cos(x-1))
In [5]: w = np.fft.fftfreq(x.size)
In [6]: phase = np.exp(-2*np.pi*1.0j*w/x.size)
In [7]: test = phase * base
In [8]: (test == shifted).all()
Out[8]: False
In [9]: shifted/base
Out[9]:
array([ 0.54030231 -0.j        ,  0.54030231 -0.23216322j,
        0.54030231 -0.47512034j,  0.54030231 -0.7417705j ,
        0.54030231 -1.05016033j,  0.54030231 -1.42919168j,
        0.54030231 -1.931478j  ,  0.54030231 -2.66788185j,
        0.54030231 -3.92462627j,  0.54030231 -6.74850534j,
        0.54030231-20.55390586j,  0.54030231+20.55390586j,
        0.54030231 +6.74850534j,  0.54030231 +3.92462627j,
        0.54030231 +2.66788185j,  0.54030231 +1.931478j  ,
        0.54030231 +1.42919168j,  0.54030231 +1.05016033j,
        0.54030231 +0.7417705j ,  0.54030231 +0.47512034j,
        0.54030231 +0.23216322j])
In [10]: np.abs(shifted/base)
Out[10]:
array([  0.54030231,   0.58807001,   0.71949004,   0.91768734,
         1.18100097,   1.52791212,   2.00562555,   2.72204338,
         3.96164334,   6.77009977,  20.56100612,  20.56100612,
         6.77009977,   3.96164334,   2.72204338,   2.00562555,
         1.52791212,   1.18100097,   0.91768734,   0.71949004,   0.58807001])

I expect that by means of shifted/base I could obtain the corresponding values of the phase factor, but as could be seen, it cannot be a phase factor, since its np.abs is >= 1!

回答1:

The problem in my code was both on input line 6, due to an incorrect (my fault) interpretation of the return value of np.fft.fftfreq(), and on the necessity to pad arrays in order to obtain sound results.

The following code works great and could be extended to multidimension.

In [1]: import numpy as np
In [2]: shift = 1
In [3]: dx = 0.5
In [4]: pad = 20
In [5]: x = np.arange(-10, 11, dx)
In [6]: y = np.cos(x)
In [7]: y = np.pad(y, (0,pad), 'constant')
In [8]: y_shift = np.cos(x-shift)
In [9]: y_fft = np.fft.fft(y)
In [10]: w = np.fft.fftfreq(y.size, dx)
In [11]: phase = np.exp(-2.0*np.pi*1.0j*w*shift)
In [12]: test = phase * y_fft
In [13]: # we use np.real since the resulting inverse fft has small imaginary part values that are zero
In [14]: inv_test = np.real(np.fft.ifft(test))
In [15]: np.allclose(y[:-pad-2],inv_test[2:-pad])
Out[15]: True

回答2:

Nice, thank you very much for sharing! I've implemented something in this lines some time ago, but couldn't really grasp the mathematics of it (I've blindly ported a whitepaper that described the algorithm). FWIW, this it it: https://github.com/creaktive/flare/blob/master/nrf905_demod.c#L376