
Waterfall diagram and fft for a vibration of an electric motor
    29 visualizaciones (últimos 30 días)
  
       Mostrar comentarios más antiguos
    
    Eduardo
 el 20 de Jun. de 2025
  
    
    
    
    
    Comentada: Mathieu NOE
      
 el 10 de Jul. de 2025
            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 comentarios
Respuesta aceptada
  Mathieu NOE
      
 el 20 de Jun. de 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 comentarios
  William Rose
      
 el 10 de Jul. de 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. 
  Mathieu NOE
      
 el 10 de Jul. de 2025
				I truly appreciate
Más respuestas (1)
  Eduardo
 el 23 de Jun. de 2025
        11 comentarios
  Mathieu NOE
      
 el 4 de Jul. de 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
Ver también
Categorías
				Más información sobre Vibration Analysis en Help Center y File Exchange.
			
	Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!






