# Subtract background from 2D signal using spectral methods

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.

• The decision on which method to go for might be slightly affected by context. Is it possible to talk a little bit more about what the signal is coming from please?
– A_A
Feb 5, 2020 at 10:02
• @A_A : Thanks for your response! I have added this information to the question. Feb 5, 2020 at 10:55
• Then, assuming that you still don't "like" the underlying signal, it sounds like you can factor it out by fitting a two dimensional cosine surface adapted to your experimental conditions (?). There definitely are ways to "kill it", but what if part of it is useful too? Another way to take that signal out would be to do a "dry" run, so that you get the background count of the specimen holder which you can then subtract from the specimen run. That should take out the background bit and leave you only with the signal that is due to the specimen. Or, am I missing something?
– A_A
Feb 5, 2020 at 11:56
• @A_A: Yes, if I can separate the two, that would be great, because the characteristics of the background could have value for different reasons. The 'dry run' is what we usually do, but there is a part of specimen-generated background as well, which we call fluorescence, and for the present purpose that is not helpful. How could I go about fitting a 2D cosine surface? Feb 5, 2020 at 13:15

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");



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");


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.

• Thanks a lot! To ensure that I've understood correctly, the idea is to identify the frequencies that belong to the background in the Fourier space, and subtract them from the main signal in the inverse-space? The discrimination of frequencies, however is still done based on human judgement (or is there a way to calculate that as well?) Feb 16, 2020 at 13:18
• Special thanks for developing intuition with the 1D example and progressing to 2D! Feb 16, 2020 at 13:19
• @AppliedAcademic Glad you found it helpful. You can mark the answer as "accepted" from the check sign control on the left of the answer and this will stop it from circulating the board as "unanswered". Yes, that is the idea but it's not so much "human judgement" as what is known (or predictable) from the experimental conditions.
– A_A
Feb 16, 2020 at 15:36