Subtract background from 2D signal using spectral methods

Question

I have a 2D signal, with intensity as a function of two angles (alpha and beta), as shown below.

This contains a background, which makes the signal's base-surface locally convex-concave. I would like to identify and subtract the background such that the bottom of the signal becomes approximately planar and further processing can be done.

My initial thoughts were to fit row-wise/column-wise baselines and individually subtract the background from each line. However, I think signal processing methods may be better suited for this task. I have calculated the FFT for the signal using SciPy/Python, but don't find any obvious frequencies that I should remove.

I'm not even very clear about whether I should be looking to remove the highest frequencies, the lowest frequencies, or both. What would be a good way to proceed from here, using Python? Is the power spectral density better to look at rather than the FFT itself?

EDIT : The signal comes from a diffraction experiment, where there is a source, a specimen and a detector. The specimen is tilted from 15 degree to 90 degree in steps of 5 degree. At each tilt, it is also rotated through 360 degrees in steps of 5 degree. The tilt and rotation are the angles alpha and beta shown in the figure. Corresponding to each alpha/beta combination, the detector receives some signal, which we are plotting as intensity.

Answer 1

There are certainly ways to "kill" the underlying signal but in this case, it is probably better to do it gracefully.

I'm not even very clear about whether I should be looking to remove the highest frequencies, the lowest frequencies, or both.

In general, completely supressing a (usually) low frequency background oscillation is done by using a high pass filter. Other ways this can be achieved is by applying a median filter to extract the variation and then subtract it from the original signal, or even peforming a low order polynomial fit to the data and then extract the surface this polynomial describes. For slightly more information on those, please see this response.

Here however, the background signal is not just any nuisance but it is related to the measurement. Therefore, if we just try to chuck it away, we may be introducing a problem to the measurements we are trying to get to. So, the better approach would be to see if it is possible to take the component out somewhat more reasonably.

More so, when you think about the unknowns, which are not many. In the "Beta" dimension, we have 1 "revolution". The same applies for the "Alpha" dimension, but there, we also have a "known" phase shift and so only observe part of the sinusoid. In the "Intensity" dimension, we have a mean background over which the counts vary. So, we know the frequency of the sinusoid as well as its phase (for the "Alpha" dimension) and we would be looking at its amplitude.

In 1 dimension, the "fit" looks like this (here in GNU Octave but easily translatable to other platforms too):

Fs=256; % Sampling frequency in Hz
T=1; % Duration of the signal in seconds
t=0:(1./Fs):(T-(1./Fs)); % Time vector
p=2.*pi.*t; % Phase vector

% Setting up the "signal"
average_count=30; % Self explanatory
variation_amp=6; % The amplitude of the sinusoidal variation
sig_amp=4; % The amplitude of what here simulates the signal.
variation_phase=pi./0.8; % The phase of the sinusoidal variation.

s = average_count + variation_amp.*sin(p+variation_phase) + sig_amp.*rand(size(p));

% At this point, we know nothing about the parameters, except
% the frequency of the sinusoid, because we know that due to 
% the experimental conditions, there has to be sinusoidal
% variation proportional to an angle that we control.

signal_mean = sum(s)
signal_amp = sum(s.*exp(-j.*p));
v = (signal_mean + signal_amp.*exp(j.*p) + conj(signal_amp).*exp(j.*-p))./length(p);
% Recovery

u = s - real(v);

% Display

subplot(311);plot(s);grid on;title("Signal");xlabel("Discrete time (sample)");ylabel("Amplitude");
subplot(312);plot(real(v));grid on;title("Recovered variation");xlabel("Discrete time (sample)");ylabel("Amplitude");
subplot(313);plot(u);grid on;title("Adjusted signal");xlabel("Discrete time (sample)");ylabel("Amplitude");

This produces:

Where, basically, what we have recovered after subtracting the predictable part is that additive signal which here looks like noise but could be any remaining "form".

This is basically Discrete Fourier Transform (DFT) but just for 1 component, the one we know (or, can predict) is in the signal. These lines:

signal_mean = sum(s)
signal_amp = sum(s.*exp(-j.*p));

Perform the forward transform. And this line:

v = (signal_mean + signal_amp.*exp(j.*p) + conj(signal_amp).*exp(j.*-p))./length(p);

Performs the inverse transform but just for one component. Notice here that signal_amp recovers both the amplitude and the phase of the sinusoidal we are after.

A similar approach would be followed in two dimensions but with the necessary adjustments:

Fs=256; % Spatial sampling frequency in Lines per unit of length
T=1; % Length of image in units of length (assumed square here for simplicity)
t=0:(1./Fs):(T-(1./Fs)); % **Space** vector
p=2.*pi.*t; % Spatial phase vector
[Xp, Yp] = meshgrid(p); % Spatial phase vectors for each dimension.

% Setting up the "signal"
average_count=30; % Self explanatory
variation_amp=6; % The amplitude of the sinusoidal variation
sig_amp=4; % The amplitude of what here simulates the signal.
% The phases of the sinusoidal variation.
variation_phase_x=pi./1.2; 
variation_phase_y=pi./0.8;

s = average_count + variation_amp.*(sin(Xp+variation_phase_x)+sin(Yp+variation_phase_y)) + sig_amp.*rand(size(Xp));

rec_average = sum(sum(s));
Ay = sum(s.*exp(-j.*Yp))
Ax = sum(s.*exp(-j.*Xp),2)

rec = rec_average./prod(size(s)) + (Ax.*exp(j.*Xp) + conj(Ax).*exp(j.*-Xp) + Ay.*exp(j.*Yp) + conj(Ay).*exp(j.*-Yp))./256;

subplot(311);title("Original");surf(s);xlabel("Dim X");ylabel("Dim Y");zlabel("Amp");
subplot(312);title("Component");surf(rec);xlabel("Dim X");ylabel("Dim Y");zlabel("Amp");
subplot(313);title("Adjusted");surf(s-rec);xlabel("Dim X");ylabel("Dim Y");zlabel("Amp");

This produces:

Where as you can see, the variation is gone and only the additive component has remained. The two snippets have an almost one-to-one correspondence to show both the similarities but also adjustments you have to apply when considering higher dimensions.

In your case, since you are dealing with "counts" that must remain positive, you can ommit the DC (the first component that is simply the average across the signal) and only remove the variation.

Hope this helps.