Removing step-function like noise from signal, knowing when the step functions occur

Question

So I currently have a system which is fit for very short time measurements. Whenever a measurement is long, however, there is a need to continuously update the reference, which, each time, will produce a small step-function noise on the signal. The problem is explained in the followin paragraphs, but as a summary, the noise I want to remove consists of a series of step functions, which I know exactly when occur, but have no knowledge of the amplitude.

For reference, see the following pictures.

Whenever we update it, we accumulate any previous error. If we call $\hat{s_{xy}}$ the measurement at instant y, relative to some reference taken at instant x, starting the measurement at instant 0 and stopping at instant 3, with an update on 2 would represent the following $$ \hat{s_{01}}= s_{01}+n_{01} $$ $$ \hat{s_{02}}= s_{02}+n_{02} $$ update!! $$ \hat{s_{03}}=(s_{23}+s_{02})+n_{02}+n_{23} $$

From here on out, the term $n_{02$ will always contaminate the signal as a step-function, which will propagate as 1/f noise. (see Fig.3)

As such, the real measurement will look like the following:

Obviously, we know exactly when the signal was updated, and we can somehow control that (update more, or less often, with a few trade-offs). The question is whetherthere is any way, through post-processing, to estimate the step-function noise (depicted in the next picture), using only the real acquisition and knowledge of when the updates take place. (in the given code, this would be sig_total, and refUpdates).

Ideally, we would like to be able to do this with just that information.

Matlab code for the basic simulation presented here:

nsamp = 1000; %number of samples
t = 1:nsamp;
[s,sraw,nraw,sig_total,noise_update,refUpdates] = simulate1fnoise(nsamp,4,0.01,0.85);


refUpdates(refUpdates == 0) = NaN;
figure, subplot(1,2,1), plot(t,sraw,t,s(1,:),'r'), 
        title('Signal with only additive noise, no updates'),
        xlabel('time')
        legend({'ideal measurement (short measurements)','physical signal'})
             subplot(1,2,2), loglog(t./nsamp,abs(fft(sraw))),
        xlim([0,0.5]),
        xlabel('frequency')

figure,subplot(1,2,1),plot(t,noise_update),       title('Update noise')
        hold on, stem(refUpdates.*noise_update),
        xlabel('time'), 
        subplot(1,2,2), loglog(t./nsamp,abs(fft(noise_update))),
        xlim([0,0.5]),
        xlabel('frequency')

figure,subplot(1,2,1)
        plot(t,sig_total',t,sraw,t,noise_update,'r'),       title('Measured signal'), xlabel('time')
        legend({'real measurement', 'ideal measurement (short measurements)','noise due to updates'})
        subplot(1,2,2), loglog(t./nsamp,abs(fft(sig_total))),xlim([0,0.5]),xlabel('frequency'), hold on
        loglog(t./nsamp,abs(fft(sraw))),xlim([0,0.5]),xlabel('frequency')



function [signal,acquisition_noUpdates,noise,acquisition,noise_update,refUpdates] = simulate1fnoise(nsamp,amp,freq,probabilityUpdate)
% nsamp -> number of samples
% amp   -> perturbation amplitude
% freq   -> frequency of perturbation
% probabilityUpdate -> starting probability of no update, decreases 5%
%                       every time it does not update
%
% Outputs: signal -> physical signal we want to measure
%          acquisition_noUpdates -> measurement in the case of no updates
%                                  (ideal case)
%          noise -> additive, unavoidable noise
%          acquisition -> data obtained from a typical measurement (with
%                         updates)
%          noise_update -> noise due to the updates (essentially, what we
%                          would want to find out from acqusition and refUpdates)
%          refUpdates   -> Vector signaling when there was an update.

nWindows =1;
t = 1:nsamp;
signal = amp*sin(2*pi*freq*t).*[zeros(1,250) hamming(nsamp-500)' zeros(1,250)] ;
noise = randn(nWindows,nsamp);
acquisition_noUpdates = signal + noise; % isto e o que se mediria sem actualizar a referencia

noise_update = zeros(nWindows,nsamp); % isto e o ruido resultante das actualizaçoes de referencia                           % (isto e o que queremos descobrir)
acquisition = zeros(nWindows,nsamp); % isto e o ruido que medimos com updates de referencia

refUpdates = zeros(nWindows,nsamp);

for wnd = 1:nWindows
    updateError = 0;
    countNoUpdates = -1;
    for idx = 1:nsamp
        countNoUpdates = countNoUpdates+1;
        acquisition(wnd,idx) = acquisition_noUpdates(wnd,idx) + updateError;
        noise_update(wnd,idx) = updateError;
        if rand > probabilityUpdate-(countNoUpdates*0.05)
            refUpdates(wnd,idx) = 1;
            updateError = updateError+noise(wnd,idx);
            countNoUpdates = 0;
        end
    end
end

end

EDIT: Just to add, some of the things I tried already include differentiating and median filtering (and then integrating the signal), which is not of much use in the presence of such low SNR. I tried this approach restricted to the update points as well, but the update density in the specific case is really high (this is intendend to be a worst case scenario. In real life, updates may be more infrequent.)

EDIT I've also been trying some things using more than one sensor. This meaning, having more than one redundant signal, in which the updates (step functions) are not necessarily at the same time. However, even in that scenario, I am unable to correct my signal. This is not ideal, as it implies sensor redundancy, but if it is the only possible solution is a good enough trade-off. Currently, I attempted to average the points at which there was an update with the average of the sensor measurements that didn't update, but for some reason it ends up being worse than simply averaging the signals in the real case, which is not good enough. I can post some simulation if anyone is interested. In that case (let's say 5 sensors, measuring the same signal, with independent updates), the result would be the following. Obviously, I can just average them. But is there a smarter way of going about it, since I also have the information from where the updates happen?

EDIT As kindly suggested by Peter K., I tried to use a recursive DC blocker. While it does improve the signal visually by attenuating the very low frequency components, it seems it also affects the signal amplitude quite heavily. Something I failed to mention in the post, is that we may be interested in even lower frequency signals, or higher frequency, with smaller amplitudes. As such, regular filtering by FIR or IIR which generally assume frequency information of the signal is not exactly what I was looking for. Instead, I was thinking of exploring the knowledge of the time-instants which contribute to 1/f noise, and potentially some redundacy in measurement to reduce in-band noise, in a more effective way than simply averaging.

These results were obtained using the filter presented in the following page, with a R constant of 0.95 : DC Blocker

Code to generate similar data: Code

EDIT Ok, I had this idea but maybe I need some help to develop it further. Piggybacking on Peter K.'s idea of using a DC blocker, I thought that maybe I could apply it only to the update points. This way, hopefully I would minimize the attenuation on the signal.

for p = 1:(nsamp-1)
    for wnd = 1:nWindows
        if refUpdates(wnd,p) == 1
                expectedval = corrected_measurement(wnd,p+1) - corrected_measurement(wnd,p) + R*corrected_measurement(wnd,p);              
            error = expectedval-corrected_measurement(wnd,p+1);

            corrected_measurement(wnd,p+1:end) = corrected_measurement(wnd,p+1:end)+error;

I am, at the same time, trying something maybe more adaptive, since this is a recursive filter. For instance, if I have a series of successive updates every sample (worst case scenario), it's probably a bad idea to use a DC blocker strongly, since it would eat the signal up. So I thought maybe R could depend on the number of samples it has been since the last step. (Let's say, R starts at 1 if the last step was on the previous sample, which then decreases 0.01 per sample).

Now, it does seem to improve a bit, but it does affect the signal. However, in the case of 10 redundant sensors, when averaging there seems to be considerable improvement. On that note, one other thing I tried, was to use the term to be multiplied by R in the DC blocker filter the average of the measurement by other sensors, which were not being updated at the time (in a way, try to have each sensor correct each other). I am not sure how good my implementation is, or if any of this assumptions are not smart at all, but i got the following results.

Keep in mind, these came from a bit more lengthy simulations, for a more realistic approach (1khz sampling rate, 10second acquisition, 1 Hz perturbation)

Up to now, this idea of using other sensors redundant information to correct the "update" step-like noise seems like the most promising solution to avoid damaging the low frequency response. I wonder if there are any additional thoughts.

EDIT

Okay, the thing I am looking into right now is to improve the measurement using a 10 sensor redundancy. Taking the average measurement at the non-updated sensors, and the mean of a small section of points of the updated sensor following the update should enable some sort of correction. The way I see it, assuming that the estimation of the mean level of the step is perfect (it won't be), then having the average of 9 good sensor measurements would yield a step error with an amplitude $\sqrt{9}=3$ times smaller. Averaging shouldn't be a big problem since we are very oversampled(1kHz sampling, interested in perturbations up to 20hz), although there will be some signal dependent bias on the compensation. I wonder if there would be other downsides to such an approach. When simulated, it seems to attain a very small improvement in some cases, and no improvement in others, so I am not convinced by the method at all.

Answer 1

OK, this is still not perfect, but what I'd suggest is first filtering using a notch filter:

function [ y ] = notch_filter(alpha,  omega, x )
%notch_filter Notch filter

b = [1 -2*cos(omega) 1];
a = [1 -2*alpha*cos(omega) alpha*alpha];

y = filter(b, a, x);

end

followed by a DC blocker:

function [ y ] = dc_blocker( alpha, x )
%dc_blocker BLock DC using coefficient alpha.

y = filter([1 -alpha], [1 -1], x);

end

Like this:

processed = notch_filter(alpha_notch, fc, acquisition_signal(1,:));
processed = dc_blocker(alpha_dc, processed);

where, for the single example below, I've selected alpha_notch = alpha_dc = 0.999;

Here, I've added the processed plot in the green overlaying your pure_signal in the top plot. This example was one of the better ones; there were others that are not great. I'll think about the problem a bit more and see if anything else comes to mind.

Answer 2

I do not have access to matlab, but I suggest a sliding window with length of several periods of the lowest frequency you are interested in. Of this window, calculate the average and subtract that value from center sample. When I have some time, I can try it in python.