大唐——直轴

采取鲁棒性较强的二阶傅里叶拟合，对数据进行拟合。拟合代码如下

%%
x1=Sheet.(2);
y1=Sheet.(1);

x2=Sheet.(5);
y2=Sheet.(4);

x3=Sheet.(8);
y3=Sheet.(7);
% %将上面的语句改为循环
% x=zeros(333,3);
% for i=2:3:8
%     x()
% %拟合方式的选择
% cftool(x1,y1);

%%
% 第一个点
[a_1,b_1]=createFits(x1, y1,x2,y2);
a_1{5,1}.a0
FUNCTION=a_1{5,1}.a0 + a_1{5,1}.a1*cos(x1*a_1{5,1}.w)+ a_1{5,1}.b1*sin(x1*a_1{5,1}.w) +a_1{5,1}.a2*cos(2*x1*a_1{5,1}.w) + a_1{5,1}.b2*sin(2*x1*a_1{5,1}.w);
[pks_po,locs_po]=findpeaks(FUNCTION);
[pks_pa,locs_pa]=findpeaks(-FUNCTION);
plot(FUNCTION);
hold on
plot(locs_po,pks_po,'o');
[pks1,locs1]=findpeaks(-FUNCTION);
plot(locs_pa,-pks1,'*');
hold on;
if length(pks_po)>1
    locs_po=locs_po(1);
    pks_po=pks_po(1);
end
if length(pks_pa)>1
    locs_pa=locs_pa(1);
    pks_pa=pks_pa(1);
end
if locs_po>100
    a=locs_po;
    b=locs_pa;
    c=pks_po;
    d=pks_pa;
    locs_po=b;
    locs_pa=a;
    pks_pa=d;
    pks_po=c;
end
disp('第一点的相角与极大值')
disp(locs_po),disp(pks_po)
disp('第一点的相角与极小值')
disp(locs_pa),disp(-pks_pa)

%第二个点
%第二组数据很平，这时候就随便取一个当作最小跳动量
[a_2,b_2]=createFits(x2, y2,x3,y3);
a_2{5,1}.a0;
FUNCTION_2=a_2{5,1}.a0 + a_2{5,1}.a1*cos(x2*a_2{5,1}.w)+ a_2{5,1}.b1*sin(x2*a_2{5,1}.w) +a_2{5,1}.a2*cos(2*x2*a_2{5,1}.w) + a_2{5,1}.b2*sin(2*x2*a_2{5,1}.w);
[pks_po_2,locs_po_2]=findpeaks(FUNCTION_2);
[pks_pa_2,locs_pa_2]=findpeaks(-FUNCTION_2);
plot(FUNCTION_2);
hold on;
plot(locs_po_2,pks_po_2,'o');
[pks_2,locs_2]=findpeaks(-FUNCTION_2);
plot(locs_pa_2,-pks_2,'*');
hold on;
if length(pks_po_2)>1
    locs_po_2=locs_po_2(1);
    pks_po_2=pks_po_2(1);
end
if length(pks_pa_2)>1
    locs_pa_2=locs_pa_2(1);
    pks_pa_2=pks_pa_2(1);
end
if locs_po_2>100
    a=locs_po_2;
    b=locs_pa_2;
    c=pks_po_2;
    d=pks_pa_2;
    locs_po_2=b;
    locs_pa_2=a;
    pks_pa_2=d;
    pks_po_2=c;
end
disp('第二点的相角与极大值')
disp(locs_po_2),disp(pks_po_2)
disp('第二点的相角与极小值')
disp(locs_pa_2),disp(-pks_pa_2)

%第三个点
[a_3,b_3]=createFits(x3, y3,x1,y1);
a_3{5,1}.a0;
FUNCTION_2=a_3{5,1}.a0 + a_3{5,1}.a1*cos(x3*a_3{5,1}.w)+ a_3{5,1}.b1*sin(x3*a_3{5,1}.w) +a_3{5,1}.a2*cos(2*x3*a_3{5,1}.w) + a_3{5,1}.b2*sin(2*x3*a_3{5,1}.w);
[pks_po_3,locs_po_3]=findpeaks(FUNCTION_2);
[pks_pa_3,locs_pa_3]=findpeaks(-FUNCTION_2);
plot(FUNCTION_2);
hold on
plot(locs_po_3,pks_po_3,'o');
[pks_2,locs_2]=findpeaks(-FUNCTION_2);
plot(locs_pa_3,-pks_2,'*');
hold on;
if length(pks_po_3)>1
    locs_po_3=locs_po_3(1);
    pks_po_3=pks_po_3(1);
end
if length(pks_pa_3)>1
    locs_pa_3=locs_pa_3(1);
    pks_pa_3=pks_pa_3(1);
end
if locs_po_3>100
    a=locs_po_3;
    b=locs_pa_3;
    c=pks_po_3;
    d=pks_pa_3;
    locs_po_3=b;
    locs_pa_3=a;
    pks_pa_3=d;
    pks_po_3=c;
end
disp('第三点的相角与极大值')
disp(locs_po_3),disp(pks_po_3)
disp('第三点的相角与极小值')
disp(locs_pa_3),disp(-pks_pa_3)

%% %至此，上面第一步已经完成，输出了每个点的极值以及相位；下面开始第二部，进行曲面拟合
position=[locs_pa;locs_pa_2;locs_pa_3];
% sub1= locs_pa_2-locs_pa
% sub2= locs_pa_3-locs_pa_2
disp('下压点为2点');
disp('旋转的角度')
theta=(locs_pa_2+locs_pa_3+locs_pa)/3
pks_pa=pks_pa*cos(theta-270);
pks_pa_2=pks_pa_2*cos(theta-270);
pks_pa_3=pks_pa_3*cos(theta-270);

%% %第三步，输出其下压量
disp('下压量');
pks_pa_2-(pks_pa+pks_pa_3)/2
%%
% function [fitresult, gof] = createFit(x1, y1)
% [xData, yData] = prepareCurveData( x1, y1 );
% ft = fittype( 'fourier2' );
% opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
% opts.Display = 'Off';
% opts.StartPoint = [0 0 0 0 0 0.017459112224018];
% % Fit model to data.
% [fitresult, gof] = fit( xData, yData, ft, opts );
% % Plot fit with data.
% figure( 'Name', 'untitled fit 1' );
% h = plot( fitresult, xData, yData );
% legend( h, 'y vs. x', 'fit curve', 'Location', 'NorthEast', 'Interpreter', 'none' );
% xlabel( 'x', 'Interpreter', 'none' );
% ylabel( 'y', 'Interpreter', 'none' );
% grid on
% end

%%
function [fitresult, gof] = createFits(x1, y1, x2, y2)
% 初始化.
% 初始化数组以存储拟合度和拟合优度。
fitresult = cell( 5, 1 );
gof = struct( 'sse', cell( 5, 1 ), ...
    'rsquare', [], 'dfe', [], 'adjrsquare', [], 'rmse', [] );

[xData, yData] = prepareCurveData( x1, y1 );
%设置拟合类型和精度.
ft = fittype( 'poly1' );
% 使模型符合数据.
[fitresult{1}, gof(1)] = fit( xData, yData, ft );
% % 画图.
% figure( 'Name', 'untitled fit 1' );
% h = plot( fitresult{1}, xData, yData );
% legend( h, 'y1 vs. x1', 'untitled fit 1', 'Location', 'NorthEast', 'Interpreter', 'none' );
% % Label axes
% xlabel( 'x1', 'Interpreter', 'none' );
% ylabel( 'y1', 'Interpreter', 'none' );
% grid on

[xData, yData] = prepareCurveData( x2, y2 );
ft = fittype( 'poly1' );
[fitresult{2}, gof(2)] = fit( xData, yData, ft );
% figure( 'Name', 'untitled fit 2' );
% h = plot( fitresult{2}, xData, yData );
% legend( h, 'y2 vs. x2', 'untitled fit 2', 'Location', 'NorthEast', 'Interpreter', 'none' );
% % Label axes
% xlabel( 'x2', 'Interpreter', 'none' );
% ylabel( 'y2', 'Interpreter', 'none' );
% grid on

[xData, yData] = prepareCurveData( x2, y2 );
ft = fittype( 'fourier2' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.Robust = 'Bisquare';
opts.StartPoint = [0 0 0 0 0 0.00875317323775607];
[fitresult{3}, gof(3)] = fit( xData, yData, ft, opts );
figure( 'Name', 'untitled fit 3' );
h = plot( fitresult{3}, xData, yData );
legend( h, 'y2 vs. x2', 'untitled fit 3', 'Location', 'NorthEast', 'Interpreter', 'none' );
% Label axes
xlabel( 'x2', 'Interpreter', 'none' );
ylabel( 'y2', 'Interpreter', 'none' );
grid on

[xData, yData] = prepareCurveData( x1, y1 );
ft = fittype( 'fourier2' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.Robust = 'Bisquare';
opts.StartPoint = [0 0 0 0 0 0.017459112224018];
[fitresult{4}, gof(4)] = fit( xData, yData, ft, opts );
figure( 'Name', 'untitled fit 4' );
h = plot( fitresult{4}, xData, yData );
legend( h, 'y1 vs. x1', 'untitled fit 4', 'Location', 'NorthEast', 'Interpreter', 'none' );
% Label axes
xlabel( 'x1', 'Interpreter', 'none' );
ylabel( 'y1', 'Interpreter', 'none' );
grid on

[xData, yData] = prepareCurveData( x1, y1 );
ft = fittype( 'fourier2' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.Robust = 'Bisquare';
opts.StartPoint = [0 0 0 0 0 0.017459112224018];
[fitresult{5}, gof(5)] = fit( xData, yData, ft, opts );
figure( 'Name', 'untitled fit 5' );
h = plot( fitresult{5}, xData, yData );
legend( h, 'y1 vs. x1', 'untitled fit 5', 'Location', 'NorthEast', 'Interpreter', 'none' );
% Label axes
xlabel( 'x1', 'Interpreter', 'none' );
ylabel( 'y1', 'Interpreter', 'none' );
grid on
end

可能在第二步计算抛面的时候存在一些问题，直接用的平均值法，有点过于简单粗暴了；

然后还有第一步中极值的输出，当有多个极值时，如何选择合适的极值是算法的关键。