
Waterfall diagram and fft for a vibration of an electric motor
    28 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!









