MLS(moving least squares) algorithm problem... help me!!!

23 次查看（过去 30 天）
종영
종영 2024-11-7，11:34
编辑： William Rose 2024-11-7，20:55
Ran in:
During MLS, it was conducted as follows. I'm curious about this, I designated weighting for five adjacent points for one point, but I don't know how I can express it diagonally for 100 points as shown in the picture.
I don't know if the code just expressed this as sum or not.
The result value also doesn't come up with an MLS... Help!
clc; clear all; close all
%% give the nodes
x = linspace(0, 1, 100);
y = sin(2 * pi * x) + 0.5 * cos(6 * pi * x + pi / 4) ;
% y_noisy = y;
y_noisy = awgn(y,20);
figure(1)
plot(x,y)
hold on
plot(x,y_noisy,'.r')
%% basis function
ps = [ones(1,length(x)) ;x ;y;x.^2 ; (x.*y);y.^2].';
%% cubic spline function
distances = NaN(length(x), 5);
% weighting = zeros(size(relative_d));
for i = 1 : length(x)
point_index = i;
xi = x(point_index);
yi = y_noisy(point_index);
if point_index == 1
d = sqrt((x(point_index+1:point_index+5) - xi).^2 + (y_noisy(point_index+1:point_index+5) - yi).^2);
distances(i, 1:length(d)) = d;
elseif point_index == 2
distances1 = sqrt((x(point_index-1) - xi).^2 + (y_noisy(point_index-1) - yi).^2);
distances2 = sqrt((x(point_index+1:point_index+4) - xi).^2 + (y_noisy(point_index+1:point_index+4) - yi).^2);
d = [distances1, distances2];
distances(i, 1:length(d)) = d;
elseif point_index > 2 && point_index < length(x) - 2
distances1 = sqrt((x(point_index-2:point_index-1) - xi).^2 + (y_noisy(point_index-2:point_index-1) - yi).^2);
distances2 = sqrt((x(point_index+1:point_index+3) - xi).^2 + (y_noisy(point_index+1:point_index+3) - yi).^2);
d = [distances1, distances2];
distances(i, 1:length(d)) = d;
elseif point_index == length(x) - 2
distances1 = sqrt((x(point_index-3:point_index-1) - xi).^2 + (y_noisy(point_index-3:point_index-1) - yi).^2);
distances2 = sqrt((x(point_index+1:point_index+2) - xi).^2 + (y_noisy(point_index+1:point_index+2) - yi).^2);
d = [distances1, distances2];
distances(i, 1:length(d)) = d;
elseif point_index == length(x) - 1
distances1 = sqrt((x(point_index-4:point_index-1) - xi).^2 + (y_noisy(point_index-4:point_index-1) - yi).^2);
distances2 = sqrt((x(point_index+1) - xi).^2 + (y_noisy(point_index+1) - yi).^2);
d = [distances1, distances2];
distances(i, 1:length(d)) = d;
elseif point_index == length(x)
d = sqrt((x(point_index-5:point_index-1) - xi).^2 + (y_noisy(point_index-5:point_index-1) - yi).^2);
distances(i, 1:length(d)) = d;
end
dmi = mean(distances, 2); %질문!
% dmi(1:100) = 10;
relative_d(i,:) = abs(distances(i,:))./dmi(i);
for j = 1:5
if relative_d(i,j) < 0.5
weighting(i,j) = 2/3 - 4*relative_d(i,j).^2 + 4*relative_d(i,j).^3;
elseif relative_d(i,j) > 0.5 && relative_d(i,j) < 1
weighting(i,j) = 4/3 - 4*relative_d(i,j) + 4*relative_d(i,j).^2 - (4/3)*relative_d(i,j).^3;
else
weighting(i,j) = 0;
end
end
weight_mean = sum(weighting, 2) /5; % 질문!
weight_fin = diag(weight_mean);
p_x(:,i) = [1; xi; yi; xi.^2; (xi .* yi); yi.^2];
end
A = ps.' * weight_fin * ps;
B = ps.' * weight_fin;
phi= p_x.' * (A^-1 * B) * y_noisy.';
% for i = 1 : 100
% figure(2)
% plot(1:5,weighting(i,:))
% hold on
% end
%% Obtain the matrixed
%
figure(2)
plot(x,y)
hold on
plot(x,y_noisy)
plot(x,phi)
legend('orginal','noise signal','MLS signal')

请先登录，再回答此问题。

回答（1 个）

William Rose
William Rose 2024-11-7，20:42
编辑：William Rose 2024-11-7，20:55
Ran in:
@종영,
In your comment you refer to "cubic spline function". Cubic splines for smoothing can be done with csaps(). The code below shows spline smoothing with three levels of smoothing: no smoothing (p=1); intermediate smoothing (p=.9999948); max smoothing (p=0).
n=100; % number of data points
x = linspace(0, 1, n);
y = sin(2 * pi * x) + 0.5 * cos(6 * pi * x + pi / 4) ;
yNoisy = awgn(y,10); % add Gaussian white noise, snr=10 dB
figure
subplot(211), plot(x,y,'-k',x,yNoisy,'.k')
legend('No noise','With noise'); grid on
Fit cubic splines to the data, with varying amounts of smoothing:
h=x(2)-x(1); % h=delta x
% p=1//(1+h^3/6); % initial guess for smoothing parameter p
% make p smaller/larger for more/less smoothing
ys1=csaps(x,yNoisy,0,x); % p=0: straight line
p=1/(1+h^3/0.2); % smoothing parameter p
ys2=csaps(x,yNoisy,p,x); % p=smoothing paramater
ys3=csaps(x,yNoisy,1,x); % p=1: cubic spline thorugh every data point
Plot results
subplot(212), plot(x,yNoisy,'.k',x,ys1,'-b',x,ys2,'-g',x,ys3,'-r')
legend('Data','p=0',['p=',num2str(p,'%.7f')],'p=1'); grid on
I realize this does not address why your code doesn't work. But the solution above does get the job done, and it is probably easier than debugging your code.

类别

在 Help Center 和 File Exchange 中查找有关 Smoothing 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by