you can use the following code which uses fft() function
fft_voltage = fft(voltage);
power_spectrum = abs(fft_voltage).^2;
dt = mean(diff(timestamps)); % time resolution
df = 1 / (length(voltage) * dt); % frequency resolution
f = (0:length(voltage)-1) * df; % frequency axis
plot(f, power_spectrum);
xlabel('Frequency (Hz)');
ylabel('PSD');
You can read about this in this documentation Power Spectral Density Estimates Using FFT