论文中的滚动轴承故障信号仿真，自己不太熟悉，有没有给力的来做下
前期我做的轴承故障仿真程序，供参考：
function da=da(t,a)
m1=5;
m2=50; %质量
c1=1000;
c2=2000; %阻尼
kb=7.055*10^9;
k=2.5*10^9; % 刚度
Fa=20; % x方向力
Z=9; % 滚动体个数
g=9.8;
d=9.5e-3;
D=46.5e-3;
wi=500*2*pi/60;
wc=0.5*(wi*(1-d/D)); % 保持架角速度
theta1=wc*t+0;
theta2=wc*t+2*pi/Z;
theta3=wc*t+2*pi*2/Z;
theta4=wc*t+2*pi*3/Z;
theta5=wc*t+2*pi*4/Z;
theta6=wc*t+2*pi*5/Z;
theta7=wc*t+2*pi*6/Z;
theta8=wc*t+2*pi*7/Z;
theta9=wc*t+2*pi*8/Z;
theta_ball=[theta1',theta2',theta3',theta4',theta5',theta6',theta7',theta8',theta9'];
%%%赫兹接触力
u = 1e-4; % 径向游隙
sigma_ball= (a(1)-a(2)).*cos(theta_ball)+(a(3)-a(4)).*sin(theta_ball) - u; %%%各滚珠的法向接触变形量
Heaviside_1 = heaviside(sigma_ball); % heaviside 函数用于筛选 sigma_ball>0的变形量
useful_sigma = sigma_ball.*Heaviside_1;
sigma_1=useful_sigma.^3;
sigma = kb * sqrt(sigma_1);
fix_bearing = (sigma.*cos(theta_ball)); % 轴承支撑力的x方向分量
fiy_bearing = (sigma.*sin(theta_ball)); % 轴承支撑力的y方向分量
Fx= sum(fix_bearing);
Fy= sum(fiy_bearing);
%%% 总公式
da=zeros(8,1);
da(1)=a(5);%x1'
da(2)=a(6);%x2'
da(3)=a(7);%y1'
da(4)=a(8);%y2'
da(5)=(Fa+m1*g-Fx-c1*( a(5)-a(6) ))/m1; %x1''
da(7)=(-Fy-c1*( a(7)-a(8) ))/m1;%y1''
da(6)=(m2*g+Fx-k*a(2)-c2*a(6)+c1*(a(5)-a(6)))/m2; %x2''
da(8)=(Fy-k*a(4)+c1*(a(7)-a(8))-c2*a(8))/m2;%y2''
end
clc;
clear;
dt=0.0001;%步长
t=0:dt:1.5;%时间
a=zeros(1,8);
[t,a]=ode45(@da,t,a);
figure(1)
plot(t,a(:,1))
title('x1方向位移')
figure(2)
plot(t,a(:,2))
title('x2方向位移')
figure(3)
plot(t,a(:,3))
title('y1方向速度')
figure(4)
plot(t,a(:,4))
title('y2方向速度')
% 频谱分析
X = a(end-1000:end,1) ;
X = X - mean( X );
[f,X_m,X_fai] = DFT(X,dt);
figure
plot( f,X_m )
xlim([0 300])
xlabel('频率/Hz')
ylabel('幅值/mm')