Re: [Numpy-discussion] Interpolation via Fourier transform
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-boun...@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 nad...@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@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion winmail.dat___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Interpolation via Fourier transform
Hi Nadav You can also read the interesting discussion at http://projects.scipy.org/numpy/ticket/748 which also contains some padding code. I still disagree with the conclusion, but oh well :) Cheers Stéfan 2009/3/6 Nadav Horesh nad...@visionsense.com: 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. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Interpolation via Fourier transform
It was one of the first things I tried, without success Nadav. -הודעה מקורית- מאת: numpy-discussion-boun...@scipy.org בשם Anne Archibald נשלח: ה 05-מרץ-09 22:06 אל: Discussion of Numerical Python נושא: Re: [Numpy-discussion] Interpolation via Fourier transform 2009/3/5 M Trumpis mtrum...@berkeley.edu: 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! There's also the hit-it-with-a-hammer approach: Just downsample in x then in y, using the one-dimensional transforms. Anne ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion winmail.dat___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Interpolation via Fourier transform
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@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
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 nad...@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@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Interpolation via Fourier transform
2009/3/5 M Trumpis mtrum...@berkeley.edu: 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! There's also the hit-it-with-a-hammer approach: Just downsample in x then in y, using the one-dimensional transforms. Anne ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
