在噪声信号中找到模型函数的稳健拟合

Question

我有一个噪声信号和一个模型函数，例如：

x=linspace(0,20);
w=[2 6 -4 5];
y=w(1)*besselj(0,x)+w(2)*besselj(1,x)+w(3)*besselj(2,x)+w(4)*besselj(3,x);
y(randi(length(y),[1 10]))=10*rand(1,10)-5;
plot(x,y,'x')

我想使用 RANSAC 在我的模型中查找w，因为这种方法在查找线时对噪声具有鲁棒性。然而，这不是一个线性问题，我无法得到合适的拟合，可能是因为我试图拟合的函数的振荡性质。

我看到 matlab 有一个 fitPolynomialRansac 函数，但是对于 a+b*x+c*x^2+d*x^3 简单案例（介于 -1 和 1 之间），即使这样也失败了。

知道如何驯服 RANSAC 吗？还是采用不同的抗噪方法？

Answer 1

MATLAB 中的ransac-implementation 似乎只能与预定义的拟合函数一起使用。

但是，您可以通过将拟合函数包装在函数中来创建解决方法。该函数应该是 *Computer Vision Toolbox" 或 Automated Driving System Toolbox 的一部分。

x = linspace(0,20).';
w = [2 6 -4 5];
y = w(1)*besselj(0,x)+w(2)*besselj(1,x)+w(3)*besselj(2,x)+w(4)*besselj(3,x);
y(randi(length(y),[1 10])) = 10*rand(1,10)-5;
plot(x,y,'x')



sampleSize = 5; % number of points to sample per trial
maxDistance = 2; % max allowable distance for inliers

% creating the function you want to obtain
bFnc = @(x,w) w(1)*besselj(0,x)+w(2)*besselj(1,x)+w(3)*besselj(2,x)+w(4)*besselj(3,x);

% wrapping this function into a cost function with two inputs (the points
% and the "model"
cstFmin = @(xy,w) sum((xy(:,2) - bFnc(xy(:,1),w)).^2);

% creating a function handle that fits the function + returns a single
% object as a model
fitFnc = @(xy) {fminsearch(@(w)cstFmin(xy,w),ones(4,1))}; % returns the "model" as a cell

% build a function that determines a distance measure
evlFnc = @(model,xy)(( xy(:,2) - bFnc(xy(:,1),model{1}) ).^2);

% call RANSAC
[modelRANSAC, inlierIdx] = ransac([x,y],fitFnc,evlFnc, sampleSize,maxDistance);

% plot results
% plot
xIN = x(inlierIdx);
yIN = bFnc(xIN,modelRANSAC);

hold on
plot(xIN,yIN,'r-')
hold off
title(num2str(sampleSize,'sampleSize = %d'))

请注意，fminsearch-算法总是从ones(4,1).. 开始。您还可以在此处集成不同的优化算法。

Answer 2

我已经修改了成本函数来测试一个硬（呃）L1 假设，并且与 David 的答案相比，得到了更稳健的拟合（我正在使用更高的噪声元素进行测试，请参阅我对 v3 的补充）：

x = linspace(0,20);

% model function
yFun=@(w) w(1)*besselj(0,x)+w(2)*besselj(1,x)+w(3)*besselj(2,x)+w(4)*besselj(3,x);

% generate training data
N = 50; % number of noisy elements
w=[2 6 -4 5]; % true paramater values
y = yFun(w); % evaluate true model function
y(randi(length(y),[1 N])) = 10*rand(1,N)-5; % add noise

% standard loss function for least sqaure
d = @(w) yFun(w)-y;

v1 = lsqnonlin(d,[1 1 1 1]); % normal least squares
v2 = lsqnonlin(@(w) sqrt(2*(sqrt(1+d(w).^2)-1)),[1 1 1 1]) % soft L1 loss
v3 = lsqnonlin(@(w) sqrt(2*(sqrt(1+abs(d(w)))-1)) ,[1 1 1 1]) %  the new cost function


plot(x,y,'x',x,yFun(v1),x,yFun(v2),x,yFun(v3))
legend('data','v1','v2','v3')

Answer 3

这只是实现@mikuszefski 的评论以使用软 L1 损失函数。它似乎确实更耐噪音：

x = linspace(0,20);

% model function
yFun=@(w) w(1)*besselj(0,x)+w(2)*besselj(1,x)+w(3)*besselj(2,x)+w(4)*besselj(3,x);

% generate training data
N = 20; % number of noisy elements
w=[2 6 -4 5]; % true paramater values
y = yFun(w); % evaluate true model function
y(randi(length(y),[1 N])) = 10*rand(1,N)-5; % add noise

% standard loss function for least sqaure
d = @(w) yFun(w)-y;

v1 = lsqnonlin(d,[1 1 1 1]); % normal least squares
v2 = lsqnonlin(@(w) sqrt(2*(sqrt(1+d(w).^2)-1)),[1 1 1 1]) % soft L1 loss