角度逼近圆弧插补法MATLAB与C语言实现插补与仿真

插补算法

插补（Interpolation），即机床数控系统依照一定方法确定刀具运动轨迹的过程。也可以说，已知曲线上的某些数据，按照某种算法计算已知点之间的中间点的方法，也称为“数据点的密化”；数控装置根据输入的零件程序的信息，将程序段所描述的曲线的起点、终点之间的空间进行数据密化，从而形成要求的轮廓轨迹，这种“数据密化”机能就称为“插补”。

圆弧插补的基本概念

如图2所示，圆弧P0PeP_0P_eP0Pe是圆心在坐标原点的逆圆弧段，P0(x0,y0)P_0(x_0,y_0)P0(x0,y0) 是圆弧起点，Pe(xe,ye)P_e(x_e,y_e)Pe(xe,ye)是圆弧终点，RRR是圆弧半径。FFF是刀具沿圆弧段的移动速度。圆弧插补的要求是根据编程指定的刀具移动速度F在圆弧P0PeP_0P_eP0Pe上找出若干个插补点P1,P2,P3,…,Pi,…P_1,P_2,P_3,…,P_i,…P1,P2,P3,…,Pi,…用若干个微小直线段来代替圆弧，使刀具在每个插补周期TTT内
沿小直线段从一个插补点移动到另一个插补点。
起始角：
θ0=tan−1y0x0\theta_0 =tan^{-1}\frac{y_0}{x_0} θ0=tan−1x0y0
终止角：
θe=tan−1yexe\theta_e =tan^{-1}\frac{y_e}{x_e} θe=tan−1xeye
每个插补周期半径矢量移动角:
Δθe=T⋅FR\Delta \theta_e =\frac{T \cdot F}{R} Δθe=RT⋅F
每个插补周期半径矢量与x轴正向夹角:
θi=θi−1+Δθ\theta_{i}=\theta_{i-1}+\Delta \theta θi=θi−1+Δθ
则每个插补点的计算为:
x=Rcosθy=Rcosθx=Rcos\theta\\y=Rcos\theta x=Rcosθy=Rcosθ

角度逼近圆弧插补算法

为了求cos⁡θcos⁡θcos⁡θ和sin⁡θsin⁡θsin⁡θ，可以先求
tanθ=sin⁡θcos⁡θtan\theta =\frac{sin⁡\theta}{cos⁡\theta} tanθ=cos⁡θsin⁡θ
根据角度逼近定理，将角度θ近似分解成代数和的形式
θ≈a0tan⁡−1120+a1tan⁡−1121+⋯+antan⁡−112n=θn\theta \approx a_{0} \tan ^{-1} \frac{1}{2^{0}}+a_{1} \tan ^{-1} \frac{1}{2^{1}}+\cdots+a_{n} \tan ^{-1} \frac{1}{2^{n}}=\theta_{n} θ≈a0tan−1201+a1tan−1211+⋯+antan−12n1=θn
xn2+yn2=(xn−1−an2nyn−1)2+(yn−1+an2nxn−1)2=(xn−12+yn−12)(1+122n)xn2+yn2=(xn−22+yn−22)[1+122(n−1)](1+122n)=…=(x02+y02)(1+122)(1+122×2)…[1+122(n−1)](1+122n)令 Kn=(1+122)(1+122×2)…[1+122(n−1)](1+122n),为迭代常数，与 n有关. \begin{aligned} &x_{n}^{2}+y_{n}^{2}=\left(x_{n-1}-\frac{a_{n}}{2^{n}} y_{n-1}\right)^{2}+\left(y_{n-1}+\frac{a_{n}}{2^{n}} x_{n-1}\right)^{2}=\left(x_{n-1}^{2}+y_{n-1}^{2}\right)\left(1+\frac{1}{2^{2 n}}\right)\\ &x_{n}^{2}+y_{n}^{2}=\left(x_{n-2}^{2}+y_{n-2}^{2}\right)\left[1+\frac{1}{2^{2(n-1)}}\right]\left(1+\frac{1}{2^{2 n}}\right)=\ldots=\\ &\left(x_{0}^{2}+y_{0}^{2}\right)\left(1+\frac{1}{2^{2}}\right)\left(1+\frac{1}{2^{2 \times 2}}\right) \ldots\left[1+\frac{1}{2^{2(n-1)}}\right]\left(1+\frac{1}{2^{2 n}}\right)\\ &\text { 令 } K_{n}=\left(1+\frac{1}{2^{2}}\right)\left(1+\frac{1}{2^{2 \times 2}}\right) \ldots\left[1+\frac{1}{2^{2(n-1)}}\right]\left(1+\frac{1}{2^{2 n}}\right),\text { 为迭代常数，与 } \mathrm{n} \text { 有关. } \end{aligned}xn2+yn2=(xn−1−2nanyn−1)2+(yn−1+2nanxn−1)2=(xn−12+yn−12)(1+22n1)xn2+yn2=(xn−22+yn−22)[1+22(n−1)1](1+22n1)=…=(x02+y02)(1+221)(1+22×21)…[1+22(n−1)1](1+22n1) 令 Kn=(1+221)(1+22×21)…[1+22(n−1)1](1+22n1), 为迭代常数，与 n 有关.

算法插补点的核心求法:

KaTeX parse error: Can't use function '$' in math mode at position 3: $̲\left.x_{n}^{2}…
\begin{array}{c}
r=\sqrt{\left(x_{0}^{2}+y_{0}{2}\right) K_{n}} \
x=R \cos \theta \approx R \cos \theta_{n}=R k x_{n}, y=R \sin \theta \approx R \sin \theta_{n}=R k y_{n}
\end{array}
$$

插补算法的比较

除了角度逼近圆弧插补算法外，通常使用的还有EDDA算法（数值积分法）和直接函数法。

表中提到角度逼近算法的速度误差易控制是指该误差有明确的误差公式，且可以通过限制与输入参考量无关的迭代次数n来控制。表中没有直接比较三种算法的运算速度，这是因为三种算法所用的计算方法不同，所以在不同的硬件支持条件下，它们的运算速度是不同的。例如在有快速算术协处理器的系统中，EDDA或直接函数法运算速度要快一些，而在没有硬件乘除部件的系统中，则角度逼近法运算速度要快一些。

在象限、整圆、终点判别处理中，角度逼近法较其它两种算法要方便一些。这是因为角度逼近法是以角度为插补参量，而其它两种算法是以坐标(x，y)为插补参量的。

程序的局限性及展望

目前该程序存在一定的局限性，它的插补圆心目前是固定在原点，并不能随着插补起点终点和半径自动算出，还有改程序需要输入角度，也与实际数控系统的输入有所不同，算法内的迭代次数也不能根据实际的弧长自动选择合适的迭代次数。目前的改进思路是根据现有的起点坐标，输入的终点坐标和圆弧半径，代入圆的函数解析式 (x−a)2+(y−b)2=R2\ (x-a)^2 + (y-b)^2 =R^2 (x−a)2+(y−b)2=R2
，求出圆心坐标(a,b)(a,b)(a,b)，再进行坐标变换，计算插补圆弧的角度；同时可以建立一个圆弧角度与迭代次数的关系，让迭代次数伴随着圆弧角度自动适配。

参考文献

[1] 除书法,周建来,唐学飞. 基于角度逼近法的圆弧插补算法[J]. 现代制造工程,2002,(08):13-14.
[2] 张力. 一种新的圆弧插补算法——角度逼近法[J]. 北京航空航天大学学报,1991,(01):77-84.

/***********************************************************************
角度逼近圆弧数控插补算法Editor：Ezekiel Mok/JyLi/Tiger chenFunction:          void AXY_interpretion(double xs, double ys, double xe, double ye, double x0, double y0,double ANG)Description:       在X-Y轴所构成的平面上，利用角度逼近原理计算圆弧的插补值系列Parameters:        xs, ys     起点 xe ye   终点x0,y0    圆心坐标ANG  插补圆弧角度Return Value:      大于或等余0      给予此运动命令的编码小于0            失败，传回值的意义可参考错误信息代码
***************************************************************************/
#include "stdio.h"
#include "math.h"
#include "conio.h"
#define n 1000//算法迭代次数
#define F 0.25//进给速度
## 参考文献：#define T 0.5// 插补周期
//角度均以弧度制表示
/*************************************************************************
**************************************************************************/
void AXY_interpretion(double xs, double ys, double xe, double ye, double x0, double y0, double ANG);int main()
{ double xs;
/**********************xs=atan(1/(2^0));xs=atan(1/pow(2,0));xs=1;ys=1;xs=atan(ys/xs);xs=atan(1/(2^0));double k=2,Kn=1;xs=(Kn*(1+1/pow(2,2*k)));**********************/xs=pow(2,8);printf("int is:%f\n",xs);AXY_interpretion(1.5,0,-1.5,0,0,0,180);//xs, ys(起点)xe ye(终点) x0,y0(圆心坐标) ANG(角度)return 0;
}
void AXY_interpretion(double xs, double ys, double xe, double ye, double x0, double y0, double ANG)
{   FILE *f1;double i=0,j=0,flag=1;double R; double theta_i,theta_S,theta_E,deta_theta,R0,S_0,S_i,deta_x,deta_y,A_Judge,A_Judge1,x_calulator,y_calulator;//角度逼近圆弧算法参数double x_temp,y_temp,S_j;//存放i-1的数据 double x_inter,y_inter,x_i,y_i,deta_xi,deta_yi;//插补操作参数double Kn=1,k=1;//计算逼近常数
/*******************************************************************
初始化参数
********************************************************************/R=(sqrt(pow(xs,2)+pow(ys,2)));S_0=atan(1/pow(2,0));theta_S=atan(ys/xs);printf("theta_S:%f\n",theta_S);ANG=3.14159*ANG/180;theta_E=ANG;
//  theta_E=atan(ye/xe);theta_i=theta_S;deta_theta=(F*T)/R;printf("deta_theta:%f\n",deta_theta);// f1=fopen("D:\My Documents\桌面\INTERP\c.txt","w+");f1=fopen("C:\\c.txt","w+");for(k=1;k<=n;k++){Kn=(Kn*(1+1/pow(2,2*k)));}printf("Kn:%f\n",Kn);R0=R/(sqrt(2*Kn));S_0=atan((1/pow(2,0)));printf("i:%f\n",theta_i);printf("E:%f\n",theta_E);   i=1;while(theta_i<=theta_E&&flag==1)//避免死循环 {theta_i=theta_i+deta_theta;if(theta_i>theta_E){theta_i=theta_E;flag=0;}printf("thi:%f\n",theta_i);if(theta_i<=1.5708&&theta_i>0){S_0=atan((1/pow(2,0)));A_Judge=theta_i-S_0; x_inter=R0;y_inter=R0;j=1;while(j<=n){ S_j=atan((1/pow(2,j)));x_calulator=(1/pow(2,j))*y_inter;y_calulator=(1/pow(2,j))*x_inter;if(A_Judge>=0){y_inter=y_inter+y_calulator;x_inter=x_inter-x_calulator;A_Judge=A_Judge-S_j;}else{y_inter=y_inter-y_calulator;x_inter=x_inter+x_calulator;A_Judge=A_Judge+S_j;}j++;  }x_i=x_inter;y_i=y_inter;deta_xi=x_i-x_temp;x_temp=x_i;deta_yi=y_i-y_temp;y_temp=y_i;fprintf(f1,"%f,",x_i);fprintf(f1,"%f\n",y_i);}if(theta_i>1.5708&&theta_i<=3.14159){S_0=atan((1/pow(2,0)));A_Judge=theta_i-S_0-1.5708; x_inter=-R0;y_inter=R0;j=1;while(j<=n){ S_j=atan((1/pow(2,j)));x_calulator=(1/pow(2,j))*y_inter;y_calulator=(1/pow(2,j))*x_inter;if(A_Judge>=0){y_inter=y_inter+y_calulator;x_inter=x_inter-x_calulator;A_Judge=A_Judge-S_j;}else{y_inter=y_inter-y_calulator;x_inter=x_inter+x_calulator;A_Judge=A_Judge+S_j;}j++;  }x_i=-x_inter;y_i=-y_inter;deta_xi=x_i-x_temp;x_temp=x_i;deta_yi=y_i-y_temp;y_temp=y_i;fprintf(f1,"%f,",-x_i);fprintf(f1,"%f\n",-y_i);}     if(theta_i>3.14159&&theta_i<=4.7124){S_0=atan((1/pow(2,0)));A_Judge=theta_i-S_0-3.14159; x_inter=-R0;y_inter=-R0;j=1;while(j<=n){ S_j=atan((1/pow(2,j)));x_calulator=(1/pow(2,j))*y_inter;y_calulator=(1/pow(2,j))*x_inter;if(A_Judge>=0){y_inter=y_inter+y_calulator;x_inter=x_inter-x_calulator;A_Judge=A_Judge-S_j;}else{y_inter=y_inter-y_calulator;x_inter=x_inter+x_calulator;A_Judge=A_Judge+S_j;}j++;   }x_i=x_inter;y_i=y_inter;deta_xi=x_i-x_temp;x_temp=x_i;deta_yi=y_i-y_temp;y_temp=y_i;fprintf(f1,"%f,",x_i);fprintf(f1,"%f\n",y_i);}                    if(theta_i>4.7124&&theta_i<=6.283){S_0=atan((1/pow(2,0)));A_Judge=theta_i-S_0-4.71; x_inter=R0;y_inter=-R0;j=1;while(j<=n){ S_j=atan((1/pow(2,j)));x_calulator=(1/pow(2,j))*y_inter;y_calulator=(1/pow(2,j))*x_inter;if(A_Judge>=0){y_inter=y_inter+y_calulator;x_inter=x_inter-x_calulator;A_Judge=A_Judge-S_j;}else{y_inter=y_inter-y_calulator;x_inter=x_inter+x_calulator;A_Judge=A_Judge+S_j;}j++;  }x_i=-x_inter;y_i=-y_inter;deta_xi=x_i-x_temp;x_temp=x_i;deta_yi=y_i-y_temp;y_temp=y_i;fprintf(f1,"%f,",-x_i);fprintf(f1,"%f\n",-y_i);}               i+=1;}}

clc,clear all;
plot_data=load('c:\\c.txt');
%  plot_data=load('D:\My Documents\Desktop\INTERP\My Documents桌面INTERPc.txt');
x=0:0.01:1;
y=sqrt(1-x.^2);
z=x.*0;
plot(x,y,'r');
hold on;
plot(-x,y,'r');
plot(-x,-y,'r');
plot(x,-y,'r');
grid on;
for i=1:1:size(plot_data)-1A1=plot_data(i,1);B1=plot_data(i,2);A2=plot_data(i+1,1);B2=plot_data(i+1,2);hold onif(rem(i,2)==0)line([A1,A2],[B1,B1]);pause(0.3);line([A2,A2],[B1,B2]);pause(0.3);elseline([A1,A1],[B1,B2]);pause(0.3);line([A1,A2],[B2,B2]);pause(0.3);end
end
hold on;
plot_data=load('c:\\c.txt');
plot(plot_data(:,1),plot_data(:,2),'*')

角度逼近圆弧插补法插补与仿真相关推荐

基于带约束S型加减速曲线的空间直线插补与空间圆弧插补算法（Matlab）
写在前面学习代码都记录在个人github上,欢迎关注~ Matlab机器人工具箱版本9.10 在前面的博文中: 基于抛物线过渡(梯形加减速)的空间直线插补算法与空间圆弧插补算法(Matlab) 基于 ...
两轴插补速度怎么给定_快速入门 | 篇十七：运动控制器多轴插补运动指令的使用...
以下是图文详解本节内容主要分为三大部分:插补运动.轨迹前瞻以及SP速度指令. 一插补运动插补是机床数控系统依照一定方法确定刀具运动轨迹的过程,插补是一个实时进行的数据密化的过程,不论是何种插补算 ...
插补、直线插补、联动与插补
转自:http://blog.gkong.com/218hty_27393.ashx 什么是插补? 试由直线的逐点比较工作节拍说明其插补过程: 答:插补是在组成轨迹的直线段或曲线段的起点和终点之间,按 ...
快速入门 | 篇十七：运动控制器多轴插补运动指令的使用
视频教程:<视频教程:运动控制器多轴插补运动指令的使用> 以下是图文详解本节内容主要分为三大部分:插补运动.轨迹前瞻以及SP速度指令. 一插补运动插补是机床数控系统依照一定方法确定刀 ...
数字积分法直线插补c语言程序,数控课程设计(数字积分法第二象限直线插补程序).docx...
数字积分法第二象限直线插补程序设计数字积分法是利用数字积分的方法,计算刀具沿各坐标轴的位移,使得刀具沿着所加工的轮廓曲线运动利用数字积分原理构成的插补装置称为数字积分器,又称数字微分 ...
数字积分法插补c语言程序,数字积分法插补原理
<数字积分法插补原理>由会员分享,可在线阅读,更多相关<数字积分法插补原理(19页珍藏版)>请在人人文库网上搜索. 1.3单元数字积分法插补原理,掌握数字积分法插补基本原理 ...
R语言数据缺失值处理（随机森林，多重插补）
缺失值是指数据由于种种因素导致的数据不完整,可以分为机械原因和人为原因.对于缺失值我们通常采用以下几种方法来进行插补. 1.读取数据通过read.csv函数导入文档,也可以用其他函数读入,如open ...
步进电机驱动数控十字滑台直线、圆弧插补（附斜椭圆插补代码实现）
本人大三,写个博文算是对近期课设的心得体会,供大家学习参考. 效果图: 先介绍一下此次用到的实验设备: 十字滑台,两个步进电机,一个舵机,一块STM32开发板,一支笔. 需要用到的软件是Keil5.C ...
回归插补法_用均值替换、回归插补及多重插补进行插补
##设置工作空间 #把"数据及程序" 文件夹复制到F盘下, 再用setwd设置工作空间 setwd("F: /数据及程序/chapter4/示例程序") #读取 ...

角度逼近圆弧插补法插补与仿真