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('最小二乘法回轉誤差評定')