[ÇóÖú] ¶àÏîʽÄâºÏ£¬±¨´í˵X£¬YÊý¾Ý³¤¶È²»Ò»Ö£¬µ«ÊǶ¨Ò峤¶ÈÒ»Ñù

??? Error using ==> polyfit at 48
X and Y vectors must be the same size.

Error in ==> Test_fitting at 36
p_Ft=polyfit(theta_Mid,Kt,Degree);   

Fx=load('Fxx.txt');
Fy=load('Fyy.txt');

n     =length(Fx);                                     % Sample quantity of Fx,Fy
theta =zeros(n,1);                                     % Rotation angle
Ft    =zeros(n,1);                                     % Tangential force initialization

for k=1:n
    theta(k)=pi/n*k;                                   % Discretize the rotation angle
end
for k=1:n
    Ft(k)=-cos(theta(k))*Fx(k)+sin(theta(k))*Fy(k);    % Evaluate Ft
end

% Calculate Kt,Kr
ft          =0.1;                                      % Feed per tooth
ap          =3;                                        % Axial cutting depth
delta_theta =10;                                       % Set the angle-cut length at both ends (Unit: degree)
delta       =round(delta_theta/180*n);                 % Quantity of samples deserted   
theta_Mid   =zeros(n-2*delta,1);                       % Angle data for fitting initialization
Kt          =zeros(n-2*delta,1);                       % Coefficient of tangential force
Ft_fit      =zeros(n-2*delta,1);                       % Fitting of Tangential force initialization


for k=delta+1:n-delta
%   Kt(k)=Ft(k)/(ft*ap);      
    Kt(k)=Ft(k)/(ap*ft*sin(theta(k)));                 % Evaluate Kt within the given interval
end

for k=1:n-2*delta
   theta_Mid(k)=theta(k+delta);                        % Cut angle data near 0 and pi
end

Degree=5;                                              % Set the polynomial degree
p_Ft=polyfit(theta_Mid,Kt,Degree);                     % Evaluate the coefficients
Kt_fit=polyval(p_Ft,theta_Mid);

for k=delta+1:n-delta
%   Ft_fit(k)=ap*ft*Kt_fit(k);      %
    Ft_fit(k)=ap*ft*sin(theta(k))*Kt_fit(k);           % Fitting expression of Ft
end

% Output
%-------------------------------------------------------------------------%
figure(1)
set(gcf,'color','w')
plot(theta,Fx,'k')                                     % Fx
hold on
plot(theta,Fy,'b')                                     % Fy
grid on
legend('Fx,Fy')
%-------------------------------------------------------------------------%
figure(2)
set(gcf,'color','w')
plot(theta,Kt,'r')                                     % Kt
hold on
plot(theta,Kt_fit,'g')                                 % Kt_fit
hold on
grid on
legend('Kt','Kt-fit')
%-------------------------------------------------------------------------%
figure(3)
set(gcf,'color','w')
plot(theta,Ft,'r')                                     % Ft
hold on
plot(theta,Ft_fit,'g')                                 % Ft_fit
hold on
grid on
legend('Ft','Ft-fit')

Fx=load('Fxx.txt');
Fy=load('Fyy.txt');

n     =length(Fx);                                            % Sample quantity of Fx,Fy
theta =zeros(n,1);                                            % Rotation angle
Ft    =zeros(n,1);                                            % Tangential force initialization

for k=1:n
    theta(k)=pi/n*k;                                          % Discretize the rotation angle
end
for k=1:n
    Ft(k)=-cos(theta(k))*Fx(k)+sin(theta(k))*Fy(k);           % Evaluate Ft
end

% Calculate Kt,Kr
ft          =0.1;                                             % Feed per tooth
ap          =3;                                               % Axial cutting depth
delta_theta =10;                                              % Set the angle-cut length at both ends (Unit: degree)
delta       =round(delta_theta/180*n);                        % Quantity of samples deserted   
theta_Mid   =zeros(n-2*delta,1);                              % Angle data for fitting initialization
Kt          =zeros(n-2*delta,1);                              % Coefficient of tangential force
Ft_fit      =zeros(n-2*delta,1);                              % Fitting of Tangential force initialization


for k=delta+1:n-delta
%   Kt(k)=Ft(k)/(ft*ap);      
    Kt(k-delta)=Ft(k)/(ap*ft*sin(theta(k)));                  % Evaluate Kt within the given interval
end

for k=1:n-2*delta
   theta_Mid(k)=theta(k+delta);                               % Cut angle data near 0 and pi
end

Degree=5;                                                     % Set the polynomial degree
p_Ft=polyfit(theta_Mid,Kt,Degree);                            % Evaluate the coefficients
Kt_fit=polyval(p_Ft,theta_Mid);

for k=delta+1:n-delta
%   Ft_fit(k)=ap*ft*Kt_fit(k);      %
    Ft_fit(k-delta)=ap*ft*sin(theta(k))*Kt_fit(k-delta);      % Fitting expression of Ft
end

% Output
%-------------------------------------------------------------------------%
figure(1)
set(gcf,'color','w')
plot(theta,Fx,'k')                                     % Fx
hold on
plot(theta,Fy,'b')                                     % Fy
grid on
legend('Fx,Fy')
%-------------------------------------------------------------------------%
figure(2)
set(gcf,'color','w')
plot(theta,Kt,'r')                                     % Kt
hold on
plot(theta,Kt_fit,'g')                                 % Kt_fit
hold on
grid on
legend('Kt','Kt-fit')
%-------------------------------------------------------------------------%
figure(3)
set(gcf,'color','w')
plot(theta,Ft,'r')                                     % Ft
hold on
plot(theta,Ft_fit,'g')                                 % Ft_fit
hold on
grid on
legend('Ft','Ft-fit')