Waterfall plot of an PMSM

[Lalo2202]

TL;DR

Analysis of vibrations over the full RPM range of a PMSM

Hello everyone,

Im working on analizing the spectram domain of frequencies present on an PMSM when turning alone.

I attach an image of my Waterfall diagramm results, as you can see there are several frequencies that appear in parallel pairs with the same amplitude along the whole RPM range, my assumption is that these pairs represent a single frequency located in the middle but for some reason it is being represented as in the imaged be solved in the script?

or is there a problem in my PMSM that could make my frequencies look like this. search online for possible reasons why this happens

clear all
close all

% Load data
file_name = 'exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated';
s=load(file_name);

% FFT window type: a=no window, b=Rectangular, c=Hann, d=Hamming, e=Flattop, f=Blackman-Harris, g=Nuttall, h=Chebyshev
winType = 'b'; 

length_data = 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_data);

Fs = 20000; T = 1/Fs;

fft_cycles = 3;
speed_aver_window = 1; %samller better for non constant speed
speed_aver_wind_points = speed_aver_window / T;
multipleORHz = 1;  % 1 = order, 2 = Hz

index_speed_sw = 1;
index_vibr = 2;
index_vibr_filt = 3;

speed_main = s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_speed_sw).Data(1:length_data);
vibr_main  = s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_vibr_filt).Data(1:length_data);
vibr_filt  = s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_vibr_filt).Data(1:length_data);


set(gcf,'color','white')
ax1=subplot(5,1,2);
plot(x,speed_main,'r');
title('Speed'); xlabel('Time [s]'); ylabel('Speed [rpm]'); legend('Sensor'); grid on

ax2=subplot(5,1,3);
plot(x,vibr_main,'.- r', x,vibr_filt,'b');
title('Vibrations'); xlabel('Time [s]'); ylabel('Acceleration [g]'); grid on

linkaxes([ax1,ax2],'x');

time_start_fft = 61;
point_start_fft = round(time_start_fft / T);

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);

% Window function
switch lower(winType)
    case 'a', win = ones(window_poins,1);
    case 'b', win = rectwin(window_poins);
    case 'c', win = hann(window_poins);
    case 'd', win = hamming(window_poins);
    case 'e', win = flattopwin(window_poins);
    case 'f', win = blackmanharris(window_poins);
    case 'g', win = nuttallwin(window_poins);
    case 'h', win = chebwin(window_poins, 60);
    otherwise, error('Unknown winType "%s"', winType);
end

vibr_main_aver = mean(vibr_main(point_start_fft:(point_start_fft + speed_aver_wind_points)));

time_wind = zeros(window_poins, 1);
vibr_wind = zeros(window_poins, 1);
speed_wind = zeros(window_poins, 1);

j = 1;
for i = point_start_fft:(point_start_fft + window_poins - 1)
    time_wind(j) = x(i);
    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,2);
plot(time_wind,vibr_wind,'r');
title('Acceleretion window'); xlabel('Time [s]'); ylabel('Acceleration [g]'); grid on

bx2=subplot(3,1,3);
plot(time_wind,speed_wind,'r');
title('Speed window'); xlabel('Time [s]'); ylabel('Speed [rpm]'); grid on
linkaxes([bx1,bx2],'x');

Y1=fft(vibr_wind);
L=window_poins;
P2 = abs(Y1/L);
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]')

%% === Waterfall Plot (Frequency Domain 0–5000 Hz, up to 3000 RPM) ===
Fs = 20000; % already defined earlier
rpm = speed_main;         % speed in RPM
vib = vibr_main;          % raw vibration signal

% Limit to 0–3000 RPM and exclude non-positive values
rpm_mask = rpm > 0 & rpm <= 3000;

vib_cut = vib(rpm_mask);
rpm_cut = rpm(rpm_mask);
x_cut = x(rpm_mask); % optional time axis

% Create frequency map
[map_f, freq, rpm_f, time_f, res_f] = rpmfreqmap(vib_cut, Fs, rpm_cut, 3, ...
    'Scale', 'dB', 'Window', 'hamming', 'Amplitude', 'rms');

% Limit to frequencies ≤ 2000 Hz
freq_limit = 2000;
freq_idx = freq <= freq_limit;

map_f = map_f(freq_idx, :);
freq = freq(freq_idx);


% Plot Waterfall (Frequency)
[fr, rp] = meshgrid(freq, rpm_f);
figure;
waterfall(fr, rp, map_f');
view(6, -60);
grid on;
xlabel('Frequency [Hz]');
ylabel('RPM');
zlabel('Amplitude [dB]');
title('Waterfall Plot (Frequency Domain)');
xlim([0 2000]);
ylim([0 3000]);

% Define which orders you want to overlay
ordersToPlot = [1, 2, 3];  % Adjust as needed
colors = lines(length(ordersToPlot));  % Different colors for each line

hold on
for k = 1:length(ordersToPlot)
    % Calculate the frequency line for that order (f = order * RPM / 60)
    freq_line = ordersToPlot(k) * rpm_f / 60;
    
    % Skip frequencies beyond x-axis limit (e.g., 2000 Hz)
    freq_line(freq_line > 2000) = NaN;
    
    % z-position: pick a value slightly above the noise floor
    z_line = ones(size(freq_line)) * (-20 + 5 * k);  % Adjust height if needed

    % Plot the order line
    plot3(freq_line, rpm_f, z_line, '--', 'LineWidth', 2, 'Color', colors(k,:));
    
    % Optionally, add a label
    text(freq_line(end), rpm_f(end), z_line(end), ...
        sprintf('%d×', ordersToPlot(k)), ...
        'Color', colors(k,:), 'FontWeight', 'bold', 'FontSize', 10);
end
hold off

Im usinng matlab 2021[/s][/s][/s][/s]

 

[Lalo2202]

For information the biggest 2 peaks are around 9.5 and 11.5 order, that links to the 2nd electrical frequencie of my motor with 10 poles