[SciPy-User] Shift theorem in Discrete Fourier Transform

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

[SciPy-User] Shift theorem in Discrete Fourier Transform

Domenico Delle Side
I'm trying to solve a problem with python+numpy in which I've some functions
of type $f(x-x_i,y-y_i,z)$ that I need to convolve with another function
$g(x,y,z,t)$. 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 $(x_i, y_i)$ to obtain the fft of
$f(x-x_i,y-y_i,z)$. In particular, $\mathcal{F}(f(x-x_i,y-y_i,z)) =
exp^{-2\pi j (x_i w_1 + y_i w_2)/N} \mathcal{F}(f(x,y,z))$, 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
the group 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!

Any help is greatly appreciated.
Thanks in advance,
Domenico Delle Side

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: Shift theorem in Discrete Fourier Transform

Lev Givon
Received from Domenico Delle Side on Thu, Apr 28, 2016 at 08:36:10AM EDT:
> I'm trying to solve a problem with python+numpy in which I've some functions
> of type $f(x-x_i,y-y_i,z)$ that I need to convolve with another function
> $g(x,y,z,t)$. 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.

FYI, there is an fftconvolve routine in scipy.signal that you might want to try.

> 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 $(x_i, y_i)$ to obtain the fft of
> $f(x-x_i,y-y_i,z)$. In particular, $\mathcal{F}(f(x-x_i,y-y_i,z)) =
> exp^{-2\pi j (x_i w_1 + y_i w_2)/N} \mathcal{F}(f(x,y,z))$, 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
> the group 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)
You shouldn't be dividing by x.size in the above line because of the way fftfreq
is defined; see the output of the attached modified script.

Although it isn't much of an issue in the above example, in general you should
zero-pad the signal being phase shifted so that the result of the circular
convolution effectively performed by fft is equivalent to ordinary convolution.
--
Lev Givon
Bionet Group | Neurokernel Project
http://lebedov.github.io/
http://neurokernel.github.io/

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user

fft_shift.py (600 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Shift theorem in Discrete Fourier Transform

Domenico Delle Side
Lev Givon <[hidden email]> writes:

> Received from Domenico Delle Side on Thu, Apr 28, 2016 at 08:36:10AM EDT:
>> I'm trying to solve a problem with python+numpy in which I've some
>> functions of type $f(x-x_i,y-y_i,z)$ that I need to convolve with
>> another function $g(x,y,z,t)$. 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.
>
> FYI, there is an fftconvolve routine in scipy.signal that you might
> want to try.

Thank you Lev, for your reply. I'm solving an heat equation with a
variable number (it could be >= 1000) of gaussian sources, shifted on
the x,y plane. So, I'm trying to optimize the code, that would be
otherwise slow. For this reason, I'm looking at methods that save
computing time. My question arises from this need. Instead of computing
the fft for each source, it would be easier to compute the base fft only
at start and then obtaining the desired fft trhough the shift theorem
(it cost only a multiplication).

> You shouldn't be dividing by x.size in the above line because of the
> way fftfreq is defined; see the output of the attached modified
> script.
>
> Although it isn't much of an issue in the above example, in general
> you should zero-pad the signal being phase shifted so that the result
> of the circular convolution effectively performed by fft is equivalent
> to ordinary convolution.

Your example has been very enlightening for me! Could you please be more
specific on the padding trick? I'm courious. For example, why p=16?
Otherwise have you any documentation pointer about?


Thank you again,
Nico

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user