Here is how to perform FFT and compute the phase angle values in rad:
data = readmatrix('D_Sample.xlsx'); % Sample data import
t = data(:,1); % Column 1 is time data
N = numel(t); % Number of data points
Ts=mean(diff(t)); % Data Sampling time, [s]
Fs=1/Ts; % Data Sampling frequency, [Hz]
% Measured data has some noise that is cleaned using low-pass filter
xf =lowpass(data(:,2), 175,Fs) ;
yf =lowpass(data(:,2), 175,Fs) ;
zf =lowpass(data(:,2), 175,Fs) ;
Nfft=2^nextpow2(N);
x=fft(xf,Nfft)/N;
y=fft(yf,Nfft)/N;
z=fft(zf,Nfft)/N;
f=Fs/2*linspace(0,1,Nfft/2+1);
figure(2)
plot(f,2*abs(y(1:Nfft/2+1)));
grid on
title('Single-Sided Amplitude Spectrum of acceleration: a_x(t)');
xlabel('Frequency(Hz)');
ylabel('|a_x(f)|');
phasex=angle(x);
phasey=angle(y);
phasez=angle(z);
figure(2)
subplot(2,1,1)
plot(f,2*abs(x(1:Nfft/2+1)));
grid on
title('Single-Sided Spectrum of a_x(t)');
ylabel('|a_x(f)|');
subplot(2,1,2)
plot(f,phasex(1:Nfft/2+1))
grid on
xlabel('Frequency(Hz)');
ylabel('Phase of a_x(rad)')
figure(3)
subplot(2,1,1)
plot(f,2*abs(y(1:Nfft/2+1)));
grid on
title('Single-Sided Spectrum of a_y(t)');
ylabel('|a_y(f)|');
subplot(2,1,2)
plot(f,phasey(1:Nfft/2+1))
grid on
xlabel('Frequency(Hz)');
ylabel('Phase of a_y(rad)')
figure(4)
subplot(2,1,1)
plot(f,2*abs(z(1:Nfft/2+1)));
grid on
title('Single-Sided Spectrum of a_z(t)');
ylabel('|a_z(f)|');
subplot(2,1,2)
plot(f,phasez(1:Nfft/2+1))
grid on
xlabel('Frequency(Hz)');
ylabel('Phase of a_z(rad)')
