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

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')

Answer 1

[William Rose]

@종영,

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.