[Numpy-discussion] Interpolation via Fourier transform

[Nadav Horesh]

I found the solution I needed for my peculiar case after reading your email based of the following stages:

I have a N x N frequency-domain matrix Z

1. Use fftshift to obtain a DC centered matrix
   Note: fftshift(fft(a)) replaces np.fft.fft(np.power(-1,np.arange(64))*a)
   Zs = np.fft.fftshift(Z)
2. pad Zs with zeros
   scale = int(ceil(float(N)/M))
   MM = scale*M
   Ztemp = np.zeros((MM,MM), dtype=complex)
   Ztemp[(MM-N)//2:(N-MM)//2,(MM-N)//2:(N-MM)//2] = Zs
3. Shift back to a "normal order"
   Ztemp = np.fft.ifftshift(Ztemp)
4. Transform to the "time domain" and sub-sample 
   z = np.fft.ifft2(Ztemp)[::scale, ::scale]

I went this was since I needed the aliasing, otherwise I could just truncate Zs to size MxM.

 Thank you,

   Nadav.


-----הודעה מקורית-----
מאת: numpy-discussion-bounces at scipy.org בשם M Trumpis
נשלח: ה 05-מרץ-09 21:51
אל: Discussion of Numerical Python
נושא: Re: [Numpy-discussion] Interpolation via Fourier transform
 
Hi Nadav.. if you want a lower resolution 2d function with the same
field of view (or whatever term is appropriate to your case), then in
principle you can truncate your higher frequencies and do this:

sig = ifft2_func(sig[N/2 - M/2:N/2 + M/2, N/2 - M/2:N/2+M/2])

I like to use an fft that transforms from an array indexing
negative-to-positive freqs to an array that indexes
negative-to-positive spatial points, so in both spaces, the origin is
at (N/2,N/2). Then the expression works as-is.

The problem is if you've got different indexing in one or both spaces
(typically positive frequencies followed by negative) you can play
around with a change of variables in your DFT in one or both spaces.
If the DFT is defined as a computing frequencies from 0,N, then
putting in n' = n-N/2  leads to a term like exp(1j*pi*q) that
multiplies f[q]. Here's a toy example:

a = np.cos(2*np.pi*5*np.arange(64)/64.)

P.plot(np.fft.fft(a).real)

P.plot(np.fft.fft(np.power(-1,np.arange(64))*a).real)

The second one is centered about index N/2

Similarly, if you need to change the limits of the summation of the
DFT from 0,N to -N/2,N/2, then you can multiply exp(1j*pi*n) to the
outside of the summation.

Like I said, easy enough in principle!

Mike

On Thu, Mar 5, 2009 at 11:02 AM, Nadav Horesh <nadavh at visionsense.com> wrote:
>
> I apology for this off topic question:
>
>  I have a 2D FT of size N x N, and I would like to reconstruct the original signal with a lower sampling frequency directly (without using an interpolation procedure): Given M < N the goal is to compute a M x M "time domain" signal.
>
>  In the case of 1D signal the trick is simple --- given a length N freq. domain Sig:
>
>  sig = np.fft.ifft(Sig, M)
>
> This trick does not work in 2D:
>
>  sig = np.fft.ifft2(Sig, (M,M))
>
> is far from being the right answer.
>
>  Any ideas?
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion at scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

-------------- next part --------------
A non-text attachment was scrubbed...
Name: winmail.dat
Type: application/ms-tnef
Size: 4579 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20090306/e3388ca9/attachment.bin>