what is wrong with my FFT from a Shaker with accelerometer mounted

1 view (last 30 days)
I tried to compute the FFT of the data attached to no avail using the FFT code below. The sensor was mounted on a shaker and I applied some frequencies from the generator at 200Hz , 100Hz and 50Hz. when I plotted the FFT I could not get the frequencies applied. from the data the first column is the timestamps and the second column is the z-data acceleration.
% Read the data from the file
z_data = load('SensorData_fromxiao1.txt');
% Separate time and acceleration data
time_data = z_data(:, 1);
accel_data = z_data(:, 2);
% Remove mean to eliminate DC offset
accel_data = accel_data - mean(accel_data);
Fs = 98.3; % Sampling frequency in Hz, based on your calculation
% Apply windowing to the data
windowed_data = accel_data .* hann(length(accel_data));
% Perform FFT
L = length(windowed_data); % Length of the signal
Y = fft(windowed_data); % Perform FFT
P2 = abs(Y/L); % Two-sided spectrum
P1 = P2(1:L/2+1); % Single-sided spectrum
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L; % Frequency vector
% Plot FFT result
figure;
plot(f, P1);
title('Single-sided amplitude spectrum');
xlabel('Frequency (Hz)');
ylabel('Amplitude');

Answers (1)

Star Strider
Star Strider on 25 Aug 2023
Edited: Star Strider on 25 Aug 2023
I compared your fft code with my stock fft code and got the same result.
I believe that the problem is the sampling frequency of your instrumentation. The sampling frequency of the data is about 98 Hz, providing a Nyquist frequency of 49 Hz. This is simply too low to see specific frequenciy peaks any higher than that, so assuming that your instrumentation did not have an anti-aliasing filter, we can assume that the one peak at about 5 Hz is aliased.
The problem is not your fft code, the problem is that to see a frequency of 200 Hz (and by extension the other frequencies below it), the sampling frequency should have been at least 500 Hz. If your instrumentation has that capability, it would be best to re-do the experiment with a 500 Hz (or higher) sampling frequency.
The result —
% Read the data from the file
z_data = load('SensorData_fromxiao1.txt');
% Separate time and acceleration data
time_data = z_data(:, 1);
accel_data = z_data(:, 2);
% Remove mean to eliminate DC offset
accel_data = accel_data - mean(accel_data);
Fs = 98.3; % Sampling frequency in Hz, based on your calculation
% Apply windowing to the data
windowed_data = accel_data .* hann(length(accel_data));
% Perform FFT
L = length(windowed_data); % Length of the signal
Y = fft(windowed_data); % Perform FFT
P2 = abs(Y/L); % Two-sided spectrum
P1 = P2(1:fix(L/2)+1); % Single-sided spectrum
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L; % Frequency vector
% Plot FFT result
figure;
plot(f, P1);
title('Single-sided amplitude spectrum');
xlabel('Frequency (Hz)');
ylabel('Amplitude');
t = time_data;
s = accel_data;
% dt = diff(t);
% sts = [mean(dt) std(dt)]
Fs = round(1/mean(diff(t)))
Fs = 98
[sr,tr] = resample(s,t,Fs);
[FTs1,Fv] = FFT1(sr,tr);
figure
plot(Fv, abs(FTs1))
grid
% Illustrating the effects of subsampling on aliasing —
Fs = 1000; % Simulate Signal
TL = 120; % Simulate Signal
t = linspace(0, TL*Fs, TL*Fs+1)/Fs; % Simulate Signal
s = sum(sin([50; 100; 150] * 2*pi*t)); % Actual Signal
sc = s(1:2:end); % Correctly Sampled Signal
tc = t(1:2:end); % Correctly Sampled Signal
Fsc = 1/(tc(2)-tc(1))
Fsc = 500
ss = s(1:10:end); % Subsampled Signal
ts = t(1:10:end); % Subsampled Signal
Fss = 1/(ts(2)-ts(1))
Fss = 100
[FTs1,Fv] = FFT1(s,t);
figure
plot(Fv, abs(FTs1))
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Original Signal')
[FTsc1,Fvc] = FFT1(sc,tc);
figure
plot(Fvc, abs(FTsc1))
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Correctly Sampled Signal')
[FTss1,Fvs] = FFT1(ss,ts);
figure
plot(Fvs, abs(FTss1))
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Incorrectly Sampled Signal')
function [FTs1,Fv] = FFT1(s,t)
s = s(:);
t = t(:);
L = numel(t);
Fs = 1/mean(diff(t));
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTs = fft((s - mean(s)).*hann(L), NFFT)/sum(hann(L));
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
FTs1 = FTs(Iv);
end
EDIT — (25 Aug 2023 at 23:03)
Added the secton illustrating the effect of subsampling and aliasing.
.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!