Month: October 2024
How to resolve “license file selected seems to be invalid” error
Post Content Post Content license MATLAB Answers — New Questions
How can I communicate between a Pixhawk controller and an actuator module using UAVCAN?
I am now building a UAV whose controller is Pixhawk 4X or 6X.
UAVCAN communication is available to drive the main motors.
According to the following page, UAVCAN is used to read GPS data.
https://jp.mathworks.com/help/uav/px4/ref/read-gps-uavcan-example.html
However, I cannot find the detailed configuration of the UAVCAN.
Would you tell me how to use the UAVCAN communication using some blocks in UAV Toolbox Support Package for PX4 Autopilot such as "uORB Read" block?I am now building a UAV whose controller is Pixhawk 4X or 6X.
UAVCAN communication is available to drive the main motors.
According to the following page, UAVCAN is used to read GPS data.
https://jp.mathworks.com/help/uav/px4/ref/read-gps-uavcan-example.html
However, I cannot find the detailed configuration of the UAVCAN.
Would you tell me how to use the UAVCAN communication using some blocks in UAV Toolbox Support Package for PX4 Autopilot such as "uORB Read" block? I am now building a UAV whose controller is Pixhawk 4X or 6X.
UAVCAN communication is available to drive the main motors.
According to the following page, UAVCAN is used to read GPS data.
https://jp.mathworks.com/help/uav/px4/ref/read-gps-uavcan-example.html
However, I cannot find the detailed configuration of the UAVCAN.
Would you tell me how to use the UAVCAN communication using some blocks in UAV Toolbox Support Package for PX4 Autopilot such as "uORB Read" block? uavcan, uav toolbox support package, pixhawk MATLAB Answers — New Questions
How reference one file to show it in multiple folders
Hello everyone,
Can I reference in a simple way one file to multiple folder? I have two docx file that should appear in multiple folders, I would like to avoid to update that files many times in each folders, how can I do?
Hello everyone,Can I reference in a simple way one file to multiple folder? I have two docx file that should appear in multiple folders, I would like to avoid to update that files many times in each folders, how can I do? Read More
How can i group the list command using SPFx ListView Command Set extension via manifest file
How can i group the list command using SPFx ListView Command Set extension via manifest file for SharePoint Add-in? Attached the image how i want. Below is my manifest file
{
“$schema”: “https://developer.microsoft.com/json-schemas/spfx/command-set-extension-manifest.schema.json“,
“id”: “e0947d89-adb7-48cc-bb9b-40a6a676b6d3”,
“alias”: “DocumentLibraryCommandSet”,
“componentType”: “Extension”,
“extensionType”: “ListViewCommandSet”,
“version”: “*”,
“manifestVersion”: 2,
“requiresCustomScript”: false,
“items”: {
“COMMAND_1”: {
“title”: { “default”: “Command One” },
“iconImageUrl”: “icons/request.png”,
“type”: “command”,
“commandGroup”: “CUSTOM_GROUP”,
“children”: [
{
“title”: “Child Command 1”,
“command”: “COMMAND_CHILD_1”
},
{
“title”: “Child Command 2”,
“command”: “COMMAND_CHILD_2”
}
]
},
“COMMAND_2”: {
“title”: { “default”: “Command Two” },
“iconImageUrl”: “icons/cancel.png”,
“type”: “command”
}
}
}
How can i group the list command using SPFx ListView Command Set extension via manifest file for SharePoint Add-in? Attached the image how i want. Below is my manifest file{“$schema”: “https://developer.microsoft.com/json-schemas/spfx/command-set-extension-manifest.schema.json”,”id”: “e0947d89-adb7-48cc-bb9b-40a6a676b6d3″,”alias”: “DocumentLibraryCommandSet”,”componentType”: “Extension”,”extensionType”: “ListViewCommandSet”,”version”: “*”,”manifestVersion”: 2,”requiresCustomScript”: false,”items”: {“COMMAND_1”: {“title”: { “default”: “Command One” },”iconImageUrl”: “icons/request.png”,”type”: “command”,”commandGroup”: “CUSTOM_GROUP”,”children”: [{“title”: “Child Command 1″,”command”: “COMMAND_CHILD_1”},{“title”: “Child Command 2″,”command”: “COMMAND_CHILD_2″}]},”COMMAND_2”: {“title”: { “default”: “Command Two” },”iconImageUrl”: “icons/cancel.png”,”type”: “command”}}} Read More
Trajectory of a 3d body
Problem Statement: –
I have a 3d body which i have launched from cruising aircraft travelling at 200m/sec at 10km altitude. The launch speed is 50m/sec and it is launched at 45degrees angle pointing downwards. I want to create the whole trajectory of the body using everything like change in pressure, density, drag, lift etc. I created a code but it is not giving the accurate answer.
the code: –
function aircraft_launch_trajectory_6DOF_advanced()
% Constants
g = 9.81; % Acceleration due to gravity (m/s^2)
% Convert given units to SI
m = 49.964 * 0.0283495; % Mass (kg) (1 ouncemass = 0.0283495 kg)
A = 68.166 * 0.00064516; % Area (m^2) (1 in^2 = 0.00064516 m^2)
Izz = 19.58 * 0.0283495 * 0.0254^2; % Moment of inertia (kg*m^2)
% Body dimensions
length = 1.80 * 0.0254; % Length (m)
width = 3.30 * 0.0254; % Width (m)
height = 8.00 * 0.0254; % Height (m)
% Center of mass location (assuming from the base of the body)
cm_z = 1.708 * 0.0254; % (m)
% Assumed center of pressure location (you may need to adjust this)
cp_z = 4.5 * 0.0254; % (m)
% Distance between center of mass and center of pressure
d_cp_cm = cp_z – cm_z;
% Initial conditions
h0 = 10000; % Initial height (m)
v0_aircraft = 200; % Aircraft velocity (m/s)
v0_launch = 50; % Launch velocity (m/s)
angle = -45; % Launch angle (degrees)
omega0 = 0; % Initial angular velocity (rad/s)
% Calculate initial velocity components
v0x = v0_aircraft – v0_launch * cosd(angle);
v0y = -v0_launch * sind(angle);
% Initial state vector [x, y, vx, vy, theta, omega]
y0 = [0; h0; v0x; v0y; deg2rad(angle); omega0];
% Time settings
tspan = [0 300]; % Time span for integration
% ODE solver options
options = odeset(‘Events’, @ground_impact, ‘RelTol’, 1e-6, ‘AbsTol’, 1e-6);
% Solve ODE using ode45 (4th order Runge-Kutta)
[t, y] = ode45(@(t, y) equations_of_motion_6DOF_advanced(t, y, m, Izz, A, d_cp_cm, length), tspan, y0, options);
% Extract results
x = y(:, 1);
h = y(:, 2);
vx = y(:, 3);
vy = y(:, 4);
theta = y(:, 5);
omega = y(:, 6);
% Calculate velocity magnitude
v = sqrt(vx.^2 + vy.^2);
% Plot trajectory
figure;
plot(x/1000, h/1000);
xlabel(‘Horizontal distance (km)’);
ylabel(‘Altitude (km)’);
title(‘2D Trajectory of Body (Advanced 6DOF Simulation)’);
grid on;
% Plot velocity
figure;
plot(t, v);
xlabel(‘Time (s)’);
ylabel(‘Velocity (m/s)’);
title(‘Velocity vs Time’);
grid on;
% Plot angular position and velocity
figure;
subplot(2,1,1);
plot(t, rad2deg(theta));
ylabel(‘Angular Position (degrees)’);
title(‘Angular Motion’);
grid on;
subplot(2,1,2);
plot(t, rad2deg(omega));
xlabel(‘Time (s)’);
ylabel(‘Angular Velocity (deg/s)’);
grid on;
% Display results
disp([‘Time of flight: ‘, num2str(t(end)), ‘ s’]);
disp([‘Horizontal distance traveled: ‘, num2str(x(end)/1000), ‘ km’]);
disp([‘Impact velocity: ‘, num2str(v(end)), ‘ m/s’]);
disp([‘Final angle: ‘, num2str(rad2deg(theta(end))), ‘ degrees’]);
disp([‘Max altitude: ‘, num2str(max(h)/1000), ‘ km’]);
disp([‘Max velocity: ‘, num2str(max(v)), ‘ m/s’]);
% Output state at specific times for debugging
debug_times = [0, 5, 10, 30, 60, 120];
for debug_time = debug_times
idx = find(t >= debug_time, 1);
if ~isempty(idx)
disp([‘State at t = ‘, num2str(t(idx)), ‘ s:’]);
disp([‘ Position: (‘, num2str(x(idx)/1000), ‘ km, ‘, num2str(h(idx)/1000), ‘ km)’]);
disp([‘ Velocity: ‘, num2str(v(idx)), ‘ m/s’]);
disp([‘ Angle: ‘, num2str(rad2deg(theta(idx))), ‘ degrees’]);
end
end
end
function dydt = equations_of_motion_6DOF_advanced(t, y, m, Izz, A, d_cp_cm, length)
% Unpack state vector
x = y(1);
h = y(2);
vx = y(3);
vy = y(4);
theta = y(5);
omega = y(6);
% Calculate velocity magnitude and angle of attack
v = sqrt(vx^2 + vy^2);
alpha = atan2(vy, vx) – theta;
% Calculate air density and speed of sound using improved atmospheric model
[rho, a] = atmosphere_model(h);
% Calculate Mach number
M = v / a;
% Get aerodynamic coefficients based on angle of attack and Mach number
[Cd, Cl, Cm] = aero_coefficients(alpha, M);
% Calculate aerodynamic forces and moment
q = 0.5 * rho * v^2; % Dynamic pressure
Fd = q * Cd * A; % Drag force
Fl = q * Cl * A; % Lift force
M_aero = q * Cm * A * length; % Aerodynamic moment
% Add moment due to center of pressure offset
M_total = M_aero + Fl * d_cp_cm;
% Apply initial boost for first 7 seconds
boost_duration = 7; % seconds
if t <= boost_duration
boost_force = m * 200 / boost_duration; % Force to achieve 200 m/s in 7 seconds
boost_fx = boost_force * cos(theta);
boost_fy = boost_force * sin(theta);
else
boost_fx = 0;
boost_fy = 0;
end
% Calculate accelerations
ax = (-Fd * cos(theta) + Fl * sin(theta) + boost_fx) / m;
ay = -9.81 + (-Fd * sin(theta) – Fl * cos(theta) + boost_fy) / m;
alpha_dot = M_total / Izz; % Angular acceleration
% Return state derivatives
dydt = [vx; vy; ax; ay; omega; alpha_dot];
end
function [rho, a] = atmosphere_model(h)
% Constants
R = 287.05; % Specific gas constant for air (J/(kg*K))
g0 = 9.80665; % Standard gravity (m/s^2)
gamma = 1.4; % Ratio of specific heats for air
% Sea level conditions
T0 = 288.15; % Temperature at sea level (K)
P0 = 101325; % Pressure at sea level (Pa)
if h < 11000 % Troposphere
T = T0 – 0.0065 * h;
P = P0 * (T / T0)^(g0 / (R * 0.0065));
elseif h < 20000 % Lower Stratosphere
T = 216.65;
P = 22632.1 * exp(-g0 * (h – 11000) / (R * T));
else % Upper Stratosphere
T = 216.65 + 0.001 * (h – 20000);
P = 5474.89 * (216.65 / T)^(34.1632);
end
rho = P / (R * T);
a = sqrt(gamma * R * T); % Speed of sound
end
function [Cd, Cl, Cm] = aero_coefficients(alpha, M)
% This function returns aerodynamic coefficients based on angle of attack and Mach number
% Note: These are simplified models and should be replaced with actual data or more complex models
% Convert alpha to degrees for easier calculations
alpha_deg = rad2deg(alpha);
% Drag coefficient
Cd0 = 0.1 + 0.1 * M; % Base drag increases with Mach number
Cdi = 0.01 * alpha_deg^2; % Induced drag
Cd = Cd0 + Cdi;
% Lift coefficient
Cl_slope = 0.1; % Lift curve slope
Cl = Cl_slope * alpha_deg * (1 – M^2/4)^(-0.5); % Prandtl-Glauert correction for compressibility
% Moment coefficient
Cm0 = -0.05; % Zero-lift moment coefficient
Cm_alpha = -0.01; % Moment curve slope
Cm = Cm0 + Cm_alpha * alpha_deg;
% Limit coefficients to realistic ranges
Cd = max(0, min(Cd, 2));
Cl = max(-2, min(Cl, 2));
Cm = max(-1, min(Cm, 1));
end
function [position, isterminal, direction] = ground_impact(t, y)
position = y(2); % Height
isterminal = 1; % Stop integration when ground is reached
direction = -1; % Negative direction (approaching ground)
endProblem Statement: –
I have a 3d body which i have launched from cruising aircraft travelling at 200m/sec at 10km altitude. The launch speed is 50m/sec and it is launched at 45degrees angle pointing downwards. I want to create the whole trajectory of the body using everything like change in pressure, density, drag, lift etc. I created a code but it is not giving the accurate answer.
the code: –
function aircraft_launch_trajectory_6DOF_advanced()
% Constants
g = 9.81; % Acceleration due to gravity (m/s^2)
% Convert given units to SI
m = 49.964 * 0.0283495; % Mass (kg) (1 ouncemass = 0.0283495 kg)
A = 68.166 * 0.00064516; % Area (m^2) (1 in^2 = 0.00064516 m^2)
Izz = 19.58 * 0.0283495 * 0.0254^2; % Moment of inertia (kg*m^2)
% Body dimensions
length = 1.80 * 0.0254; % Length (m)
width = 3.30 * 0.0254; % Width (m)
height = 8.00 * 0.0254; % Height (m)
% Center of mass location (assuming from the base of the body)
cm_z = 1.708 * 0.0254; % (m)
% Assumed center of pressure location (you may need to adjust this)
cp_z = 4.5 * 0.0254; % (m)
% Distance between center of mass and center of pressure
d_cp_cm = cp_z – cm_z;
% Initial conditions
h0 = 10000; % Initial height (m)
v0_aircraft = 200; % Aircraft velocity (m/s)
v0_launch = 50; % Launch velocity (m/s)
angle = -45; % Launch angle (degrees)
omega0 = 0; % Initial angular velocity (rad/s)
% Calculate initial velocity components
v0x = v0_aircraft – v0_launch * cosd(angle);
v0y = -v0_launch * sind(angle);
% Initial state vector [x, y, vx, vy, theta, omega]
y0 = [0; h0; v0x; v0y; deg2rad(angle); omega0];
% Time settings
tspan = [0 300]; % Time span for integration
% ODE solver options
options = odeset(‘Events’, @ground_impact, ‘RelTol’, 1e-6, ‘AbsTol’, 1e-6);
% Solve ODE using ode45 (4th order Runge-Kutta)
[t, y] = ode45(@(t, y) equations_of_motion_6DOF_advanced(t, y, m, Izz, A, d_cp_cm, length), tspan, y0, options);
% Extract results
x = y(:, 1);
h = y(:, 2);
vx = y(:, 3);
vy = y(:, 4);
theta = y(:, 5);
omega = y(:, 6);
% Calculate velocity magnitude
v = sqrt(vx.^2 + vy.^2);
% Plot trajectory
figure;
plot(x/1000, h/1000);
xlabel(‘Horizontal distance (km)’);
ylabel(‘Altitude (km)’);
title(‘2D Trajectory of Body (Advanced 6DOF Simulation)’);
grid on;
% Plot velocity
figure;
plot(t, v);
xlabel(‘Time (s)’);
ylabel(‘Velocity (m/s)’);
title(‘Velocity vs Time’);
grid on;
% Plot angular position and velocity
figure;
subplot(2,1,1);
plot(t, rad2deg(theta));
ylabel(‘Angular Position (degrees)’);
title(‘Angular Motion’);
grid on;
subplot(2,1,2);
plot(t, rad2deg(omega));
xlabel(‘Time (s)’);
ylabel(‘Angular Velocity (deg/s)’);
grid on;
% Display results
disp([‘Time of flight: ‘, num2str(t(end)), ‘ s’]);
disp([‘Horizontal distance traveled: ‘, num2str(x(end)/1000), ‘ km’]);
disp([‘Impact velocity: ‘, num2str(v(end)), ‘ m/s’]);
disp([‘Final angle: ‘, num2str(rad2deg(theta(end))), ‘ degrees’]);
disp([‘Max altitude: ‘, num2str(max(h)/1000), ‘ km’]);
disp([‘Max velocity: ‘, num2str(max(v)), ‘ m/s’]);
% Output state at specific times for debugging
debug_times = [0, 5, 10, 30, 60, 120];
for debug_time = debug_times
idx = find(t >= debug_time, 1);
if ~isempty(idx)
disp([‘State at t = ‘, num2str(t(idx)), ‘ s:’]);
disp([‘ Position: (‘, num2str(x(idx)/1000), ‘ km, ‘, num2str(h(idx)/1000), ‘ km)’]);
disp([‘ Velocity: ‘, num2str(v(idx)), ‘ m/s’]);
disp([‘ Angle: ‘, num2str(rad2deg(theta(idx))), ‘ degrees’]);
end
end
end
function dydt = equations_of_motion_6DOF_advanced(t, y, m, Izz, A, d_cp_cm, length)
% Unpack state vector
x = y(1);
h = y(2);
vx = y(3);
vy = y(4);
theta = y(5);
omega = y(6);
% Calculate velocity magnitude and angle of attack
v = sqrt(vx^2 + vy^2);
alpha = atan2(vy, vx) – theta;
% Calculate air density and speed of sound using improved atmospheric model
[rho, a] = atmosphere_model(h);
% Calculate Mach number
M = v / a;
% Get aerodynamic coefficients based on angle of attack and Mach number
[Cd, Cl, Cm] = aero_coefficients(alpha, M);
% Calculate aerodynamic forces and moment
q = 0.5 * rho * v^2; % Dynamic pressure
Fd = q * Cd * A; % Drag force
Fl = q * Cl * A; % Lift force
M_aero = q * Cm * A * length; % Aerodynamic moment
% Add moment due to center of pressure offset
M_total = M_aero + Fl * d_cp_cm;
% Apply initial boost for first 7 seconds
boost_duration = 7; % seconds
if t <= boost_duration
boost_force = m * 200 / boost_duration; % Force to achieve 200 m/s in 7 seconds
boost_fx = boost_force * cos(theta);
boost_fy = boost_force * sin(theta);
else
boost_fx = 0;
boost_fy = 0;
end
% Calculate accelerations
ax = (-Fd * cos(theta) + Fl * sin(theta) + boost_fx) / m;
ay = -9.81 + (-Fd * sin(theta) – Fl * cos(theta) + boost_fy) / m;
alpha_dot = M_total / Izz; % Angular acceleration
% Return state derivatives
dydt = [vx; vy; ax; ay; omega; alpha_dot];
end
function [rho, a] = atmosphere_model(h)
% Constants
R = 287.05; % Specific gas constant for air (J/(kg*K))
g0 = 9.80665; % Standard gravity (m/s^2)
gamma = 1.4; % Ratio of specific heats for air
% Sea level conditions
T0 = 288.15; % Temperature at sea level (K)
P0 = 101325; % Pressure at sea level (Pa)
if h < 11000 % Troposphere
T = T0 – 0.0065 * h;
P = P0 * (T / T0)^(g0 / (R * 0.0065));
elseif h < 20000 % Lower Stratosphere
T = 216.65;
P = 22632.1 * exp(-g0 * (h – 11000) / (R * T));
else % Upper Stratosphere
T = 216.65 + 0.001 * (h – 20000);
P = 5474.89 * (216.65 / T)^(34.1632);
end
rho = P / (R * T);
a = sqrt(gamma * R * T); % Speed of sound
end
function [Cd, Cl, Cm] = aero_coefficients(alpha, M)
% This function returns aerodynamic coefficients based on angle of attack and Mach number
% Note: These are simplified models and should be replaced with actual data or more complex models
% Convert alpha to degrees for easier calculations
alpha_deg = rad2deg(alpha);
% Drag coefficient
Cd0 = 0.1 + 0.1 * M; % Base drag increases with Mach number
Cdi = 0.01 * alpha_deg^2; % Induced drag
Cd = Cd0 + Cdi;
% Lift coefficient
Cl_slope = 0.1; % Lift curve slope
Cl = Cl_slope * alpha_deg * (1 – M^2/4)^(-0.5); % Prandtl-Glauert correction for compressibility
% Moment coefficient
Cm0 = -0.05; % Zero-lift moment coefficient
Cm_alpha = -0.01; % Moment curve slope
Cm = Cm0 + Cm_alpha * alpha_deg;
% Limit coefficients to realistic ranges
Cd = max(0, min(Cd, 2));
Cl = max(-2, min(Cl, 2));
Cm = max(-1, min(Cm, 1));
end
function [position, isterminal, direction] = ground_impact(t, y)
position = y(2); % Height
isterminal = 1; % Stop integration when ground is reached
direction = -1; % Negative direction (approaching ground)
end Problem Statement: –
I have a 3d body which i have launched from cruising aircraft travelling at 200m/sec at 10km altitude. The launch speed is 50m/sec and it is launched at 45degrees angle pointing downwards. I want to create the whole trajectory of the body using everything like change in pressure, density, drag, lift etc. I created a code but it is not giving the accurate answer.
the code: –
function aircraft_launch_trajectory_6DOF_advanced()
% Constants
g = 9.81; % Acceleration due to gravity (m/s^2)
% Convert given units to SI
m = 49.964 * 0.0283495; % Mass (kg) (1 ouncemass = 0.0283495 kg)
A = 68.166 * 0.00064516; % Area (m^2) (1 in^2 = 0.00064516 m^2)
Izz = 19.58 * 0.0283495 * 0.0254^2; % Moment of inertia (kg*m^2)
% Body dimensions
length = 1.80 * 0.0254; % Length (m)
width = 3.30 * 0.0254; % Width (m)
height = 8.00 * 0.0254; % Height (m)
% Center of mass location (assuming from the base of the body)
cm_z = 1.708 * 0.0254; % (m)
% Assumed center of pressure location (you may need to adjust this)
cp_z = 4.5 * 0.0254; % (m)
% Distance between center of mass and center of pressure
d_cp_cm = cp_z – cm_z;
% Initial conditions
h0 = 10000; % Initial height (m)
v0_aircraft = 200; % Aircraft velocity (m/s)
v0_launch = 50; % Launch velocity (m/s)
angle = -45; % Launch angle (degrees)
omega0 = 0; % Initial angular velocity (rad/s)
% Calculate initial velocity components
v0x = v0_aircraft – v0_launch * cosd(angle);
v0y = -v0_launch * sind(angle);
% Initial state vector [x, y, vx, vy, theta, omega]
y0 = [0; h0; v0x; v0y; deg2rad(angle); omega0];
% Time settings
tspan = [0 300]; % Time span for integration
% ODE solver options
options = odeset(‘Events’, @ground_impact, ‘RelTol’, 1e-6, ‘AbsTol’, 1e-6);
% Solve ODE using ode45 (4th order Runge-Kutta)
[t, y] = ode45(@(t, y) equations_of_motion_6DOF_advanced(t, y, m, Izz, A, d_cp_cm, length), tspan, y0, options);
% Extract results
x = y(:, 1);
h = y(:, 2);
vx = y(:, 3);
vy = y(:, 4);
theta = y(:, 5);
omega = y(:, 6);
% Calculate velocity magnitude
v = sqrt(vx.^2 + vy.^2);
% Plot trajectory
figure;
plot(x/1000, h/1000);
xlabel(‘Horizontal distance (km)’);
ylabel(‘Altitude (km)’);
title(‘2D Trajectory of Body (Advanced 6DOF Simulation)’);
grid on;
% Plot velocity
figure;
plot(t, v);
xlabel(‘Time (s)’);
ylabel(‘Velocity (m/s)’);
title(‘Velocity vs Time’);
grid on;
% Plot angular position and velocity
figure;
subplot(2,1,1);
plot(t, rad2deg(theta));
ylabel(‘Angular Position (degrees)’);
title(‘Angular Motion’);
grid on;
subplot(2,1,2);
plot(t, rad2deg(omega));
xlabel(‘Time (s)’);
ylabel(‘Angular Velocity (deg/s)’);
grid on;
% Display results
disp([‘Time of flight: ‘, num2str(t(end)), ‘ s’]);
disp([‘Horizontal distance traveled: ‘, num2str(x(end)/1000), ‘ km’]);
disp([‘Impact velocity: ‘, num2str(v(end)), ‘ m/s’]);
disp([‘Final angle: ‘, num2str(rad2deg(theta(end))), ‘ degrees’]);
disp([‘Max altitude: ‘, num2str(max(h)/1000), ‘ km’]);
disp([‘Max velocity: ‘, num2str(max(v)), ‘ m/s’]);
% Output state at specific times for debugging
debug_times = [0, 5, 10, 30, 60, 120];
for debug_time = debug_times
idx = find(t >= debug_time, 1);
if ~isempty(idx)
disp([‘State at t = ‘, num2str(t(idx)), ‘ s:’]);
disp([‘ Position: (‘, num2str(x(idx)/1000), ‘ km, ‘, num2str(h(idx)/1000), ‘ km)’]);
disp([‘ Velocity: ‘, num2str(v(idx)), ‘ m/s’]);
disp([‘ Angle: ‘, num2str(rad2deg(theta(idx))), ‘ degrees’]);
end
end
end
function dydt = equations_of_motion_6DOF_advanced(t, y, m, Izz, A, d_cp_cm, length)
% Unpack state vector
x = y(1);
h = y(2);
vx = y(3);
vy = y(4);
theta = y(5);
omega = y(6);
% Calculate velocity magnitude and angle of attack
v = sqrt(vx^2 + vy^2);
alpha = atan2(vy, vx) – theta;
% Calculate air density and speed of sound using improved atmospheric model
[rho, a] = atmosphere_model(h);
% Calculate Mach number
M = v / a;
% Get aerodynamic coefficients based on angle of attack and Mach number
[Cd, Cl, Cm] = aero_coefficients(alpha, M);
% Calculate aerodynamic forces and moment
q = 0.5 * rho * v^2; % Dynamic pressure
Fd = q * Cd * A; % Drag force
Fl = q * Cl * A; % Lift force
M_aero = q * Cm * A * length; % Aerodynamic moment
% Add moment due to center of pressure offset
M_total = M_aero + Fl * d_cp_cm;
% Apply initial boost for first 7 seconds
boost_duration = 7; % seconds
if t <= boost_duration
boost_force = m * 200 / boost_duration; % Force to achieve 200 m/s in 7 seconds
boost_fx = boost_force * cos(theta);
boost_fy = boost_force * sin(theta);
else
boost_fx = 0;
boost_fy = 0;
end
% Calculate accelerations
ax = (-Fd * cos(theta) + Fl * sin(theta) + boost_fx) / m;
ay = -9.81 + (-Fd * sin(theta) – Fl * cos(theta) + boost_fy) / m;
alpha_dot = M_total / Izz; % Angular acceleration
% Return state derivatives
dydt = [vx; vy; ax; ay; omega; alpha_dot];
end
function [rho, a] = atmosphere_model(h)
% Constants
R = 287.05; % Specific gas constant for air (J/(kg*K))
g0 = 9.80665; % Standard gravity (m/s^2)
gamma = 1.4; % Ratio of specific heats for air
% Sea level conditions
T0 = 288.15; % Temperature at sea level (K)
P0 = 101325; % Pressure at sea level (Pa)
if h < 11000 % Troposphere
T = T0 – 0.0065 * h;
P = P0 * (T / T0)^(g0 / (R * 0.0065));
elseif h < 20000 % Lower Stratosphere
T = 216.65;
P = 22632.1 * exp(-g0 * (h – 11000) / (R * T));
else % Upper Stratosphere
T = 216.65 + 0.001 * (h – 20000);
P = 5474.89 * (216.65 / T)^(34.1632);
end
rho = P / (R * T);
a = sqrt(gamma * R * T); % Speed of sound
end
function [Cd, Cl, Cm] = aero_coefficients(alpha, M)
% This function returns aerodynamic coefficients based on angle of attack and Mach number
% Note: These are simplified models and should be replaced with actual data or more complex models
% Convert alpha to degrees for easier calculations
alpha_deg = rad2deg(alpha);
% Drag coefficient
Cd0 = 0.1 + 0.1 * M; % Base drag increases with Mach number
Cdi = 0.01 * alpha_deg^2; % Induced drag
Cd = Cd0 + Cdi;
% Lift coefficient
Cl_slope = 0.1; % Lift curve slope
Cl = Cl_slope * alpha_deg * (1 – M^2/4)^(-0.5); % Prandtl-Glauert correction for compressibility
% Moment coefficient
Cm0 = -0.05; % Zero-lift moment coefficient
Cm_alpha = -0.01; % Moment curve slope
Cm = Cm0 + Cm_alpha * alpha_deg;
% Limit coefficients to realistic ranges
Cd = max(0, min(Cd, 2));
Cl = max(-2, min(Cl, 2));
Cm = max(-1, min(Cm, 1));
end
function [position, isterminal, direction] = ground_impact(t, y)
position = y(2); % Height
isterminal = 1; % Stop integration when ground is reached
direction = -1; % Negative direction (approaching ground)
end trajectory, aerospace MATLAB Answers — New Questions
Initial position detection of sensor less BLDC motor drive
my research is sensor less speed control of BLDC motor drive, it’s possible to detect ZCD of line voltage. but initial starting speed zero time cont detect ZCD points, if start open loop v/f mean rotor will rotate wrong direction at initial.
any body help me to detect the initial rotor position detection and speed control starting correct direction of sensor less speed control of BLDC motor drive?
if you able pls contact me..
thank you,
Reg,
R.Manikandan,M.E.,(Ph.D).,
AP / Department of EEE,
The Kavery Engineering College,
Mecheri (po), Salem(dt), pin:636453
Tamil Nadu, India.
Cell: 9944270473
Email: electricmani@yahoo.co.inmy research is sensor less speed control of BLDC motor drive, it’s possible to detect ZCD of line voltage. but initial starting speed zero time cont detect ZCD points, if start open loop v/f mean rotor will rotate wrong direction at initial.
any body help me to detect the initial rotor position detection and speed control starting correct direction of sensor less speed control of BLDC motor drive?
if you able pls contact me..
thank you,
Reg,
R.Manikandan,M.E.,(Ph.D).,
AP / Department of EEE,
The Kavery Engineering College,
Mecheri (po), Salem(dt), pin:636453
Tamil Nadu, India.
Cell: 9944270473
Email: electricmani@yahoo.co.in my research is sensor less speed control of BLDC motor drive, it’s possible to detect ZCD of line voltage. but initial starting speed zero time cont detect ZCD points, if start open loop v/f mean rotor will rotate wrong direction at initial.
any body help me to detect the initial rotor position detection and speed control starting correct direction of sensor less speed control of BLDC motor drive?
if you able pls contact me..
thank you,
Reg,
R.Manikandan,M.E.,(Ph.D).,
AP / Department of EEE,
The Kavery Engineering College,
Mecheri (po), Salem(dt), pin:636453
Tamil Nadu, India.
Cell: 9944270473
Email: electricmani@yahoo.co.in bldc, simpowersystems, sensor less drive, power_electronics_control, electric_motor_control MATLAB Answers — New Questions
How to do feature extraction using wavelet scattering and then perform neural network classification?
% I want to do feature extraction for signals and then apply ANN classification on them,, can you help please?
Fs = 50;
sf = waveletScattering(‘SignalLength’, 4950, ‘SamplingFrequency’, Fs) ;
s = cell(120,1);
parfor j = 1:120;
r = featureMatrix(sf, f(j, :)) ;
s{j, 1} = r;
end
h = cell2mat(s) ;
b = cell2table(s) ;
c = table2dataset(b);
%the single signal was in one row, but after this it becomes 202rows×10columns, and I can’t apply ANN like this
Any ideas?% I want to do feature extraction for signals and then apply ANN classification on them,, can you help please?
Fs = 50;
sf = waveletScattering(‘SignalLength’, 4950, ‘SamplingFrequency’, Fs) ;
s = cell(120,1);
parfor j = 1:120;
r = featureMatrix(sf, f(j, :)) ;
s{j, 1} = r;
end
h = cell2mat(s) ;
b = cell2table(s) ;
c = table2dataset(b);
%the single signal was in one row, but after this it becomes 202rows×10columns, and I can’t apply ANN like this
Any ideas? % I want to do feature extraction for signals and then apply ANN classification on them,, can you help please?
Fs = 50;
sf = waveletScattering(‘SignalLength’, 4950, ‘SamplingFrequency’, Fs) ;
s = cell(120,1);
parfor j = 1:120;
r = featureMatrix(sf, f(j, :)) ;
s{j, 1} = r;
end
h = cell2mat(s) ;
b = cell2table(s) ;
c = table2dataset(b);
%the single signal was in one row, but after this it becomes 202rows×10columns, and I can’t apply ANN like this
Any ideas? feature extraction, neural network, classification, wavelet, wavelet scattering, ecg signals, ecg recognition, ecg as signature, deep learning, machine learning MATLAB Answers — New Questions
Explore DFT Resolution, need help with answer interpretation
Hello there, I’m relatively new to Matlab and Digital Signal Processing. I’m trying to analyze the frequency resolution of the DFT as a function of the DFT size assuming a rectangular window. Using the a sum of sinusoids sampled at 5 kHz vary the frequency spacing
of the sinusoids to show the frequency resolution for at least 4 different DFT lengths, 256, 512, 1024, and 4096.
I got two versions of code from chatgpt. But I don’t know how to interpret them.
1st version
% MATLAB code to analyze the frequency resolution of the DFT using dftmtx and a rectangular window
fs = 5000; % Sampling frequency (5 kHz)
t_duration = 1; % Signal duration in seconds
t = 0:1/fs:t_duration-1/fs; % Time vector
% DFT sizes to analyze
N_dft = [256, 512, 1024, 4096];
% Frequencies of the sum of sinusoids (vary frequency spacing)
f1 = 1000; % Frequency of the first sinusoid (1 kHz)
f_spacing = [1, 5, 10, 20]; % Frequency spacing between the two sinusoids
f_end = f1 + f_spacing; % Frequency of the second sinusoid
% Prepare figure
figure;
hold on;
for N = N_dft
figure; % Create a new figure for each DFT size
hold on;
for spacing = f_spacing
f2 = f1 + spacing; % Second sinusoid frequency
% Generate the sum of two sinusoids with frequencies f1 and f2
x = sin(2*pi*f1*t) + sin(2*pi*f2*t);
% Apply rectangular window (by taking the first N samples)
x_windowed = x(1:N); % Select the first N samples for DFT calculation
% Generate DFT matrix for size N using dftmtx
DFT_matrix = dftmtx(N);
% Manually compute the DFT using the DFT matrix
X = DFT_matrix * x_windowed(:); % Compute DFT of the windowed signal
% Compute the frequency axis for the current DFT
freq_axis = (0:N-1)*(fs/N);
% Plot the magnitude of the DFT
plot(freq_axis, abs(X), ‘DisplayName’, [‘Spacing = ‘, num2str(spacing), ‘ Hz’]);
end
% Add labels and legend
xlabel(‘Frequency (Hz)’);
ylabel(‘Magnitude’);
title([‘DFT Magnitude Spectrum (N = ‘, num2str(N), ‘)’]);
legend(‘show’);
grid on;
hold off;
end
which gives me these resulted figures.
All the peaks are so similar that I can’t distinguish them at all.
2nd version of code is here:
% Parameters
fs = 5000; % Sampling frequency (Hz)
T = 1; % Signal duration (seconds)
t = 0:1/fs:T-1/fs; % Time vector
% Sinusoids: Frequencies and amplitudes
f1 = 400; % Frequency of first sinusoid (Hz)
f2 = 450; % Frequency of second sinusoid (Hz)
f3 = 500; % Frequency of third sinusoid (Hz)
% Create the signal: sum of three sinusoids
x = sin(2*pi*f1*t) + sin(2*pi*f2*t) + sin(2*pi*f3*t);
% Define DFT lengths
DFT_lengths = [256, 512, 1024, 4096];
% Plot frequency resolution for each DFT length using a rectangular window
figure;
for i = 1:length(DFT_lengths)
N = DFT_lengths(i); % DFT size
rectangular_window = rectwin(N)’;
% Apply rectangular window (implicitly applied since no windowing function)
% Select only the first N samples for each DFT calculation
x_windowed = x(1:N).*rectangular_window;
% Manually compute the DFT
% Compute the N-point DFT
X = (dftmtx(N)*x_windowed(:)).’; % Compute the N-point DFT
X_mag = abs(X); % Magnitude of DFT
% Frequency vector
f = (0:N-1) * (fs / N);
% Plot the DFT magnitude (first half)
subplot(length(DFT_lengths), 1, i);
plot(f(1:N/2), X_mag(1:N/2));
title([‘DFT Magnitude Spectrum with Rectangular Window, N = ‘, num2str(N)]);
xlabel(‘Frequency (Hz)’);
ylabel(‘|X(f)|’);
grid on;
end
% Adjust plot layout
sgtitle(‘Frequency Resolution for Different DFT Sizes Using Rectangular Window’);
%Part 5
% MATLAB code to analyze zero-padding in the DFT using a rectangular window
fs = 5000; % Sampling frequency (5 kHz)
t_duration = 1; % Signal duration in seconds
t = 0:1/fs:t_duration-1/fs; % Time vector
% Window sizes to analyze
window_sizes = [256, 512];
% Zero-padded DFT sizes to analyze
N_dft = [1024, 2048, 4096];
% Frequencies of the sum of sinusoids (vary frequency spacing)
f1 = 1000; % Frequency of the first sinusoid (1 kHz)
f_spacing = [5, 10]; % Frequency spacing between the two sinusoids
f_end = f1 + f_spacing; % Frequency of the second sinusoid
% Prepare figure
for window_size = window_sizes
figure; % Create a new figure for each window size
hold on;
for N = N_dft
for spacing = f_spacing
f2 = f1 + spacing; % Second sinusoid frequency
% Generate the sum of two sinusoids with frequencies f1 and f2
x = sin(2*pi*f1*t) + sin(2*pi*f2*t);
% Apply rectangular window (by taking the first window_size samples)
x_windowed = x(1:window_size); % Select the first window_size samples
% Zero-pad the signal if DFT size is larger than window size
x_padded = [x_windowed, zeros(1, N – window_size)];
% Generate DFT matrix for size N using dftmtx
DFT_matrix = dftmtx(N);
% Manually compute the DFT using the DFT matrix
X = DFT_matrix * x_padded(:); % Compute DFT of the windowed and zero-padded signal
% Compute the frequency axis for the current DFT
freq_axis = (0:N-1)*(fs/N);
% Plot the magnitude of the DFT
plot(freq_axis, abs(X), ‘DisplayName’, [‘Spacing = ‘, num2str(spacing), ‘ Hz, N = ‘, num2str(N)]);
end
end
% Add labels and legend
xlabel(‘Frequency (Hz)’);
ylabel(‘Magnitude’);
title([‘Zero-Padded DFT Magnitude Spectrum (Window Size = ‘, num2str(window_size), ‘)’]);
legend(‘show’);
grid on;
hold off;
end
with the second version of code I can see the differences among different sample sizes, as more sampling yields more sharp signals, but why is that and what differences it will make?Hello there, I’m relatively new to Matlab and Digital Signal Processing. I’m trying to analyze the frequency resolution of the DFT as a function of the DFT size assuming a rectangular window. Using the a sum of sinusoids sampled at 5 kHz vary the frequency spacing
of the sinusoids to show the frequency resolution for at least 4 different DFT lengths, 256, 512, 1024, and 4096.
I got two versions of code from chatgpt. But I don’t know how to interpret them.
1st version
% MATLAB code to analyze the frequency resolution of the DFT using dftmtx and a rectangular window
fs = 5000; % Sampling frequency (5 kHz)
t_duration = 1; % Signal duration in seconds
t = 0:1/fs:t_duration-1/fs; % Time vector
% DFT sizes to analyze
N_dft = [256, 512, 1024, 4096];
% Frequencies of the sum of sinusoids (vary frequency spacing)
f1 = 1000; % Frequency of the first sinusoid (1 kHz)
f_spacing = [1, 5, 10, 20]; % Frequency spacing between the two sinusoids
f_end = f1 + f_spacing; % Frequency of the second sinusoid
% Prepare figure
figure;
hold on;
for N = N_dft
figure; % Create a new figure for each DFT size
hold on;
for spacing = f_spacing
f2 = f1 + spacing; % Second sinusoid frequency
% Generate the sum of two sinusoids with frequencies f1 and f2
x = sin(2*pi*f1*t) + sin(2*pi*f2*t);
% Apply rectangular window (by taking the first N samples)
x_windowed = x(1:N); % Select the first N samples for DFT calculation
% Generate DFT matrix for size N using dftmtx
DFT_matrix = dftmtx(N);
% Manually compute the DFT using the DFT matrix
X = DFT_matrix * x_windowed(:); % Compute DFT of the windowed signal
% Compute the frequency axis for the current DFT
freq_axis = (0:N-1)*(fs/N);
% Plot the magnitude of the DFT
plot(freq_axis, abs(X), ‘DisplayName’, [‘Spacing = ‘, num2str(spacing), ‘ Hz’]);
end
% Add labels and legend
xlabel(‘Frequency (Hz)’);
ylabel(‘Magnitude’);
title([‘DFT Magnitude Spectrum (N = ‘, num2str(N), ‘)’]);
legend(‘show’);
grid on;
hold off;
end
which gives me these resulted figures.
All the peaks are so similar that I can’t distinguish them at all.
2nd version of code is here:
% Parameters
fs = 5000; % Sampling frequency (Hz)
T = 1; % Signal duration (seconds)
t = 0:1/fs:T-1/fs; % Time vector
% Sinusoids: Frequencies and amplitudes
f1 = 400; % Frequency of first sinusoid (Hz)
f2 = 450; % Frequency of second sinusoid (Hz)
f3 = 500; % Frequency of third sinusoid (Hz)
% Create the signal: sum of three sinusoids
x = sin(2*pi*f1*t) + sin(2*pi*f2*t) + sin(2*pi*f3*t);
% Define DFT lengths
DFT_lengths = [256, 512, 1024, 4096];
% Plot frequency resolution for each DFT length using a rectangular window
figure;
for i = 1:length(DFT_lengths)
N = DFT_lengths(i); % DFT size
rectangular_window = rectwin(N)’;
% Apply rectangular window (implicitly applied since no windowing function)
% Select only the first N samples for each DFT calculation
x_windowed = x(1:N).*rectangular_window;
% Manually compute the DFT
% Compute the N-point DFT
X = (dftmtx(N)*x_windowed(:)).’; % Compute the N-point DFT
X_mag = abs(X); % Magnitude of DFT
% Frequency vector
f = (0:N-1) * (fs / N);
% Plot the DFT magnitude (first half)
subplot(length(DFT_lengths), 1, i);
plot(f(1:N/2), X_mag(1:N/2));
title([‘DFT Magnitude Spectrum with Rectangular Window, N = ‘, num2str(N)]);
xlabel(‘Frequency (Hz)’);
ylabel(‘|X(f)|’);
grid on;
end
% Adjust plot layout
sgtitle(‘Frequency Resolution for Different DFT Sizes Using Rectangular Window’);
%Part 5
% MATLAB code to analyze zero-padding in the DFT using a rectangular window
fs = 5000; % Sampling frequency (5 kHz)
t_duration = 1; % Signal duration in seconds
t = 0:1/fs:t_duration-1/fs; % Time vector
% Window sizes to analyze
window_sizes = [256, 512];
% Zero-padded DFT sizes to analyze
N_dft = [1024, 2048, 4096];
% Frequencies of the sum of sinusoids (vary frequency spacing)
f1 = 1000; % Frequency of the first sinusoid (1 kHz)
f_spacing = [5, 10]; % Frequency spacing between the two sinusoids
f_end = f1 + f_spacing; % Frequency of the second sinusoid
% Prepare figure
for window_size = window_sizes
figure; % Create a new figure for each window size
hold on;
for N = N_dft
for spacing = f_spacing
f2 = f1 + spacing; % Second sinusoid frequency
% Generate the sum of two sinusoids with frequencies f1 and f2
x = sin(2*pi*f1*t) + sin(2*pi*f2*t);
% Apply rectangular window (by taking the first window_size samples)
x_windowed = x(1:window_size); % Select the first window_size samples
% Zero-pad the signal if DFT size is larger than window size
x_padded = [x_windowed, zeros(1, N – window_size)];
% Generate DFT matrix for size N using dftmtx
DFT_matrix = dftmtx(N);
% Manually compute the DFT using the DFT matrix
X = DFT_matrix * x_padded(:); % Compute DFT of the windowed and zero-padded signal
% Compute the frequency axis for the current DFT
freq_axis = (0:N-1)*(fs/N);
% Plot the magnitude of the DFT
plot(freq_axis, abs(X), ‘DisplayName’, [‘Spacing = ‘, num2str(spacing), ‘ Hz, N = ‘, num2str(N)]);
end
end
% Add labels and legend
xlabel(‘Frequency (Hz)’);
ylabel(‘Magnitude’);
title([‘Zero-Padded DFT Magnitude Spectrum (Window Size = ‘, num2str(window_size), ‘)’]);
legend(‘show’);
grid on;
hold off;
end
with the second version of code I can see the differences among different sample sizes, as more sampling yields more sharp signals, but why is that and what differences it will make? Hello there, I’m relatively new to Matlab and Digital Signal Processing. I’m trying to analyze the frequency resolution of the DFT as a function of the DFT size assuming a rectangular window. Using the a sum of sinusoids sampled at 5 kHz vary the frequency spacing
of the sinusoids to show the frequency resolution for at least 4 different DFT lengths, 256, 512, 1024, and 4096.
I got two versions of code from chatgpt. But I don’t know how to interpret them.
1st version
% MATLAB code to analyze the frequency resolution of the DFT using dftmtx and a rectangular window
fs = 5000; % Sampling frequency (5 kHz)
t_duration = 1; % Signal duration in seconds
t = 0:1/fs:t_duration-1/fs; % Time vector
% DFT sizes to analyze
N_dft = [256, 512, 1024, 4096];
% Frequencies of the sum of sinusoids (vary frequency spacing)
f1 = 1000; % Frequency of the first sinusoid (1 kHz)
f_spacing = [1, 5, 10, 20]; % Frequency spacing between the two sinusoids
f_end = f1 + f_spacing; % Frequency of the second sinusoid
% Prepare figure
figure;
hold on;
for N = N_dft
figure; % Create a new figure for each DFT size
hold on;
for spacing = f_spacing
f2 = f1 + spacing; % Second sinusoid frequency
% Generate the sum of two sinusoids with frequencies f1 and f2
x = sin(2*pi*f1*t) + sin(2*pi*f2*t);
% Apply rectangular window (by taking the first N samples)
x_windowed = x(1:N); % Select the first N samples for DFT calculation
% Generate DFT matrix for size N using dftmtx
DFT_matrix = dftmtx(N);
% Manually compute the DFT using the DFT matrix
X = DFT_matrix * x_windowed(:); % Compute DFT of the windowed signal
% Compute the frequency axis for the current DFT
freq_axis = (0:N-1)*(fs/N);
% Plot the magnitude of the DFT
plot(freq_axis, abs(X), ‘DisplayName’, [‘Spacing = ‘, num2str(spacing), ‘ Hz’]);
end
% Add labels and legend
xlabel(‘Frequency (Hz)’);
ylabel(‘Magnitude’);
title([‘DFT Magnitude Spectrum (N = ‘, num2str(N), ‘)’]);
legend(‘show’);
grid on;
hold off;
end
which gives me these resulted figures.
All the peaks are so similar that I can’t distinguish them at all.
2nd version of code is here:
% Parameters
fs = 5000; % Sampling frequency (Hz)
T = 1; % Signal duration (seconds)
t = 0:1/fs:T-1/fs; % Time vector
% Sinusoids: Frequencies and amplitudes
f1 = 400; % Frequency of first sinusoid (Hz)
f2 = 450; % Frequency of second sinusoid (Hz)
f3 = 500; % Frequency of third sinusoid (Hz)
% Create the signal: sum of three sinusoids
x = sin(2*pi*f1*t) + sin(2*pi*f2*t) + sin(2*pi*f3*t);
% Define DFT lengths
DFT_lengths = [256, 512, 1024, 4096];
% Plot frequency resolution for each DFT length using a rectangular window
figure;
for i = 1:length(DFT_lengths)
N = DFT_lengths(i); % DFT size
rectangular_window = rectwin(N)’;
% Apply rectangular window (implicitly applied since no windowing function)
% Select only the first N samples for each DFT calculation
x_windowed = x(1:N).*rectangular_window;
% Manually compute the DFT
% Compute the N-point DFT
X = (dftmtx(N)*x_windowed(:)).’; % Compute the N-point DFT
X_mag = abs(X); % Magnitude of DFT
% Frequency vector
f = (0:N-1) * (fs / N);
% Plot the DFT magnitude (first half)
subplot(length(DFT_lengths), 1, i);
plot(f(1:N/2), X_mag(1:N/2));
title([‘DFT Magnitude Spectrum with Rectangular Window, N = ‘, num2str(N)]);
xlabel(‘Frequency (Hz)’);
ylabel(‘|X(f)|’);
grid on;
end
% Adjust plot layout
sgtitle(‘Frequency Resolution for Different DFT Sizes Using Rectangular Window’);
%Part 5
% MATLAB code to analyze zero-padding in the DFT using a rectangular window
fs = 5000; % Sampling frequency (5 kHz)
t_duration = 1; % Signal duration in seconds
t = 0:1/fs:t_duration-1/fs; % Time vector
% Window sizes to analyze
window_sizes = [256, 512];
% Zero-padded DFT sizes to analyze
N_dft = [1024, 2048, 4096];
% Frequencies of the sum of sinusoids (vary frequency spacing)
f1 = 1000; % Frequency of the first sinusoid (1 kHz)
f_spacing = [5, 10]; % Frequency spacing between the two sinusoids
f_end = f1 + f_spacing; % Frequency of the second sinusoid
% Prepare figure
for window_size = window_sizes
figure; % Create a new figure for each window size
hold on;
for N = N_dft
for spacing = f_spacing
f2 = f1 + spacing; % Second sinusoid frequency
% Generate the sum of two sinusoids with frequencies f1 and f2
x = sin(2*pi*f1*t) + sin(2*pi*f2*t);
% Apply rectangular window (by taking the first window_size samples)
x_windowed = x(1:window_size); % Select the first window_size samples
% Zero-pad the signal if DFT size is larger than window size
x_padded = [x_windowed, zeros(1, N – window_size)];
% Generate DFT matrix for size N using dftmtx
DFT_matrix = dftmtx(N);
% Manually compute the DFT using the DFT matrix
X = DFT_matrix * x_padded(:); % Compute DFT of the windowed and zero-padded signal
% Compute the frequency axis for the current DFT
freq_axis = (0:N-1)*(fs/N);
% Plot the magnitude of the DFT
plot(freq_axis, abs(X), ‘DisplayName’, [‘Spacing = ‘, num2str(spacing), ‘ Hz, N = ‘, num2str(N)]);
end
end
% Add labels and legend
xlabel(‘Frequency (Hz)’);
ylabel(‘Magnitude’);
title([‘Zero-Padded DFT Magnitude Spectrum (Window Size = ‘, num2str(window_size), ‘)’]);
legend(‘show’);
grid on;
hold off;
end
with the second version of code I can see the differences among different sample sizes, as more sampling yields more sharp signals, but why is that and what differences it will make? dft, frequency resolution, digital signal processing, matlab MATLAB Answers — New Questions
Demystifying Log Ingestion API
In this blog, I’m going to delve into ingesting and transforming application logs to log analytics workspace using Log Ingestion API technique. Before we jump into the details, let’s explore the different types of ingestion techniques based on the application log types.
If applications and services logs information to text files instead of standard logging services such as Windows Event log or Syslog, Custom Text Logs can be leveraged to ingest such text file logs in log analytics workspace.
If applications and services logs information to JSON files instead of standard logging services such as Windows Event log or Syslog, Custom JSON Logs can be leveraged to ingest such text file logs in log analytics workspace.
If an application understands how to send data to an API then you can leverage Log Ingestion API to send the data to Log Analytics Workspace.
Custom Text/JSON Logs are out of scope for this blog, I might write a separate blog dedicated to these techniques later. In this blog, my focus will be on streaming data to log analytics workspace using Log Ingestion API and transforming the data for optimal usage.
Note: This blog aims to demonstrate how to ingest logs using the log ingestion API. To keep things straightforward, I’ll refer to our public documentation.
Overview of Log Ingestion API
The Logs Ingestion API in Azure Monitor lets you send data to a Log Analytics workspace using either a REST API call or client libraries. The API allows you to send data to supported Azure tables or to custom tables that you create.
Data can be sent to the Logs Ingestion API from any application that can make a REST API call. This may be a custom application that you create, or it may be an application or agent that understands how to send data to the API.
Note: The data sent by your application to the API must be formatted in JSON and match the structure expected by the DCR. It doesn’t necessarily need to match the structure of the target table because the DCR can include a transformation to convert the data to match the table’s structure.
Below diagram shows the flow of data from an application to data ingestion pipeline where transformation is applied, the transformed data is then ingested to the log analytics workspace.
Advantages of Log Ingestion API
Since we’ve covered the overview of Log Ingestion API, let’s understand why you should consider Log Ingestion API over traditional Data Collector API.
Since it requires a DCR-based custom table to ingest the data; Log Ingestion API supports transformations, which enable you to modify the data before it’s ingested into the destination table, including filtering and data manipulation.
You can send the data to multiple destinations.
Enables you to manage the destination table schema, including column names, and whether to add new columns to the destination table when the source data schema changes.
HTTP Data Collector API is getting deprecated and it’s replaced by Log Ingestion API so it’s a future proof solution.
Prerequisites for Log Ingestion API
There are certain prerequisites for configuring Log Ingestion API which are as follows:
App registration and secret: The application registration is for the authentication of the API call. It must be granted permissions to the DCR. The API call includes the Application (client) ID and Directory (tenant) ID of the application and the Value of an application secret.
Data collection endpoint (DCE): We will need endpoint URI of the DCE in our configuration. The endpoint URI in a data collection endpoint (DCE) is a connection where data sources send collected data for processing and ingestion into Azure Monitor.
Table in Log Analytics workspace: The table in the Log Analytics workspace must exist before you can send data to it. You can use one of the supported Azure tables or create a custom table using any of the available methods.
Data collection rule (DCR): Azure Monitor uses the Data collection rule (DCR) to understand the structure of the incoming data and what to do with it.
Since the “What” and “Why” part is covered, let’s understand the “How” part.
Let’s look at how to send data to Azure Monitor using Log Ingestion API. We have detailed public documentation on this topic where we have covered approach based on configuration with Azure Portal and ARM template.
I will be leveraging our public documentation in this blog. Let’s get started:
Create App Registration and Secret
Browse to Entra ID > App Registrations > Create a new registration as shown below.
We will need Application (client) ID and the Directory (tenant) ID at later stage.
Now we need to generate an application client secret, which is similar to creating a password to use with a username. Select Certificates & secrets > New client secret. Give the secret a name to identify its purpose and select an Expires duration.
Ensure that you record the Value of secret because you can’t recover it after you move away from this page.
Create Data Collection Endpoint
Browse to Data Collection Endpoint and create a new Data Collection Endpoint as shown below:
Copy the endpoint URI from the DCE, we will need it at later stage.
Create new DCR-based custom table in Log Analytics workspace
To avoid any duplication in the blog, I would recommend going through our public documentation for this step. The steps are very well articulated in our documentation for creating a DCR-based custom table and Data Collection Rule.
In my lab, the custom table name is “CustomLogDemo_CL” which is attached to “DEMO-DCR-UI” Data Collection Rule.
I have also applied some parsing and filtering to get better visibility of the data.
Browse to the Data Collection Rule created for this exercise and copy the immutable ID, we will need it at later stage. An immutable ID is an unique identifier for the data collection rule.
Assign Permission to DCR
The final step is to give the application permission to use the DCR. Any application that uses the correct application ID and application key can now send data to the new DCE and DCR.
Browse to Data Collection Rule > IAM > Add role assignment > Select Monitoring Metrics Publisher > Select the application that we created.
Select Review and assign the permission.
Follow the instructions in our public documentation to generate the sample data.
Once sample data is generated, send the sample data to Azure Monitor by following our public documentation.
Let’s check if it really works
Considering all the configurations are correct, the logs should appear in the log analytics workspace as shown below:
Transforming the data
Now that we are successfully sending data to DCE-based custom table using Log Ingestion API, let’s implement some transformation.
You can refer my blogs on transformation techniques for better understanding:
Workspace & DCR Transformation Simplified – Microsoft Community Hub
Save ingestion costs by splitting logs into multiple tables and opting for the basic tier! – Microsoft Community Hub
I will directly jump into the transformation which I’m going to implement in my custom table.
In my use case, I’m going to split the data into 2 tables, CustomLogDemo_CL and Split_CustomLogDemo_CL. Split_CustomLogDemo_CL table will have only those rows where the ResponseCode is 200, rest of the data will be ingested into CustomLogDemo_CL table.
Detailed steps as follows:
Browse to DCR > Choose the DCR created for this exercise > Automation > Export Template > Deploy > Edit Template.
As per below template, I’ve created 2 streams for this operation.
One stream is sending all the data except where the response code is 200 to CustomLogDemo_CL table.
Another stream is sending the data with response code 200 to Split_CustomLogDemo_CL table.
Let’s validate if it works as expected.
Browse to Log Analytics Workspace > Logs and check both the tables as shown below:
As we can see, we have split the incoming data as desired. You can optimize the table as per your needs with the help of transformation.
I would like to thank my colleague Kaustubh Dwivedi @Kausd for sharing his expertise and co-authoring this blog along with me.
References:
Logs Ingestion API in Azure Monitor – Azure Monitor | Microsoft Learn
Tutorial: Send data to Azure Monitor Logs with Logs ingestion API (Azure portal) – Azure Monitor | Microsoft Learn
Azure Monitor: How To Use Managed Identity with Log Ingestion API – Microsoft Community Hub
Azure Monitor: Logs Ingestion API Tips & Tricks – Microsoft Community Hub
Microsoft Tech Community – Latest Blogs –Read More
PostgreSQL 17 Preview on Azure Postgres Flexible Server
We recently announced the 𝗽𝗿𝗲𝘃𝗶𝗲𝘄 𝗼𝗳 𝗣𝗼𝘀𝘁𝗴𝗿𝗲𝗦𝗤𝗟 𝟭𝟳 on Azure Database for PostgreSQL – 𝗙𝗹𝗲𝘅𝗶𝗯𝗹𝗲 𝗦𝗲𝗿𝘃𝗲𝗿! :party_popper:
This release brings exciting new features like improved query performance, dynamic logical replication, enhanced JSON functions, and more—all backed by Azure’s reliable managed services.
Try out the preview now and share your feedback! For details, read the complete blog post👉 https://techcommunity.microsoft.com/t5/azure-database-for-postgresql/postgresql-17-preview-on-azure-postgres-flexible-server/bc-p/4263877#M474
#PostgreSQL #AzurePostgres #PGConfNYC #Database #OpenSource
We recently announced the 𝗽𝗿𝗲𝘃𝗶𝗲𝘄 𝗼𝗳 𝗣𝗼𝘀𝘁𝗴𝗿𝗲𝗦𝗤𝗟 𝟭𝟳 on Azure Database for PostgreSQL – 𝗙𝗹𝗲𝘅𝗶𝗯𝗹𝗲 𝗦𝗲𝗿𝘃𝗲𝗿! :party_popper:This release brings exciting new features like improved query performance, dynamic logical replication, enhanced JSON functions, and more—all backed by Azure’s reliable managed services.Try out the preview now and share your feedback! For details, read the complete blog post👉 https://techcommunity.microsoft.com/t5/azure-database-for-postgresql/postgresql-17-preview-on-azure-postgres-flexible-server/bc-p/4263877#M474#PostgreSQL #AzurePostgres #PGConfNYC #Database #OpenSource Read More
How to find the sensitivity of the model output if it is a variable for repeated assignment?
Hi,
I want to find the sensitivity of model output which is a tissue mass for models parameters. I received an red flag that I can’t do that for repeaded variable.
How can I handle that?
Thank you.Hi,
I want to find the sensitivity of model output which is a tissue mass for models parameters. I received an red flag that I can’t do that for repeaded variable.
How can I handle that?
Thank you. Hi,
I want to find the sensitivity of model output which is a tissue mass for models parameters. I received an red flag that I can’t do that for repeaded variable.
How can I handle that?
Thank you. sensitivity MATLAB Answers — New Questions
Error using delaunayTriangulation when calling mesh() on a pcbStack
Hi! this is the first time I ask a question here I hope im doing everthing right. Well, I’m trying to simulate the EH fields of a 2 layer pcb coil, I import the files and make the pcbStack with the following code:
P = gerberRead("gerber/6L_planar_transformerv3-F_Cu.gbr", …
"gerber/6L_planar_transformerv3-In1_Cu.gbr", …
"gerber/6L_planar_transformerv3-PTH.drl");
pb2 = pcbStack(P);
f_cu=pb2.Layers{2};
in1_cu=pb2.Layers{4};
p = pcbStack;
p.BoardThickness=0.203e-3*1+1e-3;
p.ViaDiameter=0.1e-3;
d1 = dielectric("FR4");
d1.Thickness = 0.203e-3;
air = dielectric("Air");
air.Thickness = 1e-3;
p.BoardShape = antenna.Rectangle(Length=60e-3, Width=60e-3,Center=[0.1187 -0.0775]);
%p.Layers = {f_cu,d1,in1_cu,d1,in2_cu,d1,in3_cu,d1,d1,in4_cu,d1,b_cu};
p.Layers={air,f_cu,d1,in1_cu,air}
p.FeedDiameter = 0.1e-3;
p.FeedLocations = [146.1e-3 -72e-3 2;]
p.ViaLocations = [143e-3 -77.5e-3 2 4;
117.5e-3 -77.5e-3 2 4]
show(p)
The reason I’m importing the pcbStack like that is that I’m planning on simulating a planar transformer EH fields and the gerberRead() only allows for 2 layered pcb, this way i can add more layers by importing more gerber files. However when I run
mesh(p)
After some 20 minutes of the program running, it returns the following errors:
Error using delaunayTriangulation
Enforcing Delaunay constraints leads to unstable repeated constraint intersections.
This may be due to near coincident intersecting constraints or near duplicate points.
Error in em.internal.meshprinting.inter6_full_mesh
Error in em.internal.meshprinting.multiLayerMetalImprint
Error in em.PCBStructures/generateMeshForStripFeedAndVia
Error in pcbStack/meshGenerator (line 1340)
Mesh = generateMeshForStripFeedAndVia(obj,mesherInput);
Error in em.MeshGeometry/updateMeshForPcbStack
Error in em.MeshGeometry/protectedmesh
Error in em.MeshGeometryAnalysis/mesh (line 66)
protectedmesh(obj,varargin{:});
And since some functions of the RF toolbox like current() EHFields() first call mesh() I cant run simulations of my design. I don’t know if the size could be the reason of the delaunayTriangulation failure or something is wrong with my setupHi! this is the first time I ask a question here I hope im doing everthing right. Well, I’m trying to simulate the EH fields of a 2 layer pcb coil, I import the files and make the pcbStack with the following code:
P = gerberRead("gerber/6L_planar_transformerv3-F_Cu.gbr", …
"gerber/6L_planar_transformerv3-In1_Cu.gbr", …
"gerber/6L_planar_transformerv3-PTH.drl");
pb2 = pcbStack(P);
f_cu=pb2.Layers{2};
in1_cu=pb2.Layers{4};
p = pcbStack;
p.BoardThickness=0.203e-3*1+1e-3;
p.ViaDiameter=0.1e-3;
d1 = dielectric("FR4");
d1.Thickness = 0.203e-3;
air = dielectric("Air");
air.Thickness = 1e-3;
p.BoardShape = antenna.Rectangle(Length=60e-3, Width=60e-3,Center=[0.1187 -0.0775]);
%p.Layers = {f_cu,d1,in1_cu,d1,in2_cu,d1,in3_cu,d1,d1,in4_cu,d1,b_cu};
p.Layers={air,f_cu,d1,in1_cu,air}
p.FeedDiameter = 0.1e-3;
p.FeedLocations = [146.1e-3 -72e-3 2;]
p.ViaLocations = [143e-3 -77.5e-3 2 4;
117.5e-3 -77.5e-3 2 4]
show(p)
The reason I’m importing the pcbStack like that is that I’m planning on simulating a planar transformer EH fields and the gerberRead() only allows for 2 layered pcb, this way i can add more layers by importing more gerber files. However when I run
mesh(p)
After some 20 minutes of the program running, it returns the following errors:
Error using delaunayTriangulation
Enforcing Delaunay constraints leads to unstable repeated constraint intersections.
This may be due to near coincident intersecting constraints or near duplicate points.
Error in em.internal.meshprinting.inter6_full_mesh
Error in em.internal.meshprinting.multiLayerMetalImprint
Error in em.PCBStructures/generateMeshForStripFeedAndVia
Error in pcbStack/meshGenerator (line 1340)
Mesh = generateMeshForStripFeedAndVia(obj,mesherInput);
Error in em.MeshGeometry/updateMeshForPcbStack
Error in em.MeshGeometry/protectedmesh
Error in em.MeshGeometryAnalysis/mesh (line 66)
protectedmesh(obj,varargin{:});
And since some functions of the RF toolbox like current() EHFields() first call mesh() I cant run simulations of my design. I don’t know if the size could be the reason of the delaunayTriangulation failure or something is wrong with my setup Hi! this is the first time I ask a question here I hope im doing everthing right. Well, I’m trying to simulate the EH fields of a 2 layer pcb coil, I import the files and make the pcbStack with the following code:
P = gerberRead("gerber/6L_planar_transformerv3-F_Cu.gbr", …
"gerber/6L_planar_transformerv3-In1_Cu.gbr", …
"gerber/6L_planar_transformerv3-PTH.drl");
pb2 = pcbStack(P);
f_cu=pb2.Layers{2};
in1_cu=pb2.Layers{4};
p = pcbStack;
p.BoardThickness=0.203e-3*1+1e-3;
p.ViaDiameter=0.1e-3;
d1 = dielectric("FR4");
d1.Thickness = 0.203e-3;
air = dielectric("Air");
air.Thickness = 1e-3;
p.BoardShape = antenna.Rectangle(Length=60e-3, Width=60e-3,Center=[0.1187 -0.0775]);
%p.Layers = {f_cu,d1,in1_cu,d1,in2_cu,d1,in3_cu,d1,d1,in4_cu,d1,b_cu};
p.Layers={air,f_cu,d1,in1_cu,air}
p.FeedDiameter = 0.1e-3;
p.FeedLocations = [146.1e-3 -72e-3 2;]
p.ViaLocations = [143e-3 -77.5e-3 2 4;
117.5e-3 -77.5e-3 2 4]
show(p)
The reason I’m importing the pcbStack like that is that I’m planning on simulating a planar transformer EH fields and the gerberRead() only allows for 2 layered pcb, this way i can add more layers by importing more gerber files. However when I run
mesh(p)
After some 20 minutes of the program running, it returns the following errors:
Error using delaunayTriangulation
Enforcing Delaunay constraints leads to unstable repeated constraint intersections.
This may be due to near coincident intersecting constraints or near duplicate points.
Error in em.internal.meshprinting.inter6_full_mesh
Error in em.internal.meshprinting.multiLayerMetalImprint
Error in em.PCBStructures/generateMeshForStripFeedAndVia
Error in pcbStack/meshGenerator (line 1340)
Mesh = generateMeshForStripFeedAndVia(obj,mesherInput);
Error in em.MeshGeometry/updateMeshForPcbStack
Error in em.MeshGeometry/protectedmesh
Error in em.MeshGeometryAnalysis/mesh (line 66)
protectedmesh(obj,varargin{:});
And since some functions of the RF toolbox like current() EHFields() first call mesh() I cant run simulations of my design. I don’t know if the size could be the reason of the delaunayTriangulation failure or something is wrong with my setup pcb, rf-toolbox, pcbstack, delaunaytriangulation MATLAB Answers — New Questions
set the tick format of y axis
How can I set the format of the ticks of the y axis?
I want to change the upper one to the lower one, as the former is too wide in space.How can I set the format of the ticks of the y axis?
I want to change the upper one to the lower one, as the former is too wide in space. How can I set the format of the ticks of the y axis?
I want to change the upper one to the lower one, as the former is too wide in space. plotting MATLAB Answers — New Questions
Trouble with addressing a “jagged array” in Excel VBA – subscript out of range
I have been troubled with a subscript out of range issue. At the very least I am looking to be able to address an array of arrays, apparently called a jagged array.
I have a simple example (chatGPT) that works wonderfully and when I look at the structure of arr_filteredData, the example matches my production array. When I alter the simple example to pass down the arr_filteredData from the calling SUB, I end up receiving a subscript out of range error 9.
The one difference I noticed is that the sample starts at base 0 and the one I have is starting at base 1. and mine looks like arr_filteredData (1,1) and the sample looks like arr_filteredData(0)(0). I think this is the primary issue. How should I be addressing the production array?
Unfortunately, I can’t supply data to go along with this example. The code snippet is attached.
Alter the sample routine with the production array arr_filteredData
I have been troubled with a subscript out of range issue. At the very least I am looking to be able to address an array of arrays, apparently called a jagged array.I have a simple example (chatGPT) that works wonderfully and when I look at the structure of arr_filteredData, the example matches my production array. When I alter the simple example to pass down the arr_filteredData from the calling SUB, I end up receiving a subscript out of range error 9. The one difference I noticed is that the sample starts at base 0 and the one I have is starting at base 1. and mine looks like arr_filteredData (1,1) and the sample looks like arr_filteredData(0)(0). I think this is the primary issue. How should I be addressing the production array?Unfortunately, I can’t supply data to go along with this example. The code snippet is attached. Alter the sample routine with the production array arr_filteredData Read More
CoinEX Referral Code: en3gb (Claim Sign-up Bonus)
The best CoinEX referral code for 2024 is “en3gb”. If you’re looking at how to get the CoinEX invite code, ReviewExchanges is a great place to start. By registering with this code, you’ll receive a substantial $5,000 sign-up bonus. Plus, referring other friends can grant you an impressive 50% commission, allowing you to earn up to 10,000 USDT as a welcome reward. Here’s how to use your CoinEX referral code: sign up using the “en3gb” promo code and start your crypto trading journey with extra perks while inviting your friends to enjoy CoinEX’s bonuses.
ExchangeCoineXCoineX Referral Codeen3gbTrading Fee Discount 50%Welcome BonusUp to $100000 USDTKYCNo
The best CoinEX referral code for 2024 is “en3gb”. If you’re looking at how to get the CoinEX invite code, ReviewExchanges is a great place to start. By registering with this code, you’ll receive a substantial $5,000 sign-up bonus. Plus, referring other friends can grant you an impressive 50% commission, allowing you to earn up to 10,000 USDT as a welcome reward. Here’s how to use your CoinEX referral code: sign up using the “en3gb” promo code and start your crypto trading journey with extra perks while inviting your friends to enjoy CoinEX’s bonuses.ExchangeCoineXCoineX Referral Codeen3gbTrading Fee Discount 50%Welcome BonusUp to $100000 USDTKYCNo Read More
Smoothing data without span
By how many samples will the data be smoothed if no span is used?By how many samples will the data be smoothed if no span is used? By how many samples will the data be smoothed if no span is used? smooth, span MATLAB Answers — New Questions
Bulk alter SQL column data value in MS-SQL2019
Complete command to alter SQL column data value.
Hi,
A newbie here, how to bulk update data in SQL column? For example my database name is “abcDB”
I have multiple tables in “abcDb” that start with .dbo.sun_001 to .dbo.sun_020.
In every .dbo.sun_001 to dbo.sun.020 tables there is one column name “server” that I would like to change the value from Dns02 -> Dns04 for all the row.
Any advise on a single to two command to do it in the most simplest way? Appreciate very much.
Complete command to alter SQL column data value.Hi, A newbie here, how to bulk update data in SQL column? For example my database name is “abcDB”I have multiple tables in “abcDb” that start with .dbo.sun_001 to .dbo.sun_020.In every .dbo.sun_001 to dbo.sun.020 tables there is one column name “server” that I would like to change the value from Dns02 -> Dns04 for all the row. Any advise on a single to two command to do it in the most simplest way? Appreciate very much. Read More
Explore DFT frequency accuracy, got stuck
I’m trying to analyze the frequency accuracy of the DFT as a function of DFT size assuming a rectangular
window. Using single sinusoidal waveform sampled at 5 kHz vary its frequency in small
increments to show the how the frequency accuracy changes for at least 4 different DFT length,
256, 512, 1024, and 4096.
Here is my code:
%Part 1
% MATLAB code to analyze the frequency accuracy of the DFT manually
fs = 5000; % Sampling frequency 5 kHz
f_start = 1000; % Start frequency of sinusoid (1 kHz)
f_end = 1050; % End frequency of sinusoid (1.05 kHz)
f_step = 1; % Frequency step size (small increments)
t_duration = 1; % Signal duration in seconds
t = 0:1/fs:t_duration-1/fs; % Time vector
% DFT sizes to analyze
N_dft = [256, 512, 1024, 4096];
% Prepare figure
figure;
hold on;
for N = N_dft
freq_accuracy = []; % Store accuracy results for each DFT size
for f = f_start:f_step:f_end
% Generate sinusoidal waveform
x = sin(2*pi*f*t);
rectangular_window = rectwin(N)’;
% Apply rectangular window (implicitly applied since no windowing function)
% Select only the first N samples for each DFT calculation
x_windowed = x(1:N).*rectangular_window;
% Manually compute the DFT
X = x_windowed*dftmtx(N);
% Compute the frequency axis for the current DFT
freq_axis = (0:N-1)*(fs/N);
% Find the index of the peak in the magnitude of the DFT
[~, peak_idx] = max(abs(X));
% Estimated frequency from the DFT
f_estimated = freq_axis(peak_idx);
% Store the frequency error (absolute difference between true and estimated frequency)
freq_error = abs(f_estimated – f);
freq_accuracy = [freq_accuracy, freq_error];
end
% Plot the frequency accuracy for the current DFT size
plot(f_start:f_step:f_end, freq_accuracy, ‘DisplayName’, [‘N = ‘, num2str(N)]);
end
% Add labels and legend
xlabel(‘True Frequency (Hz)’);
ylabel(‘Frequency Error (Hz)’);
title(‘Frequency Accuracy of the DFT for Different DFT Sizes (Manual DFT)’);
legend(‘show’);
grid on;
hold off;
However, the result I got is really weird.
At some frequency the estimate error between true frequency can be 3000, while most of the time it’s under 1.
However, it I substitue this line: X = x_windowed*dftmtx(N);
with this line: X = fft(x_windowed, N);
Then I got this image
It seems that this is the correct figure.
I read from this website that fft and dftmtx function can be interchangeble
https://www.mathworks.com/help/signal/ref/dftmtx.html
But my results doesn’t approve it. I want to know why.
Besides, am I applying windowing function in a correct way?
The code is as following:
rectangular_window = rectwin(N)’;
% Apply rectangular window (implicitly applied since no windowing function)
% Select only the first N samples for each DFT calculation
x_windowed = x(1:N).*rectangular_window;I’m trying to analyze the frequency accuracy of the DFT as a function of DFT size assuming a rectangular
window. Using single sinusoidal waveform sampled at 5 kHz vary its frequency in small
increments to show the how the frequency accuracy changes for at least 4 different DFT length,
256, 512, 1024, and 4096.
Here is my code:
%Part 1
% MATLAB code to analyze the frequency accuracy of the DFT manually
fs = 5000; % Sampling frequency 5 kHz
f_start = 1000; % Start frequency of sinusoid (1 kHz)
f_end = 1050; % End frequency of sinusoid (1.05 kHz)
f_step = 1; % Frequency step size (small increments)
t_duration = 1; % Signal duration in seconds
t = 0:1/fs:t_duration-1/fs; % Time vector
% DFT sizes to analyze
N_dft = [256, 512, 1024, 4096];
% Prepare figure
figure;
hold on;
for N = N_dft
freq_accuracy = []; % Store accuracy results for each DFT size
for f = f_start:f_step:f_end
% Generate sinusoidal waveform
x = sin(2*pi*f*t);
rectangular_window = rectwin(N)’;
% Apply rectangular window (implicitly applied since no windowing function)
% Select only the first N samples for each DFT calculation
x_windowed = x(1:N).*rectangular_window;
% Manually compute the DFT
X = x_windowed*dftmtx(N);
% Compute the frequency axis for the current DFT
freq_axis = (0:N-1)*(fs/N);
% Find the index of the peak in the magnitude of the DFT
[~, peak_idx] = max(abs(X));
% Estimated frequency from the DFT
f_estimated = freq_axis(peak_idx);
% Store the frequency error (absolute difference between true and estimated frequency)
freq_error = abs(f_estimated – f);
freq_accuracy = [freq_accuracy, freq_error];
end
% Plot the frequency accuracy for the current DFT size
plot(f_start:f_step:f_end, freq_accuracy, ‘DisplayName’, [‘N = ‘, num2str(N)]);
end
% Add labels and legend
xlabel(‘True Frequency (Hz)’);
ylabel(‘Frequency Error (Hz)’);
title(‘Frequency Accuracy of the DFT for Different DFT Sizes (Manual DFT)’);
legend(‘show’);
grid on;
hold off;
However, the result I got is really weird.
At some frequency the estimate error between true frequency can be 3000, while most of the time it’s under 1.
However, it I substitue this line: X = x_windowed*dftmtx(N);
with this line: X = fft(x_windowed, N);
Then I got this image
It seems that this is the correct figure.
I read from this website that fft and dftmtx function can be interchangeble
https://www.mathworks.com/help/signal/ref/dftmtx.html
But my results doesn’t approve it. I want to know why.
Besides, am I applying windowing function in a correct way?
The code is as following:
rectangular_window = rectwin(N)’;
% Apply rectangular window (implicitly applied since no windowing function)
% Select only the first N samples for each DFT calculation
x_windowed = x(1:N).*rectangular_window; I’m trying to analyze the frequency accuracy of the DFT as a function of DFT size assuming a rectangular
window. Using single sinusoidal waveform sampled at 5 kHz vary its frequency in small
increments to show the how the frequency accuracy changes for at least 4 different DFT length,
256, 512, 1024, and 4096.
Here is my code:
%Part 1
% MATLAB code to analyze the frequency accuracy of the DFT manually
fs = 5000; % Sampling frequency 5 kHz
f_start = 1000; % Start frequency of sinusoid (1 kHz)
f_end = 1050; % End frequency of sinusoid (1.05 kHz)
f_step = 1; % Frequency step size (small increments)
t_duration = 1; % Signal duration in seconds
t = 0:1/fs:t_duration-1/fs; % Time vector
% DFT sizes to analyze
N_dft = [256, 512, 1024, 4096];
% Prepare figure
figure;
hold on;
for N = N_dft
freq_accuracy = []; % Store accuracy results for each DFT size
for f = f_start:f_step:f_end
% Generate sinusoidal waveform
x = sin(2*pi*f*t);
rectangular_window = rectwin(N)’;
% Apply rectangular window (implicitly applied since no windowing function)
% Select only the first N samples for each DFT calculation
x_windowed = x(1:N).*rectangular_window;
% Manually compute the DFT
X = x_windowed*dftmtx(N);
% Compute the frequency axis for the current DFT
freq_axis = (0:N-1)*(fs/N);
% Find the index of the peak in the magnitude of the DFT
[~, peak_idx] = max(abs(X));
% Estimated frequency from the DFT
f_estimated = freq_axis(peak_idx);
% Store the frequency error (absolute difference between true and estimated frequency)
freq_error = abs(f_estimated – f);
freq_accuracy = [freq_accuracy, freq_error];
end
% Plot the frequency accuracy for the current DFT size
plot(f_start:f_step:f_end, freq_accuracy, ‘DisplayName’, [‘N = ‘, num2str(N)]);
end
% Add labels and legend
xlabel(‘True Frequency (Hz)’);
ylabel(‘Frequency Error (Hz)’);
title(‘Frequency Accuracy of the DFT for Different DFT Sizes (Manual DFT)’);
legend(‘show’);
grid on;
hold off;
However, the result I got is really weird.
At some frequency the estimate error between true frequency can be 3000, while most of the time it’s under 1.
However, it I substitue this line: X = x_windowed*dftmtx(N);
with this line: X = fft(x_windowed, N);
Then I got this image
It seems that this is the correct figure.
I read from this website that fft and dftmtx function can be interchangeble
https://www.mathworks.com/help/signal/ref/dftmtx.html
But my results doesn’t approve it. I want to know why.
Besides, am I applying windowing function in a correct way?
The code is as following:
rectangular_window = rectwin(N)’;
% Apply rectangular window (implicitly applied since no windowing function)
% Select only the first N samples for each DFT calculation
x_windowed = x(1:N).*rectangular_window; dft, discrete fourier transform, matlab, digital signal processing MATLAB Answers — New Questions
How to convert from FORTRAN TO MATLAB
Dear all,
Please am working on a project and now found a FORTRAN code to speed up my research but I having been using MATLAB and so, I need to convert the newly found FORTRAN code to MATLAB. Please kindly help me out on this conversion or give a detail explanation on how best I can convert from FORTRAN to MATLAB. The FORTRAN code is in the attachment. Thanks in anticipation!Dear all,
Please am working on a project and now found a FORTRAN code to speed up my research but I having been using MATLAB and so, I need to convert the newly found FORTRAN code to MATLAB. Please kindly help me out on this conversion or give a detail explanation on how best I can convert from FORTRAN to MATLAB. The FORTRAN code is in the attachment. Thanks in anticipation! Dear all,
Please am working on a project and now found a FORTRAN code to speed up my research but I having been using MATLAB and so, I need to convert the newly found FORTRAN code to MATLAB. Please kindly help me out on this conversion or give a detail explanation on how best I can convert from FORTRAN to MATLAB. The FORTRAN code is in the attachment. Thanks in anticipation! fortran code MATLAB Answers — New Questions
Monthly calls – not taking into consideration UTC +11
Hi
it’s great you’ve got these monthly calls, but setting aside Friday Pacific time means they’re available to those of us 11 hours ahead (Australia) means they’re hosted on a Saturday. Most of us don’t work in the office on a Saturday …
Can you please consider hosting them on a Thursday your time so that we can catch them? This would be greatly appreciated 🙂
Hi it’s great you’ve got these monthly calls, but setting aside Friday Pacific time means they’re available to those of us 11 hours ahead (Australia) means they’re hosted on a Saturday. Most of us don’t work in the office on a Saturday …Can you please consider hosting them on a Thursday your time so that we can catch them? This would be greatly appreciated 🙂 Read More
