How to separate the upwards-propagating from the downwards-propagating waves?

Question

I've got a real 2D image of propagating waves in both directions. How can I separate this image into one with the upward-propagating waves and one with the downward-propagating ones?

This example has one wave in each direction. The x-scale shows time, the y-scale shows position, and intensity shows the amplitude of the wave at that position and time. The first image is my input.

However, I'm not very used to 2D signal processing, so the way I'm doing this is by just setting two quadrants of the FFT spectrum to zero and then reverse-transforming back. I got the other two images using the following algorithm:

X = fftshift(fft2( x )); % x is the input (image above)
X( 1:floor(size(X,1)/2),    1:floor(size(X,2)/2)   ) = 0; % null 1st quadrant
X( ceil(size(X,1)/2)+1:end, ceil(size(X,2)/2)+1:end) = 0; % null 3rd quadrant
y = ifft2( ifftshift(X) );
y = real( y ); % y is the output (image above)

It feels like such a barbaric method would introduce all kinds of artefacts, so I'm looking for feedback from someone more experienced.

• Is there a better way to achieve this?
• How should I "soften" the edges of the 0/1 "mask" I'm multiplying the FFT with?
• The output is complex, so I'm discarding the imaginary part. Is that correct?
• I'm also a bit worried about the fact that the noise will not be white anymore, but I suppose there isn't much to do about that except for low-pass filtering the whole image.

Answer 1

Instead of nulling the 2 quadrants, multiply them by $i$. Then the real part will be one image and the imaginary part the other image.

X = fftshift(fft2( x )); % x is the input (image above)
X( 1:floor(size(X,1)/2),    1:floor(size(X,2)/2)   ) = ...
    X( 1:floor(size(X,1)/2),    1:floor(size(X,2)/2)   ) * 1i;
X( ceil(size(X,1)/2)+1:end, ceil(size(X,2)/2)+1:end) = ...
    X( ceil(size(X,1)/2)+1:end, ceil(size(X,2)/2)+1:end) * 1i;
y = ifft2( ifftshift(X) );
imagesc(real(y)); pause;
imagesc(imag(y)); pause;

Since this is unitary operation the energy of this complex image will be the same as the original, as will the noise energy. Doesn't look very white to me though, seems there is differences between entire rows.

If you want to get rid of the artefacts, you could use a smoother function such as cosine of the angle of the spectrum:

[rows,cols] = size(x);
IM = fft2(x);
% Create polar co-ordinate of FFT
[ux, uy] = meshgrid(([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...
    ([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2)));
ux = ifftshift(ux);   % Quadrant shift to put 0 frequency at the corners
uy = ifftshift(uy);
th = atan2(uy,ux);
r = sqrt(ux.^2 + uy.^2);
% cos^2
R2 = cos(th +pi/4).^2 + 1i*cos(th +3*pi/4).^2;
R2(1,1) = 0;
I_a = real(ifft2(R2.*IM));
I_b = imag(ifft2(R2.*IM));
imagesc(I_a); pause;
imagesc(I_b); pause;

I don't think the energy remains constant in this case, but assuming white noise energy the noise of the resulting two images should be similar.

Using the polar form of the spectrum it will be easier to experiment and create smoother cutoffs. For example, you could 'narrow' the angle of the filter to isolate the wave further.

Answer 2

1. I would expect that this is (at least roughly) the best way to do this.
2. If you find that you need to soften the edges, you could use a cosine roll-off.
3. I'm not sure about which bit to discard. I would look at the data and see where most of the energy lies - compare just taking the real part to just taking the magnitude (signed magnitude I guess), and see which version accounts for the most energy in the image.
4. Are you sure the noise will no longer be white? I think it might still be. At least along the radial lines you are keeping you haven't changed the noise profile.

Your results look pretty good to me.