How do I recover continous signal from descrete measurements?

3 vues (au cours des 30 derniers jours)
Brian
Brian le 8 Juil 2021
Modifié(e) : Brian le 12 Juil 2021
I have a slowly varying signal that is always increasing or decreasing in value. I measure this with a descrete device (both signal and time value). I would like to determine the gradient of the true signal and therefore I would like to recover/estimate the true signal first. I have made a simple example which have similar behavior as what I am dealing with. Assuming that I know the signal values my measureing device can give (e.g. integer muliplum of 0.1) can I then recover the real signal with high accuracy in the gradient?
I appreciate any hint help you can give me to go in the right direction.
x = 1:0.01:10;
y = 1./x;
figure
plot(x,y)
hold on
y_discrete = round(y*10)/10;
plot(x,y_discrete,'.')
ylim([0,1])
EDIT:
Based on comments I have added a little extra description and a suggestion for a solution.
Assumptions:
  • There is no noise (it will be removed before this step)
  • Within each constant value block (the red lines in the plot above) the true value change slowly
  • The estimation of the true signal is done as postprocessing so I have the whole dataset available and the calculation time is not an issue.
%Real and discrete sigmal
x = 1:0.01:10;
y = 1./x;
figure
plot(x,y)
hold on
y_discrete = round(y*10)/10;
plot(x,y_discrete,'.')
ylim([0,1])
%Approximation of the real signal using splines assuming that when the
%value of the discrete signal changes a good approximation of the real value
%is the average between the two levels.
idx = diff(y_discrete)';
idx = find(idx ~= 0);
x_intrp = zeros(length(idx)+2,1);
y_intrp = zeros(length(idx)+2,1);
x_intrp(1) = x(1);
y_intrp(1) = y(1);
for i = 1:length(idx)
y_intrp(i+1) = (y_discrete(idx(i))+y_discrete(idx(i)+1))/2;
x_intrp(i+1) = (x(idx(i))+x(idx(i)+1))/2;
end
x_intrp(i+2) = x(end);
y_intrp(i+2) = y(end);
plot(x_intrp,y_intrp,'r*')
m = makima(x_intrp,y_intrp,x); %Akima piecewise cubic Hermite interpolation
plot(x,m)
legend('True signal', 'Discrete signal', 'True nodes estimates', 'Spline estimation')
  3 commentaires
Bjorn Gustavsson
Bjorn Gustavsson le 8 Juil 2021
If you know for certain the shape of the underlying function you might have a fighting chance by looking at some general least-square fitting between a parameterization of that function, taking into account the discretization of your measurements etc. But there will be problems with noise-sensitivity of the entire process, so you will have to temper your expectations.
Brian
Brian le 12 Juil 2021
Hi @Bjorn Gustavsson. I have added some more information about the problem and a solution. For this specific case the solution works fine. However, it might not work well for other cases and maybe there exists more robust methods cabable of something similar (Estimating true signal from discrete measurements).

Connectez-vous pour commenter.

Réponse acceptée

Bjorn Gustavsson
Bjorn Gustavsson le 12 Juil 2021
Maybe you get something thats a bit more general (and perhaps?) more robust by simple least-square fitting:
m_fun = @(pars,x) round(interp1(pars(1:2:end),pars(2:2:end),x,'pchip')*10)/10; % 2 slightly different
M_fun = @(pars,x) round(interp1(pars(1:2:end),pars(2:2:end),x,'spline')*10)/10;% rounding model-functions
% you might consider adding arbitrary bells-and-whistles to this...
% Error-function:
err_fcn = @(pars,x,y,fcn) sum((y-fcn(pars,x)).^2);
par0 = [1 1 3 1/3 6 1/6 9 1/9]; % initial guess for 4 node-points for the approximating spline
par1 = fminsearch(@(pars) err_fcn(pars,x,y,m_fun),par0);
Par1 = fminsearch(@(pars) err_fcn(pars,x,y,M_fun),par0);
% disp(Par1) returns:
% Par1 =
% 0.98183 0.82362 2.9534 0.35617 6.6921 0.14992 8.8114 0.12509
% disp(par1) returns:
% par1 =
% 1.0057 0.89591 2.8314 0.33493 6.0081 0.16794 9.7301 0.10913
clf
plot(x,y)
hold on
plot(x,y_discrete,'.')
plot(x,m_fun(par1,x),'b')
plot(x,M_fun(Par1,x),'g')
In any non-artificial case you'll at least have to consider discretization-noise at each jump in the data - regardless of what noise-removal technique you consider...
HTH

Plus de réponses (0)

Produits


Version

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by