
Waterfall diagram and fft for a vibration of an electric motor
26 views (last 30 days)
Show older comments
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
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
@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.
More Answers (1)
Eduardo
on 23 Jun 2025
11 Comments
Mathieu NOE
on 4 Jul 2025
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
Categories
Find more on Vibration Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!





