2D FFT: problems with fftshift + setting up wavenumbers

Hi,

My i/p data is an (NxN) matrix, whose elements S(x,y) describe some kind intensity over a 2D region.

I want to perform a 2D Fourier transform on this data, which I tried to do as follows:

F = fft2(S); F2 = fftshift(F); F3=abs(F3); imshow(abs(F3));

I tried to set the wavevectors as below:

N=50; kx1 =mod(1/2+(0:(N-1))/N, 1)-1/2; kx = kx1*(2*pi/delta);

ky1 = mod(1/2 + (0:(N-1))/N,1)-1/2; ky = ky1*(2*pi/delta);

[p,q]=meshgrid(kx,ky);

Now when i try to do a contour plot using h2=contourf(p,q,F3); i find that(kx,ky)=0 is not at the center. If I use fftshift a second time h2=contourf(p,q,fftshift(F3)); then get the DC component at the center. But please note that F3 is already fftshift-ed.Can someone please tell me what is going on?

Secondly, I want to keep the wavenumbers only upto the Nyquist wavenumber and throw away the symmetrically redundant ones. How does one do this?

Thirdly, how i can be sure that when i plot contourf(p,q,F3), i am plotting F3 at the corresponding wavenumbers?

Answer 1

[David Young]

You've used fftshift to put the zero-frequency component at the centre of the array F3. However, your wavevector matrices kx and ky have zero as the first component, so they do not correspond to the layout of F3. They need to have zero in the middle.

countourf then sorts the elements of all the matrices. Because kx and ky have negative values to the right of positive values, this has the effect of carrying out another fftshift on F3, before it is plotted.

You could omit the first call to fftshift. Alternatively, you could define kx and ky correctly for the centred FFT, like this:

kx = linspace(-pi/delta, (pi-2*pi/N)/delta, N);
ky = kx;
[p,q]=meshgrid(kx,ky);

I hope that answers your first question. I think it also answers your third one: contourf always plots F3(m,n) at the coordinates x=p(m,n), y=q(m,n). You just have to make sure p and q are correct.

On your second question: if your data are always real, you can discard the right half, the left half, the bottom half or the top half of your plot if you don't want to see any redundant information. It's an arbitrary choice. You just use the colon operator to select the parts of the arrays that you want to plot, for example

contourf(p(:,N/2+1:end), q(:,N/2+1:end), F3(:,N/2+1:end));

(though note that this particular option omits the Nyquist frequency itself).