在噪声信号中找到模型函数的鲁棒拟合

问题描述

我有一个嘈杂的信号和一个模型函数，例如：

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？还是具有不同的鲁棒性？

解决方法

这只是实现@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

,

我修改了测试较硬L1假设的成本函数，并且与David的答案相比拟合度更高（我测试的噪声元素数量更多，请参见 bar baz foo ab c d efgh ab c d efgh ...... A 6.583721 -1.554734 1.922187 1.100208 4.944954 -1.343831 0.939265 -3.614612 ...... B 6.138441 0.653721 -0.204472 1.890755 -0.347505 1.633708 0.392096 0.414880 ...... C 0.951489 2.695940 -1.494028 0.907464 1.905409 -1.021097 -2.399670 0.799798 ...... 的添加内容）：

v3

,

MATLAB中的ransac实现似乎只能在预定义的拟合函数中正常工作。

但是，您可以通过将拟合函数包装在函数中来创建解决方法。该功能应该是*“计算机视觉工具箱”或自动驾驶系统工具箱的一部分。

x = linspace(0,20).';
w = [2 6 -4 5];
y = w(1)*besselj(0,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);

% 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(:,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)开始。您还可以在此处集成其他优化算法。