
Waterfall diagram and fft for a vibration of an electric motor
104 views (last 30 days)
Show older comments
Eduardo
on 20 Jun 2025 at 13:12
Commented: Mathieu NOE
on 10 Jul 2025 at 14:56
Hello Everyone,
I have a data set: exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.mat
That consist of a test going slowyly from 0 to 3000 RPM of a motor, from which i want to do a waterfall diagramm to check for resonances, but i cant manage to do so. can somebody help?
I have an own code i tried to do but i dont hink its correct. here it is:
clear all
close all
% % % For WorkStation dSpace recorders!
file_name = 'exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated';
s=load(file_name);
%%% Select window type for FFT
% a=no window, b=Rectangular, c=Hann, d=Hamming, e=Flattop, f=Blackman-Harris, g=Nuttall, h=Chebyshev
winType = 'b';
length=length(s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.X(1).Data);
x=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.X(1).Data(1:length);
%%% Important base recording params
Fs=20000; %%% sampling frequency [Hz]
T=1/Fs;
%%% important base FFT params
fft_cycles = 2 ;
speed_aver_window = 1; %%% in sec
speed_aver_wind_points = speed_aver_window / T;
multipleORHz = 1; %%% 1 - multiply, 2 is Hz
%%% main indexis
index_speed_kistler = 1;
index_torque_an = 2;
index_torque_an_filt = 3;
index_vibr = 4;
index_vibr_filt = 5;
index_iq_act = 7;
index_id_act = 6;
index_speed_sw = 8;
index_speed_sew = 9;
%%% Define staret point of FFT:
time_start_fft =40;
point_start_fft = round(time_start_fft / T);
%%%% Indexes for analysis
index_main_speed = index_speed_kistler;
index_main_torque = index_torque_an_filt;
index_main_vibr = index_vibr;
speed_main = s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_main_speed).Data(1:length);
torque_main = s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_main_torque).Data(1:length);
vibr_main = s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_main_vibr).Data(1:length);
%%% Extraction of currents
Id_act_1=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_id_act).Data(1:length);
Iq_act_1=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_iq_act).Data(1:length);
%%% Extraction of speed and torque
SEW_speed=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_speed_sew).Data(1:length);
Kistler_speed=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_speed_kistler).Data(1:length);
Kistler_torque_an=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_torque_an).Data(1:length);
Kistler_torque_an_filt=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_torque_an_filt).Data(1:length);
N_pres=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_speed_sw).Data(1:length);
vibr=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_vibr).Data(1:length);
vibr_filt=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_vibr_filt).Data(1:length);
%%% plot results of the test
set(gcf,'color','white')
ax1=subplot(5,1,1);
plot(x,Kistler_torque_an,'r', x,Kistler_torque_an_filt,'b');
title('Torque Kistler');
xlabel('Time [s]')
ylabel('Torque [Nm]');
grid on
ax2=subplot(5,1,2);
plot(x,N_pres,'r',x,SEW_speed,'b',x,Kistler_speed,'g');
title('Speed');
xlabel('Time [s]')
ylabel('Speed [rpm]');
legend('Sensor','SEW','Kistler');
grid on
ax3=subplot(5,1,3);
plot(x,vibr,'.- r', x,vibr_filt,'b');
title('Vibrations');
xlabel('Time [s]')
ylabel('Acceleration [g]');
grid on
ax4=subplot(5,1,4);
plot(x,Id_act_1,'r');
title('Id act vs ref');
xlabel('Time [s]')
ylabel('Id [A]');
grid on
ax5=subplot(5,1,5);
plot(x,Iq_act_1,'r');
title('Iq act vs ref');
xlabel('Time [s]')
ylabel('Iq [A]');
grid on
linkaxes([ax1,ax2,ax3,ax4,ax5],'x');
%%% Making window to analyze FFT
speed_main_rpm_aver = abs(mean(speed_main(point_start_fft:(point_start_fft + speed_aver_wind_points))));
speed_rps = abs(speed_main_rpm_aver / 60);
period = 1 / speed_rps;
window_sec = fft_cycles * period;
window_poins = round(window_sec * Fs);
%%% Select and build desired window
switch lower(winType)
case 'a' % No window (raw data)
win = ones(window_poins,1);
case 'b' % Rectangular
win = rectwin(window_poins);
case 'c' % Hann
win = hann(window_poins);
case 'd' % Hamming
win = hamming(window_poins);
case 'e' % Flattop
win = flattopwin(window_poins);
case 'f' % Blackman-Harris
win = blackmanharris(window_poins);
case 'g' % Nuttall
win = nuttallwin(window_poins);
case 'h' % Chebyshev
win = chebwin(window_poins, 60); % 60 dB sidelobe suppression
otherwise
error('Unknown winType "%s"', winType);
end
% Compute mean values for DC removal
torque_main_aver = mean(torque_main(point_start_fft:(point_start_fft + speed_aver_wind_points)));
vibr_main_aver = mean(vibr_main(point_start_fft:(point_start_fft + speed_aver_wind_points)));
% Preallocate arrays
time_wind = zeros(window_poins, 1);
torque_wind = zeros(window_poins, 1);
vibr_wind = zeros(window_poins, 1);
speed_wind = zeros(window_poins, 1);
% Create windowed signals
j = 1;
for i = point_start_fft:(point_start_fft + window_poins - 1)
time_wind(j) = x(i);
torque_wind(j) = (torque_main(i) - torque_main_aver) * win(j);
vibr_wind(j) = (vibr_main(i) - vibr_main_aver) * win(j);
speed_wind(j) = speed_main(i);
j = j + 1;
end
figure
set(gcf,'color','white')
bx1=subplot(3,1,1);
plot(time_wind,torque_wind,'r');
title('Torque window');
xlabel('Time [s]')
ylabel('Torque [Nm]');
grid on
bx2=subplot(3,1,2);
plot(time_wind,vibr_wind,'r');
title('Acceleretion window');
xlabel('Time [s]')
ylabel('Acceleration [g]');
grid on
bx3=subplot(3,1,3);
plot(time_wind,speed_wind,'r');
title('Speed window');
xlabel('Time [s]')
ylabel('Speed [rpm]');
grid on
linkaxes([bx1,bx2,bx3],'x');
%%% FFT of Vibrations
L=window_poins;
t = (0:L-1)*T; % Time vector
Y1=fft(vibr_wind,L);
P2 = abs(Y1/L).^2;
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
if multipleORHz == 1
freq_ref = speed_main_rpm_aver/60;
freq_plot_lim = 100;
elseif multipleORHz == 2
freq_ref = 1;
freq_plot_lim = 4000;
end
f1 = (Fs*(0:(L/2))/L)/(freq_ref);
figure
set(gcf,'color','white')
bar(f1,P1)
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (multiple of mechanical frequency)')
xlim([0 freq_plot_lim])
ylabel('Absolute value of Harmonic VIBRATIONS [g]')
% --- Define cutoff RPM for data selection ---
rpm_min = 0;
rpm_max = 4200;
% Find indices corresponding to rpm_min and rpm_max in time vector x
ind_min = find(speed_main >= rpm_min, 1, 'first');
ind_max = find(speed_main <= rpm_max, 1, 'last');
% Cut signals based on these indices
time_cut = x(ind_min:ind_max);
rpm_cut = speed_main(ind_min:ind_max);
vibr_cut = vibr_main(ind_min:ind_max);
if any(rpm_cut <= 0)
warning('RPM data contains non-positive values. Fixing...');
rpm_cut(rpm_cut <= 0) = NaN; % Or interpolate, or remove
% Option 1: interpolate missing values (if feasible)
rpm_cut = fillmissing(rpm_cut, 'linear');
% Option 2: truncate all rows with NaNs (if that's acceptable)
validIdx = ~isnan(rpm_cut);
vibr_cut = vibr_cut(validIdx);
rpm_cut = rpm_cut(validIdx);
end
% Sampling frequency from original code
Fs = 4000; % Hz
% % ---- Generate order map and waterfall plot ----
% Assuming you have the rpmordermap function in your path.
% If not, I can help to replace it with an equivalent.
% Frequency-based RPM map
[map, freq, rpm_axis, time_map, res] = rpmordermap(vibr_cut, Fs, rpm_cut, 2, ...
'Scale', 'dB', 'Window', 'hann', 'Amplitude', 'rms');
[fr, rp] = meshgrid(freq, rpm_axis);
% Waterfall plot in frequency
figure;
waterfall(fr, rp, map');
view(6, -60);
grid on;
xlabel('Frequency [Hz]');
ylabel('RPM');
zlabel('Amplitude [dB]');
title('Waterfall Plot (Frequency Map)');
% --- FFT over speed steps (similar to your reference code) ---
% Define speed steps (adjust based on your rpm range)
speed_steps = rpm_min:100:rpm_max;
% Define FFT window length in seconds (e.g., 0.5 s)
fft_window_sec = 0.5;
fft_window_points = round(fft_window_sec * Fs);
% Initialize figure for FFT results
for i = 1:length(speed_steps)
% Find start time index closest to current speed step
idx_start = find(rpm_cut >= speed_steps(i), 1, 'first');
if isempty(idx_start) || (idx_start + fft_window_points - 1) > length(vibr_cut)
continue; % Skip if index invalid or window exceeds data length
end
% Extract segment
segment = vibr_cut(idx_start : idx_start + fft_window_points - 1);
time_segment = time_cut(idx_start : idx_start + fft_window_points - 1);
% Remove DC
segment = segment - mean(segment);
% Apply Hann window
winvec = hann(length(segment));
% FFT
L = length(segment);
Y = fft(segment .* winvec);
P2 = abs(Y / L);
P1 = P2(1 : floor(L/2) + 1);
P1(2:end-1) = 2 * P1(2:end-1);
f = Fs * (0:(L/2)) / L;
% Plot time-domain and FFT for this step
subplot(2,1,1);
plot(time_segment, segment);
title(sprintf('Vibration at Speed Step: %d RPM', speed_steps(i)));
xlabel('Time [s]');
ylabel('Acceleration [g]');
grid on;
hold on;
subplot(2,1,2);
bar(f, P1);
title(sprintf('FFT Spectrum at Speed Step: %d RPM', speed_steps(i)));
xlabel('Frequency [Hz]');
ylabel('Amplitude');
grid on;
hold on;
end
code should be correct until the part of the water fall.
also the index of the data is no longer 1 to 9 but 1 to 3 being speed-raw vibration- filtered vibration
I would appreciete if somebody could help me make a code i can use with different data sets, maybe some have 1 to 6 and with the option of ploting the sigbals at the beginning and the fft how in the code.
Thanks in advanced
0 Comments
Accepted Answer
Mathieu NOE
on 20 Jun 2025 at 15:09
hello again
this is a demo based on your previously posted data (could not load the one in this post)
spectrogram + time plot (with aligned time axes)

code :
%% For WorkStation dSpace recorders!
file_name = 'exp1_HV_eaux_B0_ASQ_N7_vibrations_NT3a_pos_3000rpm_20kHz_foam';
s=load(file_name);
%% use this trick to avoid using the long filename every time
StrName=fieldnames(s);
StrName=StrName{1};
%% main Code
time = s.(StrName).X.Data;
duration = time(end) - time(1);
dt = mean(diff(time));
Fs = 1/dt;
%% main indexes
index_speed_kistler = 1;
index_torque_an = 2;
index_torque_an_filt = 3;
index_vibr = 4;
index_vibr_filt = 5;
index_iq_act = 7;
index_id_act = 6;
index_speed_sw = 8;
index_speed_sew = 9;
%%%% Indexes for analysis
index_main_speed = index_speed_sw;
index_main_torque = index_torque_an_filt;
index_main_vibr = index_vibr;
vibr_main = s.(StrName).Y(index_main_vibr).Data;
%% Extraction of speed and torque
SEW_speed=s.(StrName).Y(index_speed_sw).Data;
Kistler_speed=s.(StrName).Y(index_speed_kistler).Data;
% extract time portion
time_start_fft = 0;
% ind = (time>=time_start_fft) & (time<=time_start_fft+duration);
ind = (SEW_speed>=1);
time = time(ind);
vibr_main = vibr_main(ind);
SEW_speed = SEW_speed(ind);
Kistler_speed = Kistler_speed(ind);
%% FFT plot
nfft = 1024*8; %
Overlap = 0.75; % (75% ) overlap
df = Fs/nfft;
% check signal stationnarity with spectrogram : OK
nfft = 1024; %
Overlap = 0.75; % (75% ) overlap
[S,F,T] = spectrogram(vibr_main,hanning(nfft),round(Overlap*nfft),nfft,Fs,'yaxis','minthreshold',-80); % WINDOW,NOVERLAP,NFFT,Fs
SEW_speed_interp = interp1(time,SEW_speed,T);
% freq_normalized = F*60/SEW_speed_interp;
dBrange = 80; % dB range (anything below max - dBrange will not be displayed)
S_dB= 20*log10(abs(S));
figure;
ax = gobjects(2,1);
ax(1) = subplot(2,1,1);
imagesc(T,F,S_dB);
xlabel('Time (s)');
ylabel('Frequency (Hz)');
title('Spectrogram');
colormap('jet')
hcb=colorbar('ver'); % colorbar handle
hcb.FontSize = 9;
hcb.Title.String = "dB";
hcb.Title.FontSize = 15;
set(gca,'YDir','normal');
xlim([T(1) T(end)]);
ylim([0 5000]);
caxis([max(S_dB(:))-dBrange max(S_dB(:))]);
ax(2) = subplot(2,1,2);
plot(ax(end), time,vibr_main)
grid on
xlim([T(1) T(end)]);
% Set the width and height of all axes to the min width & height
allAxPos = vertcat(ax.Position);
allAxPos(:,3:4) = min(allAxPos(:,3:4)).*ones(numel(ax),1);
set(ax,{'position'},mat2cell(allAxPos,ones(numel(ax),1),4))
11 Comments
William Rose
on 10 Jul 2025 at 14:28
@Eduardo, I accepted the answer from @Mathieu NOE, on your behalf, since it seems the problems are solved as much as they can be. It is a courtesy to accept an answer, especially when the answerer provides so much good advice.
Mathieu NOE
on 10 Jul 2025 at 14:56
I truly appreciate
More Answers (1)
Eduardo
on 23 Jun 2025 at 8:45
11 Comments
Mathieu NOE
on 4 Jul 2025 at 13:33
glad to hear that you have (almost) achieved your goal
yeap EMI can cause a lot of trouble
proper grounding , shielding small signals wires , differential inputs, it can take a while until you get clean signals inside your PC
good luck for the future and always remember to ask yourself if you can trust your screen .
better double check than go blind
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!