How to denoise a response to step forcing?

The characteristic feature of a signal that is a response to a step forcing is that it changes fast in the beginning and then it "decays to" an asymptote, i.e., it changes slower an slower. This seems to cause a problem with denoising. The problem is compounded if the corrupting noise is not white but have some persistence, say, it is an autoregressive process. I have tried a Savitzky-Golay filter (sgolay) and a Wavelet denosing algorithm. They perform similarly when the window size for sgolay is not long. But this is not a satisfactory result, as the filtered signals are closer to the corrupted signals than the uncorrupted sought-for one. This is clearly mostly the consequence of the persistence of the noise. I can extend the window size for sgolay in order to achieve a better result on the plateau, but the performance is deterriorating during the early fast-change period. The former is demonstarted by the code snippet by showing two realisations of the corrupted signal such that the smoothed curves fit closer to each other and to the true signal. The latter is demonstarted by either of the smoothed curves diverging considerably from the true signal.

In fact, we know that at zero time the signal value is zero. Is there a way to constrain the smoothing to respect this one data point, i.e., that the smoothed curve is anchored to it?

Also, given that we have two realisations of the corrupted signal, we can actually estimate the noise strength (and even the other autoregressive parameter). Is there a technique by which this information can be utilised to achieve a better smoothing?

Note: I didn't use 'impulseest' as it appears to me that it considers the noise corruption as part of the signal.

t = 0:2^8; 
N = length(t);

tau = 10; % time scale of relaxation 
x_asym = 4; % relative increase of the observable x t -> infinity

x_r2s = x_asym*(1 - exp(-t/tau)); % uncorrupted sought-for signal

sig = 0.05; % noise strength 
% Use AR(1) to corrupt the signal instead of white noise.
phi = 0.6; % AR coefficeint 

xi1 = sig*randn(1,N);
xi2 = sig*randn(1,N);
pert1 = filter(1,[1 -phi],xi1);
pert2 = filter(1,[1 -phi],xi2);

x_prd_1 = x_r2s+pert1;
x_prd_2 = x_r2s+pert2;

sys = ar(x_prd_1-x_prd_2,1);

sig_est = 2*sys.NoiseVariance; 
phi_est = -sys.A(2); 

x_sg_11 = sgolayfilt(x_prd_1,3,11);
x_sg_21 = sgolayfilt(x_prd_2,3,11);

x_sg_12 = sgolayfilt(x_prd_1,3,101);
x_sg_22 = sgolayfilt(x_prd_2,3,101);

x_wd_1  =   wdenoise(x_prd_1);
x_wd_2  =   wdenoise(x_prd_2);

figure
%semilogx(t,x_r2s,t,x_prd_1,t,x_sg_11,t,x_prd_2,t,x_sg_21)
plot(t,x_r2s,t,x_prd_1,t,x_sg_11,t,x_prd_2,t,x_sg_21)
figure
%semilogx(t,x_r2s,t,x_prd_1,t,x_sg_12,t,x_prd_2,t,x_sg_22)
plot(t,x_r2s,t,x_prd_1,t,x_sg_12,t,x_prd_2,t,x_sg_22)
figure
%semilogx(t,x_r2s,t,x_prd_1,t,x_wd_1 ,t,x_prd_2,t,x_wd_2 )
plot(t,x_r2s,t,x_prd_1,t,x_wd_1 ,t,x_prd_2,t,x_wd_2 )

Answer 1

[Mathieu NOE]

hello

I tried 2 simple methods :

• sliding window averaging : not very effective
• to fit an exp model : will almost give the exact original signal

second option seems to be the best (simple) solution . i did not test the rest of the code as I miss some Tbx

t = 0:2^8; 
N = length(t);
tau = 10; % time scale of relaxation 
x_asym = 4; % relative increase of the observable x t -> infinity
x_r2s = x_asym*(1 - exp(-t/tau)); % uncorrupted sought-for signal
sig = 0.05; % noise strength 
% Use AR(1) to corrupt the signal instead of white noise.
phi = 0.6; % AR coefficeint 
xi1 = sig*randn(1,N);
xi2 = sig*randn(1,N);
pert1 = filter(1,[1 -phi],xi1);
pert2 = filter(1,[1 -phi],xi2);
x_prd_1 = x_r2s+pert1;
x_prd_2 = x_r2s+pert2;


% denoising using sliding avg method 
N = 10;
x_prd_1savg = myslidingavg(x_prd_1, N);

% denoising using exponential fit method 
% code is giving good results with template equation : % y = a.*(1-exp(b.*(x-c)));
x = t;
y = x_prd_1;

f = @(a,b,c,x) a.*(1-exp(b.*(x-c)));
obj_fun = @(params) norm(f(params(1), params(2), params(3),x)-y);
sol = fminsearch(obj_fun, [y(end),0,0]);
a_sol = sol(1);
b_sol = sol(2);
c_sol = sol(3);
x_prd_1fit = f(a_sol, b_sol,c_sol, x);

figure
plot(t,x_r2s,'b',t,x_prd_1,'r',t,x_prd_1savg,'c',t,x_prd_1fit,'k');
legend('signal','signal+noise','after sliding avg','exp fit');

function out = myslidingavg(in, N)
%   OUTPUT_ARRAY = MYSLIDINGAVG(INPUT_ARRAY, N)
%
%  The function 'slidingavg' implements a one-dimensional filtering, applying a sliding window to a sequence. Such filtering replaces the center value in
%  the window with the average value of all the points within the window. When the sliding window is exceeding the lower or upper boundaries of the input
%  vector INPUT_ARRAY, the average is computed among the available points. Indicating with nx the length of the the input sequence, we note that for values
%  of N larger or equal to 2*(nx - 1), each value of the output data array are identical and equal to mean(in).
%
%  *  The input argument INPUT_ARRAY is the numerical data array to be processed.
%  *  The input argument N  is the number of neighboring data points to average over for each point of IN.
%
%  *  The output argument OUTPUT_ARRAY is the output data array.

if (isempty(in)) | (N<=0)                                              % If the input array is empty or N is non-positive,
 disp(sprintf('SlidingAvg: (Error) empty input data or N null.'));     % an error is reported to the standard output and the
 return;                                                               % execution of the routine is stopped.
end % if

if (N==1)                                                              % If the number of neighbouring points over which the sliding 
 out = in;                                                             % average will be performed is '1', then no average actually occur and
 return;                                                               % OUTPUT_ARRAY will be the copy of INPUT_ARRAY and the execution of the routine
end % if                                                               % is stopped.

nx   = length(in);             % The length of the input data structure is acquired to later evaluate the 'mean' over the appropriate boundaries.

if (N>=(2*(nx-1)))                % If the number of neighbouring points over which the sliding 
 out = mean(in)*ones(size(in));   % average will be performed is large enough, then the average actually covers all the points
 return;                          % of INPUT_ARRAY, for each index of OUTPUT_ARRAY and some CPU time can be gained by such an approach.
end % if                          % The execution of the routine is stopped.

out = zeros(size(in));         % In all the other situations, the initialization of the output data structure is performed.

if rem(N,2)~=1                 % When N is even, then we proceed in taking the half of it:
 m = N/2;                      % m = N     /  2.
else                           % Otherwise (N >= 3, N odd), N-1 is even ( N-1 >= 2) and we proceed taking the half of it:
 m = (N-1)/2;                  % m = (N-1) /  2.
end % if

for i=1:nx,                 % For each element (i-th) contained in the input numerical array, a check must be performed:
  dist2start = i-1;         %  index distance from current index to start index (1)
  dist2end = nx-i;          %  index distance from current index to end index (nx)
  if dist2start<m || dist2end<m        % if we are close to start / end of data, reduce the mean calculation on centered data vector reduced to available samples
      dd = min(dist2start,dist2end);    % min of the two distance (start or end)
  else
      dd = m;
  end % if
  out(i) = mean(in(i-dd:i+dd));     % mean of centered data , reduced to available samples at both ends of the data vector

end % for i

end

[Tamas Bodai]

Hi Mathieu. Thanks so much for your answer! I see that you use your own sliding window averaging function instead of Matlab's 'movmean'. I run a quick test verifying that they don't do exactly the same when the window reaches beyond the data span. Regarding the other solution, which performs well on the example, unfortunately this is not what would work for me. I used the simple exponential signal form just as a test example. In the actual data that i would like to analyse we don't know the functional form of the signal. It is like an exponential, or rather a sum of two exponentials with fairly large time scale separation (for which reason it doesn't have a plateau but keeps increasing in the observation period), but it also has some further features, some wiggles on shorter time scales, perhaps somewhat longer than the characteristic time scale of the corrupting noise process. I would like to resolve as much unknown features of the signal as possible. sgolay seems like a good approach, but lets me down at the fast changing initial period and with the anchoring. Perhaps a variable window size version of sgolay with the possibility of prescribing some constraints would be best, but i'm not sure that it implies a well-defined algorithm, and i certainly cannot easily find a function like this on the internet. Or, perhaps the wavelet denoising can also be adapted using the knowledge about the noise that the two realisations provide us.

[Mathieu NOE]

hello Tamas

if your signal could be modelized as a combination of exponentials , why not trying to fit a more representative model ?

on the other side I was also thinking about an adaptative window algorithm, but that requires also some a priori knowledge or hypothesis.

do you have "true" measured data to test ?

[Tamas Bodai]

Hi Mathieu, thanks a lot for your reply. The answers to your questions:

As i say, it is not exactly a sum of exponentials, and i would like to recover as much of the true signal as possible.

Yes, i can share my measured data; please see it attached. That is, i have two realisations of the noise corrupted signal. The underlying corrupted signal is the same, without any doubt, only the realisations of the noise process are different. As for the noise, the generating process is also identical without any doubt, and the relaisations are completely independent.

I would really appreciate any further help.

[Mathieu NOE]

I let you test the code below and we 'll discuss it later on

load('x_prd_12.mat')
% data = x_prd_1;
data = x_prd_2;

samples = length(data);
data = double(data);

x = 1:samples;

x_log = logspace(0,log10(samples),samples);
data_int = interp1(x,data,x_log,'linear');

N = 40;
data_avg = movmean(data_int,N);
data_avg2 = interp1(x_log,data_avg,x,'linear');
data_avg2(1) = data(1);

figure(1),plot(x,data,'b',x,data_avg2,'r');

[Tamas Bodai]

Thanks, Mathieu, this looks rather promising! The same idea of resampling x should also work with the sgolay filter. I will try that later, and also apply to some other observed response signals that i have. Really appreciate the help!

[Mathieu NOE]

you're welcome !

[Mathieu NOE]

hello Tamas

tried now with s-golay filter - good results with low order filter

load('x_prd_12.mat')
% data = x_prd_1;
data = x_prd_2;

samples = length(data);
data = double(data);

x = 1:samples;

x_log = logspace(0,log10(samples),samples);
data_int = interp1(x,data,x_log,'pchip');

% N = 80;
% data_avg = movmean(data_int,N);
% data_avg2 = interp1(x_log,data_avg,x,'spline');
% data_avg2(1) = data(1);

data_avg = sgolayfilt(data_int,1,81); 
data_avg2 = interp1(x_log,data_avg,x,'spline');
data_avg2(1) = data(1);    

figure(1),plot(x,data,'b',x,data_avg2,'r');

[Tamas Bodai]

Thanks a lot, Mathieu! In addition to your resampling trick, we need to optimise the 1. window size and the 2. order of the sgolay filter. The first one can be done by 1. predicting the noise variance from parting the ensemble into two, and picking the window size with which the variance of residuals for the complete ensemble matches the predicted noise variance the closest. The latter can be done by 2. minimizing the difference of the smoothed curves for the two subensembles.

[Mathieu NOE]

ok , glad to see that this is alive !