大唐项目存稿

%% 将表格导入，并将行转换成列
AngleData1=AngleData';
PositionData1=PositionData';
%下面是三个点的情况。若点的个数变化需要
x1=AngleData1(:,1);
y1=PositionData1(:,1);

x2=AngleData1(:,2);
y2=PositionData1(:,2);

x3=AngleData1(:,3);
y3=PositionData1(:,3);

%将数组的第一个点给去掉，因为是零，会极大影响拟合效果
x1(1, :) = [];
y1(1, :) = [];

x2(1, :) = [];
y2(1, :) = [];

x3(1, :) = [];
y3(1, :) = [];

%筛选出非零的行（列）数，即有效测量点的数量
n=0;
a=size(PositionData);
for u=find(PositionData(1:a(1))~=0)
    n=n+1;
end
%n是有效测量点的数量
%%

%对第一个测量点进行筛选
[valae_x1_min,position_x1_min]=min(x1);%valae_x1_min是数组x_1中的最小值，position_x1_min是数组x_1中的序数，即第几个数
x1(position_x1_min, :) = [];
y1(position_x1_min, :) = [];%这两行代码是去掉数组x_1，y_1中对应序列的那个数

%对第2个测量点进行筛选
[valae_x2_min,position_x2_min]=min(x2);
x2(position_x2_min, :) = [];
y2(position_x2_min, :) = [];

%对第3个测量点进行筛选
[valae_x3_min,position_x3_min]=min(x3);
x3(position_x3_min, :) = [];
y3(position_x3_min, :) = [];

%% 对数据进行拟合，采取的是傅里叶拟合。根据方差可知拟合效果较好（如果想看数据点和拟合曲线的话可以将注释的plot取消注释）
%第一个点
[a_1,~]=createFits(x1,y1,x2,y2);%调用数据拟合函数
FUNCTION_1=a_1{5,1}.a0 + a_1{5,1}.a1*cos(x3*a_1{5,1}.w)+ a_1{5,1}.b1*sin(x3*a_1{5,1}.w) +a_1{5,1}.a2*cos(2*x3*a_1{5,1}.w) + a_1{5,1}.b2*sin(2*x3*a_1{5,1}.w);
%FUNCTION_1是第一个测量点的拟合曲线
[pks_po_1,locs_po_1]=findpeaks(FUNCTION_1);%pks_po_1，locs_po_1分别为拟合曲线的极大值点和极大值点所对应的相角（注意：极值点一般有多个）
[pks_pa_1,locs_pa_1]=findpeaks(-FUNCTION_1);%pks_pa_1，locs_pa_1分别为拟合曲线的极小值点和极小值点所对应的相角（注意：极值点一般有多个）
%下面的plot一直到hold on是画拟合曲线，并且极大值点用○圈起来，极小值点用※圈起来
figure( 'Name', 'untitled fit 1' );
plot(FUNCTION_1);
hold on
plot(locs_po_1,pks_po_1,'o');
[pks_1,~]=findpeaks(-FUNCTION_1);
plot(locs_pa_1,-pks_1,'*');
hold on;
%因为极值点的个数一般有多个，所以要进行判断。当极大值个数大于一时，最大的极大值就是最大值；当极大值个数等于一时，极大值就是最大值
if length(pks_po_1)>1
[pks_po_1,weizhi_max_1]=max(pks_po_1);%极大值会生成一个数组，pks_po_1是这些极大值中的最大值，然后weizhi_max_1是这个数组中最大值所对应数组中的序数
locs_po_1=locs_po_1(weizhi_max_1);%已知最大值点的序数，就能求出对应的相角
end
if length(pks_pa_1)>1
 [pks_pa_1,weizhi_min_1]=max(-pks_pa_1);%这里采取的是求极小值的最大值，目的是去除凹槽的影响
locs_pa_1=locs_pa_1(weizhi_min_1);
end

%第二个点
[a_2,b_2]=createFits(x2, y2,x3,y3);
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);
figure( 'Name', 'untitled fit 2' );
plot(FUNCTION_2);
hold on;
plot(locs_po_2,pks_po_2,'o');
[pks_2,~]=findpeaks(-FUNCTION_2);
plot(locs_pa_2,-pks_2,'*');
hold on;
if length(pks_po_2)>1
    [pks_po_2,weizhi_max_2]=max(pks_po_2);
     locs_po_2=locs_po_2(weizhi_max_2);
end
if length(pks_pa_2)>1
    [pks_pa_2,weizhi_min_2]=max(-pks_pa_2);
     locs_pa_2=locs_pa_2(weizhi_min_2);
end

%第三个点
[a_3,b_3]=createFits(x3, y3,x1,y1);
FUNCTION_3=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_3);
[pks_pa_3,locs_pa_3]=findpeaks(-FUNCTION_3);
figure( 'Name', 'untitled fit 3' );
plot(FUNCTION_3);
hold on
plot(locs_po_3,pks_po_3,'o');
[pks_2,locs_2]=findpeaks(-FUNCTION_3);
plot(locs_pa_3,-pks_2,'*');
hold on;
if length(pks_po_3)>1
    [pks_po_3,weizhi_max_3]=max(pks_po_3);
     locs_po_3=locs_po_3(weizhi_max_3);
end
if length(pks_pa_3)>1
    [pks_pa_3,weizhi_min_3]=max(-pks_pa_3);
     locs_pa_3=locs_pa_3(weizhi_min_3);
end

%% %至此，上面第一步已经完成，输出了每个点的极值以及相位；下面引入输出变量out1,out2,out3和out4
locs_pa=[locs_pa_1,locs_pa_2,locs_pa_3];
pks_pa=[pks_pa_1,pks_pa_2,pks_pa_3];
pks_po=[pks_po_1,pks_po_2,pks_po_3];

out1=[1];
out2=[locs_pa(1)];
out3=[abs(pks_pa(1))];

for i=1:1:n
out1=i;
out2=locs_pa(i);
out3=abs(pks_pa(i));

Out1(i)=out1(1);
Out2(i)=out2(1);
Out3(i)=out3(1);
end
disp('测量点序号为');
disp(Out1);
disp('测量点相位角');
 disp(Out2);
disp('测量点最值');
disp(Out3);

%% 调用的数据拟合的函数（傅里叶拟合）
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