Help with computing FFT

*Please help me with my issue, it may seem that my question is repeated, but it is not. My code has been updating each week and now I face a different issue.

I am modelling a helical gear pair with transmission error and want to see the frequency content of the transmission error, which should include the gear meshing frequencies and rotational speed. Hence, I need to find the FFT of 'TE'. I have read https://uk.mathworks.com/help/matlab/ref/fft.html. and try to do the FFT in my function call.

My main code-

function [yp,TE,KS] = spur1(t,y,Torque)
yp = zeros(12,1);                           %vector of 12 equations
beta = 32*(pi/180);                      %Helix angle (rad) /13
P = 18.5*(pi/180);                         %Pressure angle (rad)
speed = 1000/60; 
Freq = 1000*20/60;


n_a = 22;                           % no of driver gear teeth  
n_p = 59;                           % no of driven gear teeth  
m = 2e-3*cos(beta);
R_a = (m*n_a)/2;                         %Radius
R_p = (m*n_p)/2;                         %Radius
i_r = n_p/n_a;                        % gear ratio




% Driver gear  
m_a = 5.3e-2;                      %mass of driver gear (kg)
c_ay=6.3e2;                       %Damping of driver gear in y direction (Ns/m)
c_az=1.3e2;                       %Damping of driver gear in z direction (Ns/m)
k_ay= 1.9e8;                      %Stiffness of driver gear in y direction (N/m)
k_az= 8.3e6;                      %Stiffness of driver gear in z direction (N/m)
I_a = 1.7e-5;          %Moment of inertia of driver gear (Kg.m3)

% Driven gear  
m_p = 2.5e-1;                                  %mass of driven gear (kg)
c_py=3.6e3;                                   %Damping of driven gear in y direction (Ns/m)
c_pz=4.3e2;                                   %Damping of driven gear in z direction (Ns/m)
k_py = 8.3e8;                                %Stiffness of driven gear in y direction (N/m)
k_pz = 1.2e7;                                  %Stiffness of driven gear in z direction (N/m)
I_p = 3.9e-4;                      %Moment of inertia of driven gear (Kg.m3)


Torque = Torque(t)/1000;
Opp_Torque = 68.853/1000;                        % Average torque value

% e = (pi*2*(R_a + R_p)*tan(beta))/(4*cos(P));
e= 20e-6 + 20e-8*sin(2*pi*Freq*t);  

ks = 0.9e8 + 0.9e6*sin(2*pi*Freq*t);                   %Shear stiffness 
% Cs = 0.1 + 0.0001*sin(2*pi*Freq*t);                     %Shear damping
Cs = 2*0.1*sqrt((ks*I_a*I_p)/(I_a*R_a + I_p*R_p)) + 2e-2*0.1*sqrt((ks*I_a*I_p)/(I_a*R_a + I_p*R_p))*sin(2*pi*Freq*t);
TE = -y(11)*R_p + y(5)*R_a ;                   %transmission error
b = 20e-6;                                       %Backlash
if(TE>b)
    h = TE - b;
        KS = h*ks;
elseif(TE < -1*b)
        h = TE + b;
        KS = h*ks;
else 
    h = 0;
   KS = h*ks; 
end

Fy = cos(beta)*(Cs*(cos(beta)*(y(2) - y(8) + R_a*y(6) - R_p*y(12)) + sin(beta)*(y(4) - y(10))) + ks*(cos(beta)*(y(1) - y(7) + R_a*y(5) - R_p*y(11)) + sin(beta)*(y(3) - y(9)))); %circumferential excitation force- formula from paper
Fz = sin(beta)*(Cs*(cos(beta)*(y(2) - y(8) + R_a*y(6) - R_p*y(12)) + sin(beta)*(y(4) - y(10))) + ks*(cos(beta)*(y(1) - y(7) + R_a*y(5) - R_p*y(11)) + sin(beta)*(y(3) - y(9)))); %circumferential excitation force- formula from paper

%Driver gear equations
yp(1) = y(2); 
yp(2) = (-Fy - k_ay*y(1) - c_ay*y(2))/m_a;  %acceleration in y (up and down) direction (m/s2)
yp(3) = y(4);  
yp(4) = (-Fz - k_az*y(3) - c_az*y(4))/m_a;  %acceleration in z (side to side) direction (m/s2)
yp(5) = y(6);  
yp(6) = (Torque - Fy*R_a)/I_a; %angular acceleration
%Driven gear equations
yp(7) = y(8);  
yp(8) = (Fy - k_py*y(7) - c_py*y(8))/m_p;            %acceleration in y (up and down) direction (m/s2)
yp(9) = y(10);  
yp(10) = (Fz - k_pz*y(9) - c_pz*y(10))/m_p;          %acceleration in y (side to side) direction (m/s2)
yp(11) = y(12);  
yp(12) = (-Opp_Torque*i_r + Fy*R_p)/I_p;                % angular accelration (rad/s2) 
end

My function call where I try to compute FFT-

A = load("torque_for_Sid_60s.mat");
T_a = A.torque_60s(1:5:600001);                  %Torque (Nm) chnages with time step
speed = 1000/60; 
tspan=0:0.00001*5:6;    %can try 4 or 5
y0 = [0;0;0;0;0;104.719;0;0;0;0;0;104.719/3]; 

tic
%These are constant for every call to reduced2 so compute them here and
%pass them into the modified function reduced3
time = 0:0.00001*5:6;
Torque = griddedInterpolant(time, T_a);
[t,y] = ode45(@(t,y) spur1(t,y,Torque),tspan,y0); 
TE= zeros(size(t));
KS = zeros(size(t));
for i=1:numel(t)
  [~,TE(i),KS(i)] = spur1(t(i),y(i,:),Torque);
end
nexttile
plot(t,TE,"magenta")
xlabel('Time (s)')
ylabel('Transmission Error TE (meter)')
title('Transmission Error vs time')
axis padded
grid on

Fs = 1/(t(2)-t(1)); % Sampling frequency
L = length(t); % Length of signal
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(TE, NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

nexttile;
plot(f, 2*abs(Y(1:NFFT/2+1)))
xlabel('Frequency (Hz)')
ylabel('Amplitude')
title('FFT of TE vs frequency')
axis padded
grid on

% nexttile
% plot(t,KS,"red")
% xlabel('Time (s)')
% ylabel('KS (N)')
% title('KS force vs time')
% axis padded
% grid on
newtime = toc


% %Driver gear graphs
% nexttile  
% plot(t,y(:,2))
% ylabel('(y) up and down velocity (m/s)')
% xlabel('Time')
% title('Driver gear velocity in y direction vs time')
% axis padded  
% grid on

% nexttile  
% plot(t,y(:,4))  
% ylabel('(z) side to side velocity (m/s)')
% xlabel('Time')
% title('Driver gear velocity in z direction vs time')  
% axis padded  
% grid on

% nexttile  
% plot(t,y(:,6))  
% ylabel('angular displacement (rad)')
% xlabel('Time')
% title('Driver gear angular displacment vs time')  
% axis padded  
% grid on


% % Driven gear graphs
% nexttile  
% plot(t,y(:,8),"green")  
% ylabel('(y) up and down velocity (m/s)')
% xlabel('Time')
% title('Driven gear velocity in y direction vs time')  
% axis padded  
% grid on

% nexttile  
% plot(t,y(:,10),"green")  
% ylabel('(z) side to side velocity (m/s)')
% xlabel('Time')
% title('Driven gear velocity in z direction vs time')  
% axis padded  
% grid on

% nexttile  
% plot(t,y(:,12),"green")  
% ylabel('angular displacement (rad)')
% xlabel('Time')
% title('Driven gear angular displacement vs time')  
% axis padded  
% grid on

My torque file - torque_for_Sid_60s.mat

Answer 1

[Mathieu NOE]

hello

here some suggestions and mods I did in your code

• remove the initial transient (t< 0.1 s)
• detrend data (so you don't have that large DC output i the fft spectrum)
• display fft magnitude in log scale

Now we see better the fft peaks

(you can use findpeaks to get them)

clearvars
clc
close all

A = load("torque_for_Sid_60s.mat");
T_a = A.torque_60s(1:5:600001);                  %Torque (Nm) chnages with time step
speed = 1000/60; 
tspan=0:0.00001*5:6;    %can try 4 or 5
y0 = [0;0;0;0;0;104.719;0;0;0;0;0;104.719/3]; 
tic
%These are constant for every call to reduced2 so compute them here and
%pass them into the modified function reduced3
time = 0:0.00001*5:6;
Torque = griddedInterpolant(time, T_a);
[t,y] = ode45(@(t,y) spur1(t,y,Torque),tspan,y0); 
TE= zeros(size(t));
KS = zeros(size(t));
for i=1:numel(t)
  [~,TE(i),KS(i)] = spur1(t(i),y(i,:),Torque);
end

% MN : remove first 0.1 second (transient)
idx = (t<0.1);
t(idx) = [];
TE(idx) = [];


nexttile
plot(t,TE,"magenta")
xlabel('Time (s)')
ylabel('Transmission Error TE (meter)')
title('Transmission Error vs time')
axis padded
grid on
Fs = 1/(t(2)-t(1)); % Sampling frequency
L = length(t); % Length of signal
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(detrend(TE), NFFT)/L;  % detrend the signal before fft to remove the large DC (and near DC) fft output
f = Fs/2*linspace(0,1,NFFT/2+1);
nexttile;
semilogy(f, 2*abs(Y(1:NFFT/2+1)))
xlabel('Frequency (Hz)')
ylabel('Amplitude')
title('FFT of TE vs frequency')
axis padded
grid on
% nexttile
% plot(t,KS,"red")
% xlabel('Time (s)')
% ylabel('KS (N)')
% title('KS force vs time')
% axis padded
% grid on
newtime = toc


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [yp,TE,KS] = spur1(t,y,Torque)
yp = zeros(12,1);                           %vector of 12 equations
beta = 32*(pi/180);                      %Helix angle (rad) /13
P = 18.5*(pi/180);                         %Pressure angle (rad)
speed = 1000/60; 
Freq = 1000*20/60;
n_a = 22;                           % no of driver gear teeth  
n_p = 59;                           % no of driven gear teeth  
m = 2e-3*cos(beta);
R_a = (m*n_a)/2;                         %Radius
R_p = (m*n_p)/2;                         %Radius
i_r = n_p/n_a;                        % gear ratio
% Driver gear  
m_a = 5.3e-2;                      %mass of driver gear (kg)
c_ay=6.3e2;                       %Damping of driver gear in y direction (Ns/m)
c_az=1.3e2;                       %Damping of driver gear in z direction (Ns/m)
k_ay= 1.9e8;                      %Stiffness of driver gear in y direction (N/m)
k_az= 8.3e6;                      %Stiffness of driver gear in z direction (N/m)
I_a = 1.7e-5;          %Moment of inertia of driver gear (Kg.m3)
% Driven gear  
m_p = 2.5e-1;                                  %mass of driven gear (kg)
c_py=3.6e3;                                   %Damping of driven gear in y direction (Ns/m)
c_pz=4.3e2;                                   %Damping of driven gear in z direction (Ns/m)
k_py = 8.3e8;                                %Stiffness of driven gear in y direction (N/m)
k_pz = 1.2e7;                                  %Stiffness of driven gear in z direction (N/m)
I_p = 3.9e-4;                      %Moment of inertia of driven gear (Kg.m3)
Torque = Torque(t)/1000;
Opp_Torque = 68.853/1000;                        % Average torque value
% e = (pi*2*(R_a + R_p)*tan(beta))/(4*cos(P));
e= 20e-6 + 20e-8*sin(2*pi*Freq*t);  
ks = 0.9e8 + 0.9e6*sin(2*pi*Freq*t);                   %Shear stiffness 
% Cs = 0.1 + 0.0001*sin(2*pi*Freq*t);                     %Shear damping
Cs = 2*0.1*sqrt((ks*I_a*I_p)/(I_a*R_a + I_p*R_p)) + 2e-2*0.1*sqrt((ks*I_a*I_p)/(I_a*R_a + I_p*R_p))*sin(2*pi*Freq*t);
TE = -y(11)*R_p + y(5)*R_a ;                   %transmission error
b = 20e-6;                                       %Backlash
if(TE>b)
    h = TE - b;
        KS = h*ks;
elseif(TE < -1*b)
        h = TE + b;
        KS = h*ks;
else 
    h = 0;
   KS = h*ks; 
end
Fy = cos(beta)*(Cs*(cos(beta)*(y(2) - y(8) + R_a*y(6) - R_p*y(12)) + sin(beta)*(y(4) - y(10))) + ks*(cos(beta)*(y(1) - y(7) + R_a*y(5) - R_p*y(11)) + sin(beta)*(y(3) - y(9)))); %circumferential excitation force- formula from paper
Fz = sin(beta)*(Cs*(cos(beta)*(y(2) - y(8) + R_a*y(6) - R_p*y(12)) + sin(beta)*(y(4) - y(10))) + ks*(cos(beta)*(y(1) - y(7) + R_a*y(5) - R_p*y(11)) + sin(beta)*(y(3) - y(9)))); %circumferential excitation force- formula from paper
%Driver gear equations
yp(1) = y(2); 
yp(2) = (-Fy - k_ay*y(1) - c_ay*y(2))/m_a;  %acceleration in y (up and down) direction (m/s2)
yp(3) = y(4);  
yp(4) = (-Fz - k_az*y(3) - c_az*y(4))/m_a;  %acceleration in z (side to side) direction (m/s2)
yp(5) = y(6);  
yp(6) = (Torque - Fy*R_a)/I_a; %angular acceleration
%Driven gear equations
yp(7) = y(8);  
yp(8) = (Fy - k_py*y(7) - c_py*y(8))/m_p;            %acceleration in y (up and down) direction (m/s2)
yp(9) = y(10);  
yp(10) = (Fz - k_pz*y(9) - c_pz*y(10))/m_p;          %acceleration in y (side to side) direction (m/s2)
yp(11) = y(12);  
yp(12) = (-Opp_Torque*i_r + Fy*R_p)/I_p;                % angular accelration (rad/s2) 
end

[Siddharth Jain]

Thank you for your answer, How can I get linear scale on the y axis instead of log?

[Mathieu NOE]

hello again

use plot instead of semilogy (my suggestion) here :

semilogy(f, 2*abs(Y(1:NFFT/2+1)))

[Siddharth Jain]

Thank you!

Just one more thing, if you could also help me with computing the FFT of torque I would be very grateful. If I modify your suggestion by replacing TE by torque, its doesn't work

[Mathieu NOE]

hello again

here your are my friend

clearvars
clc
close all
A = load("torque_for_Sid_60s.mat");
T_a = A.torque_60s(1:5:600001);                  %Torque (Nm) chnages with time step
speed = 1000/60; 
tspan=0:0.00001*5:6;    %can try 4 or 5
y0 = [0;0;0;0;0;104.719;0;0;0;0;0;104.719/3]; 
tic
%These are constant for every call to reduced2 so compute them here and
%pass them into the modified function reduced3
time = 0:0.00001*5:6;
Torque = griddedInterpolant(time, T_a);
[t,y] = ode45(@(t,y) spur1(t,y,Torque),tspan,y0); 
TE= zeros(size(t));
KS = zeros(size(t));
for i=1:numel(t)
  [~,TE(i),KS(i)] = spur1(t(i),y(i,:),Torque);
end

% MN : remove first 0.1 second (transient)
idx = (t<0.1);
t(idx) = [];
TE(idx) = [];
T_a(idx) = [];

nexttile
plot(t,TE,"magenta")
xlabel('Time (s)')
ylabel('Transmission Error TE (meter)')
title('Transmission Error vs time')
axis padded
grid on
Fs = 1/(t(2)-t(1)); % Sampling frequency
L = length(t); % Length of signal
NFFT = 2^nextpow2(L); % Next power of 2 from length of y


%%
TE_fft = fft(detrend(TE), NFFT)/L;  % detrend the signal before fft to remove the large DC (and near DC) fft output
f = Fs/2*linspace(0,1,NFFT/2+1);
nexttile;
plot(f, 2*abs(TE_fft(1:NFFT/2+1)))
xlabel('Frequency (Hz)')
ylabel('Amplitude')
title('FFT of TE vs frequency')
axis padded
grid on

%%
T_a_fft = fft(detrend(T_a), NFFT)/L;  % detrend the signal before fft to remove the large DC (and near DC) fft output
nexttile;
plot(f, 2*abs(T_a_fft(1:NFFT/2+1)))
xlabel('Frequency (Hz)')
ylabel('Amplitude')
title('FFT of Torque vs frequency')
axis padded
grid on

% nexttile
% plot(t,KS,"red")
% xlabel('Time (s)')
% ylabel('KS (N)')
% title('KS force vs time')
% axis padded
% grid on
newtime = toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [yp,TE,KS] = spur1(t,y,Torque)
yp = zeros(12,1);                           %vector of 12 equations
beta = 32*(pi/180);                      %Helix angle (rad) /13
P = 18.5*(pi/180);                         %Pressure angle (rad)
speed = 1000/60; 
Freq = 1000*20/60;
n_a = 22;                           % no of driver gear teeth  
n_p = 59;                           % no of driven gear teeth  
m = 2e-3*cos(beta);
R_a = (m*n_a)/2;                         %Radius
R_p = (m*n_p)/2;                         %Radius
i_r = n_p/n_a;                        % gear ratio
% Driver gear  
m_a = 5.3e-2;                      %mass of driver gear (kg)
c_ay=6.3e2;                       %Damping of driver gear in y direction (Ns/m)
c_az=1.3e2;                       %Damping of driver gear in z direction (Ns/m)
k_ay= 1.9e8;                      %Stiffness of driver gear in y direction (N/m)
k_az= 8.3e6;                      %Stiffness of driver gear in z direction (N/m)
I_a = 1.7e-5;          %Moment of inertia of driver gear (Kg.m3)
% Driven gear  
m_p = 2.5e-1;                                  %mass of driven gear (kg)
c_py=3.6e3;                                   %Damping of driven gear in y direction (Ns/m)
c_pz=4.3e2;                                   %Damping of driven gear in z direction (Ns/m)
k_py = 8.3e8;                                %Stiffness of driven gear in y direction (N/m)
k_pz = 1.2e7;                                  %Stiffness of driven gear in z direction (N/m)
I_p = 3.9e-4;                      %Moment of inertia of driven gear (Kg.m3)
Torque = Torque(t)/1000;
Opp_Torque = 68.853/1000;                        % Average torque value
% e = (pi*2*(R_a + R_p)*tan(beta))/(4*cos(P));
e= 20e-6 + 20e-8*sin(2*pi*Freq*t);  
ks = 0.9e8 + 0.9e6*sin(2*pi*Freq*t);                   %Shear stiffness 
% Cs = 0.1 + 0.0001*sin(2*pi*Freq*t);                     %Shear damping
Cs = 2*0.1*sqrt((ks*I_a*I_p)/(I_a*R_a + I_p*R_p)) + 2e-2*0.1*sqrt((ks*I_a*I_p)/(I_a*R_a + I_p*R_p))*sin(2*pi*Freq*t);
TE = -y(11)*R_p + y(5)*R_a ;                   %transmission error
b = 20e-6;                                       %Backlash
if(TE>b)
    h = TE - b;
        KS = h*ks;
elseif(TE < -1*b)
        h = TE + b;
        KS = h*ks;
else 
    h = 0;
   KS = h*ks; 
end
Fy = cos(beta)*(Cs*(cos(beta)*(y(2) - y(8) + R_a*y(6) - R_p*y(12)) + sin(beta)*(y(4) - y(10))) + ks*(cos(beta)*(y(1) - y(7) + R_a*y(5) - R_p*y(11)) + sin(beta)*(y(3) - y(9)))); %circumferential excitation force- formula from paper
Fz = sin(beta)*(Cs*(cos(beta)*(y(2) - y(8) + R_a*y(6) - R_p*y(12)) + sin(beta)*(y(4) - y(10))) + ks*(cos(beta)*(y(1) - y(7) + R_a*y(5) - R_p*y(11)) + sin(beta)*(y(3) - y(9)))); %circumferential excitation force- formula from paper
%Driver gear equations
yp(1) = y(2); 
yp(2) = (-Fy - k_ay*y(1) - c_ay*y(2))/m_a;  %acceleration in y (up and down) direction (m/s2)
yp(3) = y(4);  
yp(4) = (-Fz - k_az*y(3) - c_az*y(4))/m_a;  %acceleration in z (side to side) direction (m/s2)
yp(5) = y(6);  
yp(6) = (Torque - Fy*R_a)/I_a; %angular acceleration
%Driven gear equations
yp(7) = y(8);  
yp(8) = (Fy - k_py*y(7) - c_py*y(8))/m_p;            %acceleration in y (up and down) direction (m/s2)
yp(9) = y(10);  
yp(10) = (Fz - k_pz*y(9) - c_pz*y(10))/m_p;          %acceleration in y (side to side) direction (m/s2)
yp(11) = y(12);  
yp(12) = (-Opp_Torque*i_r + Fy*R_p)/I_p;                % angular accelration (rad/s2) 
end

[Siddharth Jain]

Hi @Mathieu NOE, I just realised that this code is doing FFT on the unmodified input torque. Actually I want to perform the FFT on the torque that is modified in the main code -'Torque' variable

Torque = Torque(t)/1000;

[Mathieu NOE]

hello

well simply replace T_a by Torque(t)/1000 in this line and you are done

T_a_fft = fft(detrend(T_a), NFFT)/L;  % detrend the signal before fft to remove the large DC (and near DC) fft output

is now :

T_a_fft = fft(detrend(Torque(t)/1000), NFFT)/L;  % detrend the signal before fft to remove the large DC (and near DC) fft output

here the full code updated :

A = load("torque_for_Sid_60s.mat");
T_a = A.torque_60s(1:5:600001);                  %Torque (Nm) chnages with time step
speed = 1000/60; 
tspan=0:0.00001*5:6;    %can try 4 or 5
y0 = [0;0;0;0;0;104.719;0;0;0;0;0;104.719/3]; 
tic
%These are constant for every call to reduced2 so compute them here and
%pass them into the modified function reduced3
time = 0:0.00001*5:6;
Torque = griddedInterpolant(time, T_a);
[t,y] = ode45(@(t,y) spur1(t,y,Torque),tspan,y0); 
TE= zeros(size(t));
KS = zeros(size(t));
for i=1:numel(t)
  [~,TE(i),KS(i)] = spur1(t(i),y(i,:),Torque);
end
% MN : remove first 0.1 second (transient)
idx = (t<0.1);
t(idx) = [];
TE(idx) = [];
T_a(idx) = [];
nexttile
plot(t,TE,"magenta")
xlabel('Time (s)')
ylabel('Transmission Error TE (meter)')
title('Transmission Error vs time')
axis padded
grid on
Fs = 1/(t(2)-t(1)); % Sampling frequency
L = length(t); % Length of signal
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
%%
TE_fft = fft(detrend(TE), NFFT)/L;  % detrend the signal before fft to remove the large DC (and near DC) fft output
f = Fs/2*linspace(0,1,NFFT/2+1);
nexttile;
plot(f, 2*abs(TE_fft(1:NFFT/2+1)))
xlabel('Frequency (Hz)')
ylabel('Amplitude')
title('FFT of TE vs frequency')
axis padded
grid on
%%
T_a_fft = fft(detrend(Torque(t)/1000), NFFT)/L;  % detrend the signal before fft to remove the large DC (and near DC) fft output
nexttile;
plot(f, 2*abs(T_a_fft(1:NFFT/2+1)))
xlabel('Frequency (Hz)')
ylabel('Amplitude')
title('FFT of Torque vs frequency')
axis padded
grid on
% nexttile
% plot(t,KS,"red")
% xlabel('Time (s)')
% ylabel('KS (N)')
% title('KS force vs time')
% axis padded
% grid on
newtime = toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [yp,TE,KS] = spur1(t,y,Torque)
yp = zeros(12,1);                           %vector of 12 equations
beta = 32*(pi/180);                      %Helix angle (rad) /13
P = 18.5*(pi/180);                         %Pressure angle (rad)
speed = 1000/60; 
Freq = 1000*20/60;
n_a = 22;                           % no of driver gear teeth  
n_p = 59;                           % no of driven gear teeth  
m = 2e-3*cos(beta);
R_a = (m*n_a)/2;                         %Radius
R_p = (m*n_p)/2;                         %Radius
i_r = n_p/n_a;                        % gear ratio
% Driver gear  
m_a = 5.3e-2;                      %mass of driver gear (kg)
c_ay=6.3e2;                       %Damping of driver gear in y direction (Ns/m)
c_az=1.3e2;                       %Damping of driver gear in z direction (Ns/m)
k_ay= 1.9e8;                      %Stiffness of driver gear in y direction (N/m)
k_az= 8.3e6;                      %Stiffness of driver gear in z direction (N/m)
I_a = 1.7e-5;          %Moment of inertia of driver gear (Kg.m3)
% Driven gear  
m_p = 2.5e-1;                                  %mass of driven gear (kg)
c_py=3.6e3;                                   %Damping of driven gear in y direction (Ns/m)
c_pz=4.3e2;                                   %Damping of driven gear in z direction (Ns/m)
k_py = 8.3e8;                                %Stiffness of driven gear in y direction (N/m)
k_pz = 1.2e7;                                  %Stiffness of driven gear in z direction (N/m)
I_p = 3.9e-4;                      %Moment of inertia of driven gear (Kg.m3)
Torque = Torque(t)/1000;
Opp_Torque = 68.853/1000;                        % Average torque value
% e = (pi*2*(R_a + R_p)*tan(beta))/(4*cos(P));
e= 20e-6 + 20e-8*sin(2*pi*Freq*t);  
ks = 0.9e8 + 0.9e6*sin(2*pi*Freq*t);                   %Shear stiffness 
% Cs = 0.1 + 0.0001*sin(2*pi*Freq*t);                     %Shear damping
Cs = 2*0.1*sqrt((ks*I_a*I_p)/(I_a*R_a + I_p*R_p)) + 2e-2*0.1*sqrt((ks*I_a*I_p)/(I_a*R_a + I_p*R_p))*sin(2*pi*Freq*t);
TE = -y(11)*R_p + y(5)*R_a ;                   %transmission error
b = 20e-6;                                       %Backlash
if(TE>b)
    h = TE - b;
        KS = h*ks;
elseif(TE < -1*b)
        h = TE + b;
        KS = h*ks;
else 
    h = 0;
   KS = h*ks; 
end
Fy = cos(beta)*(Cs*(cos(beta)*(y(2) - y(8) + R_a*y(6) - R_p*y(12)) + sin(beta)*(y(4) - y(10))) + ks*(cos(beta)*(y(1) - y(7) + R_a*y(5) - R_p*y(11)) + sin(beta)*(y(3) - y(9)))); %circumferential excitation force- formula from paper
Fz = sin(beta)*(Cs*(cos(beta)*(y(2) - y(8) + R_a*y(6) - R_p*y(12)) + sin(beta)*(y(4) - y(10))) + ks*(cos(beta)*(y(1) - y(7) + R_a*y(5) - R_p*y(11)) + sin(beta)*(y(3) - y(9)))); %circumferential excitation force- formula from paper
%Driver gear equations
yp(1) = y(2); 
yp(2) = (-Fy - k_ay*y(1) - c_ay*y(2))/m_a;  %acceleration in y (up and down) direction (m/s2)
yp(3) = y(4);  
yp(4) = (-Fz - k_az*y(3) - c_az*y(4))/m_a;  %acceleration in z (side to side) direction (m/s2)
yp(5) = y(6);  
yp(6) = (Torque - Fy*R_a)/I_a; %angular acceleration
%Driven gear equations
yp(7) = y(8);  
yp(8) = (Fy - k_py*y(7) - c_py*y(8))/m_p;            %acceleration in y (up and down) direction (m/s2)
yp(9) = y(10);  
yp(10) = (Fz - k_pz*y(9) - c_pz*y(10))/m_p;          %acceleration in y (side to side) direction (m/s2)
yp(11) = y(12);  
yp(12) = (-Opp_Torque*i_r + Fy*R_p)/I_p;                % angular accelration (rad/s2) 
end

[Siddharth Jain]

Hi again @Mathieu NOE,

Need your help again

Is it possible to make the linear scale FFT plot with the transient part in the beginning removed so that the plot is neater and individual frequency peaks can be seen.

[Mathieu NOE]

hello again

it's already done here

% MN : remove first 0.1 second (transient)
idx = (t<0.1);
t(idx) = [];
TE(idx) = [];
T_a(idx) = [];

[Mathieu NOE]

Hello

Problem solved ?

would you mind accepting my answer ? thanks !