test_three_point_frequency_domain.m

matlab圆度误差的评定算法

clc;
clear;
%三点频域法
%注意m、p的取值必须互质，G(k)的零值需要剔除

%权函数G(k)
N=64;%采样点数
n=0:N-1;%每圈的采样点数
theta=n*2*pi/N;%采样点的相位角
m=13;p=16;%注意m、p的取值必须互质
alpha=m*2*pi/N;
bata=p*2*pi/N;%第一个传感器安装在alpa处，第二个传感器安装在（alpha＋bata）角处
rr=10;%基圆半径
k=0:N-1;
a=-sin(2*pi*(m+p)/N)/sin(2*pi*p/N);
b=sin(2*pi*m/N)/sin(2*pi*p/N);
WN=exp(-j*2*pi/N);
Gk=1+a*WN.^(-m*k)+b*WN.^(-(m+p)*k);%k=1、N－1时Gk＝0
Gk(2)=1;Gk(N)=1;%剔除零值

%读取文件 
[filename,filepath]=uigetfile('*.txt','Select Input data file'); 
file = [filepath filename]; 
fid = fopen(file,'r'); 
if fid == -1 
disp('Error opening the file') 
end 
i=1;
while 1 
nextline = fgetl(fid); %读第一行 
if ~ischar(nextline), break, end %读到最后跳出 
g(:,i) = sscanf(nextline, ' %f %f ')';%按行读取后根据你自己的需要改为按行或按列输出数据
i=i+1;
end 
fclose(fid);

%传感器输入数据进行集合平均平滑噪声信号
Long=length(g);
nn=Long/192;%传感器采集的数据圈数
reaa=reshape(g,N*3,nn);
meanaa=sum(reaa,2)/nn;%按行求和得到列向量
aa=meanaa(1:end/3)';ab=meanaa(end/3+1:end*2/3)';ac=meanaa(end*2/3+1:end)';
sprintf('A通道的跳动量为: %0.4f，C通道的跳动量为：%0.4f',td1,td2)

sa=aa-sum(aa)/N;sb=ab-sum(ab)/N;sc=ac-sum(ac)/N;%各个传感器测量值减去直流分量
sn=sa+a*sb+b*sc;%sn为信号的加权组合

%离散快速傅立叶变换 S(k)=FFT(S(n))
k=0:N-1;
Sk=fft(sn,N);
Lc=0.5;%Lc为传感器测头的量程
if Sk(1)<0.03*Lc
    Sk(2)=0;Sk(N)=0;%过滤一阶谐波分量（由主轴安装偏心造成）
else
    disp('测试信号中混入了较大的噪声信号，该次测量数据无效');
    break;
end
figure
subplot(311)
plot(n,sa)
xlabel('n');
ylabel('sa');
title('A通道时域图')
subplot(312)
plot(n,sb)
xlabel('n');
ylabel('sb');
title('B通道时域图')
subplot(313)
plot(n,sc)
xlabel('n');
ylabel('sc');
title('C通道时域图')

%傅立叶反变换 H(n)=IFFT(S(k)/G(k))
Hk=Sk./Gk;
hn=ifft(Hk,N);
hn=real(hn);
figure
subplot(211),stem(k,abs(Hk));
xlabel('k');
ylabel('abs(Hk)');
title('频域法形状误差幅频图')
subplot(212),stem(n,hn);
xlabel('n');
ylabel('real(hn)');
title('频域法形状误差时域图')
deltax=sa-hn;
deltay=(sb-[hn(m+1:N),hn(1:m)]-deltax.*cos(m*2*pi/N))./sin(m*2*pi/N);
        
figure
polar(theta,rr+hn);
title('频域法分离的形状误差图')

%figure
[angle,rh]=cart2pol(deltax,deltay);%atan2命令可以求－180～＋180的反正切
%stem(angle*180/pi,rh);
%title('分离后的回转误差')
%xlabel('回转误差相位角')
%ylabel('回转误差误差值')

figure
plot(n,angle*180/pi) 
title('回转误差相位随时间的变化图')
figure
plot(n,rh)
title('回转误差大小随时间的变化图')

%圆度误差的评定
ess_lsc=Roundness_Assessment(rr+hn,rr)
title('最小二乘法圆度误差评定')
%回转误差的评定
ess_rot=Roundness_Assessment(rr+sqrt(deltay.^2+deltax.^2),rr)
title('最小二乘法回转误差评定')