Hello,
I’m trying to determine the index for a change in gradient along the radial direction from x=0 to x=1, which is expected to be negative. My code which runs but doesn’t give correct result. Kindly have a look at the same and suggest me the possible correction(s) please.
data.variable.gradpressure = 10000 x 100 matrix
data.variable.x = 1 x 100 vector
Moreover the matchflag is also false when I try to match the columns which are both equal.
The code is as below:
% Determine the index for negative gradient from x=0 to x=1
clear grad_width;
global grad_width;
global data;

match_flags = false(size(data.variable.gradpressure,1),1); % Logical vector for matching rows

% Check if any row of pressure equals x
for i = 1:size(data.variable.x, 1)
if isequal(data.variable.gradpressure(:,i)’, data.variable.x(:,i))
match_flags(i) = true;
end
end
%disp(match_flags)

threshold_fraction = 0.1;

% Preallocate result arrays
nTimes = length(data.variable.t);
grad_foot_idx = zeros(nTimes,1);
grad_head_idx = zeros(nTimes,1);
grad_width = zeros(nTimes,1);

for i = 1:nTimes
gradp = data.variable.gradpressure(i,:); % 1D vector at time i

% Avoid division by zero
valid = gradp(1:end-1) ~= 0;

% Compute gradient change (like second derivative)
change_gp = zeros(1, length(gradp)-1);
change_gp(valid) = diff(gradp)./gradp(1:end-1); % Central diff approx

% Find steepest drop (most negative change)
[max_grad1, max_grad1_idx] = min(abs(change_gp)); % Min of |second derivative|
grad_threshold = threshold_fraction * abs(max_grad1);

% Search left (foot)
left_idx = max_grad1_idx;
while left_idx > 1 && abs(change_gp(left_idx)) > grad_threshold
left_idx = left_idx – 1;
end

% Search right (head)
right_idx = max_grad1_idx;
while right_idx < length(change_gp) && abs(change_gp(right_idx)) > grad_threshold
right_idx = right_idx + 1;
end

% Store
grad_foot_idx(i) = left_idx;
grad_head_idx(i) = right_idx;
grad_width(i) = data.variable.x(right_idx) – data.variable.x(left_idx);
end

% Display last result (or modify to plot or analyze all)
disp(grad_width(end));Hello,
I’m trying to determine the index for a change in gradient along the radial direction from x=0 to x=1, which is expected to be negative. My code which runs but doesn’t give correct result. Kindly have a look at the same and suggest me the possible correction(s) please.
data.variable.gradpressure = 10000 x 100 matrix
data.variable.x = 1 x 100 vector
Moreover the matchflag is also false when I try to match the columns which are both equal.
The code is as below:
% Determine the index for negative gradient from x=0 to x=1
clear grad_width;
global grad_width;
global data;

match_flags = false(size(data.variable.gradpressure,1),1); % Logical vector for matching rows

% Check if any row of pressure equals x
for i = 1:size(data.variable.x, 1)
if isequal(data.variable.gradpressure(:,i)’, data.variable.x(:,i))
match_flags(i) = true;
end
end
%disp(match_flags)

threshold_fraction = 0.1;

% Preallocate result arrays
nTimes = length(data.variable.t);
grad_foot_idx = zeros(nTimes,1);
grad_head_idx = zeros(nTimes,1);
grad_width = zeros(nTimes,1);

for i = 1:nTimes
gradp = data.variable.gradpressure(i,:); % 1D vector at time i

% Avoid division by zero
valid = gradp(1:end-1) ~= 0;

% Compute gradient change (like second derivative)
change_gp = zeros(1, length(gradp)-1);
change_gp(valid) = diff(gradp)./gradp(1:end-1); % Central diff approx

% Find steepest drop (most negative change)
[max_grad1, max_grad1_idx] = min(abs(change_gp)); % Min of |second derivative|
grad_threshold = threshold_fraction * abs(max_grad1);

% Search left (foot)
left_idx = max_grad1_idx;
while left_idx > 1 && abs(change_gp(left_idx)) > grad_threshold
left_idx = left_idx – 1;
end

% Search right (head)
right_idx = max_grad1_idx;
while right_idx < length(change_gp) && abs(change_gp(right_idx)) > grad_threshold
right_idx = right_idx + 1;
end

% Store
grad_foot_idx(i) = left_idx;
grad_head_idx(i) = right_idx;
grad_width(i) = data.variable.x(right_idx) – data.variable.x(left_idx);
end

% Display last result (or modify to plot or analyze all)
disp(grad_width(end)); Hello,
I’m trying to determine the index for a change in gradient along the radial direction from x=0 to x=1, which is expected to be negative. My code which runs but doesn’t give correct result. Kindly have a look at the same and suggest me the possible correction(s) please.
data.variable.gradpressure = 10000 x 100 matrix
data.variable.x = 1 x 100 vector
Moreover the matchflag is also false when I try to match the columns which are both equal.
The code is as below:
% Determine the index for negative gradient from x=0 to x=1
clear grad_width;
global grad_width;
global data;

match_flags = false(size(data.variable.gradpressure,1),1); % Logical vector for matching rows

% Check if any row of pressure equals x
for i = 1:size(data.variable.x, 1)
if isequal(data.variable.gradpressure(:,i)’, data.variable.x(:,i))
match_flags(i) = true;
end
end
%disp(match_flags)

threshold_fraction = 0.1;

% Preallocate result arrays
nTimes = length(data.variable.t);
grad_foot_idx = zeros(nTimes,1);
grad_head_idx = zeros(nTimes,1);
grad_width = zeros(nTimes,1);

for i = 1:nTimes
gradp = data.variable.gradpressure(i,:); % 1D vector at time i

% Avoid division by zero
valid = gradp(1:end-1) ~= 0;

% Compute gradient change (like second derivative)
change_gp = zeros(1, length(gradp)-1);
change_gp(valid) = diff(gradp)./gradp(1:end-1); % Central diff approx

% Find steepest drop (most negative change)
[max_grad1, max_grad1_idx] = min(abs(change_gp)); % Min of |second derivative|
grad_threshold = threshold_fraction * abs(max_grad1);

% Search left (foot)
left_idx = max_grad1_idx;
while left_idx > 1 && abs(change_gp(left_idx)) > grad_threshold
left_idx = left_idx – 1;
end

% Search right (head)
right_idx = max_grad1_idx;
while right_idx < length(change_gp) && abs(change_gp(right_idx)) > grad_threshold
right_idx = right_idx + 1;
end

% Store
grad_foot_idx(i) = left_idx;
grad_head_idx(i) = right_idx;
grad_width(i) = data.variable.x(right_idx) – data.variable.x(left_idx);
end

% Display last result (or modify to plot or analyze all)
disp(grad_width(end)); gradient MATLAB Answers — New Questions

​

Hello
the solution of the PDE:
df(x,y)/dx + df(x,y)/dy =0
should be f= f(x-y) , with f an arbitrary function of the argument which is the compound variable (x-y).
Wanting to verify this symbolic solution I would like to execute these commands
1) syms x y real
2) syms f(x-y)
and then do the calculation
3) ver=diff(f(x-y),x)+diff(f(x-y),y)
and simplifying it should give
ver=zero

but unfortunately on the second line it
gives me an error :
———————————————————————
Error using symfun.parseString (line 101)
Invalid variable name.

Error in syms (line 276)
[name, vars] = symfun.parseString(x,’allowPercent’);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
———————————————————————–

I was kindly asking
how to proceed and be able to do this
symbolic verification

thankyou very much !
ValerioHello
the solution of the PDE:
df(x,y)/dx + df(x,y)/dy =0
should be f= f(x-y) , with f an arbitrary function of the argument which is the compound variable (x-y).
Wanting to verify this symbolic solution I would like to execute these commands
1) syms x y real
2) syms f(x-y)
and then do the calculation
3) ver=diff(f(x-y),x)+diff(f(x-y),y)
and simplifying it should give
ver=zero

but unfortunately on the second line it
gives me an error :
———————————————————————
Error using symfun.parseString (line 101)
Invalid variable name.

Error in syms (line 276)
[name, vars] = symfun.parseString(x,’allowPercent’);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
———————————————————————–

I was kindly asking
how to proceed and be able to do this
symbolic verification

thankyou very much !
Valerio Hello
the solution of the PDE:
df(x,y)/dx + df(x,y)/dy =0
should be f= f(x-y) , with f an arbitrary function of the argument which is the compound variable (x-y).
Wanting to verify this symbolic solution I would like to execute these commands
1) syms x y real
2) syms f(x-y)
and then do the calculation
3) ver=diff(f(x-y),x)+diff(f(x-y),y)
and simplifying it should give
ver=zero

but unfortunately on the second line it
gives me an error :
———————————————————————
Error using symfun.parseString (line 101)
Invalid variable name.

Error in syms (line 276)
[name, vars] = symfun.parseString(x,’allowPercent’);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
———————————————————————–

I was kindly asking
how to proceed and be able to do this
symbolic verification

thankyou very much !
Valerio symbolic functions MATLAB Answers — New Questions

​

I did a simulink simulation for PV array with DC/DC Buck Converter to charged the battery. But i after try to simulate it i found that my PV voltage is close to my Output Voltage. It seems like its not a correct results. The parameter of PV Voc (19V), Isc (0.124A) ,Vmp (15.4V), Imp ( 0.1A) is set by my ownself from a solar panel datasheet. While the converter should step down the voltage to 8.4V (2s lithium ion battery). P&O MPPT already been implement in this circuit.I did a simulink simulation for PV array with DC/DC Buck Converter to charged the battery. But i after try to simulate it i found that my PV voltage is close to my Output Voltage. It seems like its not a correct results. The parameter of PV Voc (19V), Isc (0.124A) ,Vmp (15.4V), Imp ( 0.1A) is set by my ownself from a solar panel datasheet. While the converter should step down the voltage to 8.4V (2s lithium ion battery). P&O MPPT already been implement in this circuit. I did a simulink simulation for PV array with DC/DC Buck Converter to charged the battery. But i after try to simulate it i found that my PV voltage is close to my Output Voltage. It seems like its not a correct results. The parameter of PV Voc (19V), Isc (0.124A) ,Vmp (15.4V), Imp ( 0.1A) is set by my ownself from a solar panel datasheet. While the converter should step down the voltage to 8.4V (2s lithium ion battery). P&O MPPT already been implement in this circuit. mppt, dc/dc converter, pv arrays, buck converter, maximum power point tracking, charging, battery MATLAB Answers — New Questions

​

I am working on a CI pipeline using a headless version of MATLAB in a Docker® container. Building Simulink Real-Time models requires the "Simulink Real-Time Target Support Package" and "Speedgoat I/O Blockset" (mandatory since R2024a).
However, it seems that these support packages cannot be installed inside a Docker container. How can I get this CI pipeline up and running?I am working on a CI pipeline using a headless version of MATLAB in a Docker® container. Building Simulink Real-Time models requires the "Simulink Real-Time Target Support Package" and "Speedgoat I/O Blockset" (mandatory since R2024a).
However, it seems that these support packages cannot be installed inside a Docker container. How can I get this CI pipeline up and running? I am working on a CI pipeline using a headless version of MATLAB in a Docker® container. Building Simulink Real-Time models requires the "Simulink Real-Time Target Support Package" and "Speedgoat I/O Blockset" (mandatory since R2024a).
However, it seems that these support packages cannot be installed inside a Docker container. How can I get this CI pipeline up and running?  MATLAB Answers — New Questions

​

Hello, I have some data where I my relevant column (Y position) is as this yellow graph

I want to pull out just the data that is between the green arrows – I do know the actual starting and finishing Y positions (101 and 121 respectively) , but there are other values at these values – but they are constant

My data is actually in an array, and the Yposition is column 2, but I delete the 1st column, so the column becomes 1
this was my attempt to remove all rows of my array outside of these green arrows.

% Get user input to the real starting and finishing Y values
prompt = {‘Col:’,’start:’,’finish:’};
dlgtitle = ‘Keep Rows Between’;
dims = [1 35];
definput = {‘1′,’101′,’121’}; % definput = {‘Region:’,num2str(X),num2str(Y)};
answer = inputdlg(prompt,dlgtitle,dims,definput);
col=str2num(answer{1}); low=str2num(answer{2}); high=str2num(answer{3});

% Get Current table data
t=app.UITable;
d=t.Data;

d=rmmissing(d); % One wat Remove Nans

Y=table2array(d);
Y=Y(:,col);

low
head(Y)
class(Y)

% Get index of
idx1=find(Y<=low);
idx2=find(Y>=high);
size(idx1)
size(idx2)
id=[idx1;idx2];

d(id,:)=[];
t.Data=d;
However, this keeps all values not just within the greena rrows.
I have included my dataHello, I have some data where I my relevant column (Y position) is as this yellow graph

I want to pull out just the data that is between the green arrows – I do know the actual starting and finishing Y positions (101 and 121 respectively) , but there are other values at these values – but they are constant

My data is actually in an array, and the Yposition is column 2, but I delete the 1st column, so the column becomes 1
this was my attempt to remove all rows of my array outside of these green arrows.

% Get user input to the real starting and finishing Y values
prompt = {‘Col:’,’start:’,’finish:’};
dlgtitle = ‘Keep Rows Between’;
dims = [1 35];
definput = {‘1′,’101′,’121’}; % definput = {‘Region:’,num2str(X),num2str(Y)};
answer = inputdlg(prompt,dlgtitle,dims,definput);
col=str2num(answer{1}); low=str2num(answer{2}); high=str2num(answer{3});

% Get Current table data
t=app.UITable;
d=t.Data;

d=rmmissing(d); % One wat Remove Nans

Y=table2array(d);
Y=Y(:,col);

low
head(Y)
class(Y)

% Get index of
idx1=find(Y<=low);
idx2=find(Y>=high);
size(idx1)
size(idx2)
id=[idx1;idx2];

d(id,:)=[];
t.Data=d;
However, this keeps all values not just within the greena rrows.
I have included my data Hello, I have some data where I my relevant column (Y position) is as this yellow graph

I want to pull out just the data that is between the green arrows – I do know the actual starting and finishing Y positions (101 and 121 respectively) , but there are other values at these values – but they are constant

My data is actually in an array, and the Yposition is column 2, but I delete the 1st column, so the column becomes 1
this was my attempt to remove all rows of my array outside of these green arrows.

% Get user input to the real starting and finishing Y values
prompt = {‘Col:’,’start:’,’finish:’};
dlgtitle = ‘Keep Rows Between’;
dims = [1 35];
definput = {‘1′,’101′,’121’}; % definput = {‘Region:’,num2str(X),num2str(Y)};
answer = inputdlg(prompt,dlgtitle,dims,definput);
col=str2num(answer{1}); low=str2num(answer{2}); high=str2num(answer{3});

% Get Current table data
t=app.UITable;
d=t.Data;

d=rmmissing(d); % One wat Remove Nans

Y=table2array(d);
Y=Y(:,col);

low
head(Y)
class(Y)

% Get index of
idx1=find(Y<=low);
idx2=find(Y>=high);
size(idx1)
size(idx2)
id=[idx1;idx2];

d(id,:)=[];
t.Data=d;
However, this keeps all values not just within the greena rrows.
I have included my data find, table2array MATLAB Answers — New Questions

​

As the attached picture shows , when I double click a variable on workspace I cannot see the exact number. Thanks for any one who can helps.As the attached picture shows , when I double click a variable on workspace I cannot see the exact number. Thanks for any one who can helps. As the attached picture shows , when I double click a variable on workspace I cannot see the exact number. Thanks for any one who can helps. workspace, variable, variable display MATLAB Answers — New Questions

​

I’m working on modeling an advanced energy management system for electric vehicles using MATLAB/Simulink. My project involves integrating batteries and supercapacitors to manage energy efficiently, where supercapacitors handle high energy demands during acceleration and store energy recovered from braking. Additionally, I’m incorporating a backup power source using fuel and exploring the integration of renewable energy sources like solar power. I’m facing challenges in connecting and configuring these components correctly in Simulink. Can someone guide me on the best practices for setting up this hybrid system and suggest specific blocks or configurations that would work efficiently for this application?I’m working on modeling an advanced energy management system for electric vehicles using MATLAB/Simulink. My project involves integrating batteries and supercapacitors to manage energy efficiently, where supercapacitors handle high energy demands during acceleration and store energy recovered from braking. Additionally, I’m incorporating a backup power source using fuel and exploring the integration of renewable energy sources like solar power. I’m facing challenges in connecting and configuring these components correctly in Simulink. Can someone guide me on the best practices for setting up this hybrid system and suggest specific blocks or configurations that would work efficiently for this application? I’m working on modeling an advanced energy management system for electric vehicles using MATLAB/Simulink. My project involves integrating batteries and supercapacitors to manage energy efficiently, where supercapacitors handle high energy demands during acceleration and store energy recovered from braking. Additionally, I’m incorporating a backup power source using fuel and exploring the integration of renewable energy sources like solar power. I’m facing challenges in connecting and configuring these components correctly in Simulink. Can someone guide me on the best practices for setting up this hybrid system and suggest specific blocks or configurations that would work efficiently for this application? all MATLAB Answers — New Questions

​

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 advancedHello 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 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 matlab, fft, waterfalldiagramm, campbelldiagramm, electric_motor_control, pmsm MATLAB Answers — New Questions

​

Why do I receive License Manager Error -9?Why do I receive License Manager Error -9? Why do I receive License Manager Error -9? error -9 MATLAB Answers — New Questions

​

Dear staff,
I’m learning how to implement a TGI model in Simbiology using a scripting approach. For that purpose I started to implement my model (based on the example tumor_growth_fitPKPD_completed.sbproj) just with a modified growth component and I’m trying to fit a dataset. Attached you can find the dataset (synth_data.mat) and the TGI script (test_code.m). As you can see the initial tumor mass of different patients varies quite a lot so I can’t use a unique initialization for the variable w0. I would like to understand if (and how) it is possible to pass to the model the initial tumor mass for each ID, in order to avoid to fit this parameter.
Thank you very much for your help.Dear staff,
I’m learning how to implement a TGI model in Simbiology using a scripting approach. For that purpose I started to implement my model (based on the example tumor_growth_fitPKPD_completed.sbproj) just with a modified growth component and I’m trying to fit a dataset. Attached you can find the dataset (synth_data.mat) and the TGI script (test_code.m). As you can see the initial tumor mass of different patients varies quite a lot so I can’t use a unique initialization for the variable w0. I would like to understand if (and how) it is possible to pass to the model the initial tumor mass for each ID, in order to avoid to fit this parameter.
Thank you very much for your help. Dear staff,
I’m learning how to implement a TGI model in Simbiology using a scripting approach. For that purpose I started to implement my model (based on the example tumor_growth_fitPKPD_completed.sbproj) just with a modified growth component and I’m trying to fit a dataset. Attached you can find the dataset (synth_data.mat) and the TGI script (test_code.m). As you can see the initial tumor mass of different patients varies quite a lot so I can’t use a unique initialization for the variable w0. I would like to understand if (and how) it is possible to pass to the model the initial tumor mass for each ID, in order to avoid to fit this parameter.
Thank you very much for your help. simbiology MATLAB Answers — New Questions

​