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 |
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) 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 |
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 |
Free forum by Nabble | Edit this page |