clear all;clc;

load ('MOnixinhao.mat');

%处理键相信号获得键相出发时间点arrive_time

[maxtab,mintab]=peakdet(tacho,0.5);

arrive_time=t0(maxtab(:,1)');

%通过方程式拟合获取角度域对应的时间节点angle_time。以每周期采样256个数据点为准.需要去掉第一个和最后一个时间节点.

%从第二个键相时间点开始往前pi往后pi的方式开始采样。%%作为参考，一共转了140圈，所以总角度不应该大于140*2*pi=879.6459

angle_time=[];

for k=2:(length(arrive_time)-1);

    A=[1,arrive_time(k-1),arrive_time(k-1)^2;1,arrive_time(k),arrive_time(k)^2;1,arrive_time(k+1),arrive_time(k+1)^2];

    B=[0;2*pi;4*pi]+(k-2)*2*pi;

    C=A\B;

    

    N=256;n=0:255;

    angle_trin=B(1,:)+pi+2*pi*(n./N);

    if C(3)~=0

        this_time=((4*C(3).*(angle_trin-C(1))+C(2)^2).^(0.5)-C(2))/(2*C(3));

    elseif C(3)==0 

        this_time=(angle_trin-C(1))./C(2);

    end

    angle_time=[angle_time,this_time];

end

%对角域对应的时间节点进行插值计算，插值方法为多项式插值，获取角域数据序列angle_val_trin 

angle_val_trin=[];

for t=1:length(angle_time);

    [minval,minpos]=min(abs(t0-angle_time(t)));

    int_y=[1,t0(minpos-1),t0(minpos-1)^2,t0(minpos-1)^3;1,t0(minpos),t0(minpos)^2,t0(minpos)^3;

        1,t0(minpos+1),t0(minpos+1)^2,t0(minpos+1)^3;1,t0(minpos+2),t0(minpos+2)^2,t0(minpos+2)^3];

    int_t=[x(minpos-1);x(minpos);x(minpos+1);x(minpos+2)];

    int_a=int_y\int_t;

    this_val=int_a(1)+int_a(2)*angle_time(t)+int_a(3)*angle_time(t)^2+int_a(4)*angle_time(t)^3;

    angle_val_trin=[angle_val_trin,this_val];

end

%阶域的各参数。总阶数256阶，对应angle_val_trin的长度进行计算

ordertrick=abs(fft(angle_val_trin));

ordertrick=ordertrick(1:(length(ordertrick)/2));

order_num=0:256/(length(angle_val_trin)-1):128;

plot(order_num,ordertrick);

问题一大把 慢慢解决中