Lalo2202
- 2
- 0
- TL;DR Summary
- 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
Im usinng matlab 2021[/s][/s][/s][/s]
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
Code:
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
Last edited by a moderator: