How to resolve dimension error?
7 views (last 30 days)
Show older comments
%% ========================== COMPUTE & PLOT IMU vs GNSS ==========================
clear; clc; close all;
disp("=== Loading Processed Data ===");
load('ECI_Data.mat'); % GNSS Data in ECI
load('IMU_Data_ECI.mat'); % IMU Data converted to ECI
load('Synthetic_GNSS_Data.mat'); % Original GNSS Data
load('Synthetic_IMU_Data.mat'); % Original IMU Data
%% -------------------- Part A: Compute Velocity & Position for IMU --------------------
disp("=== Computing IMU Velocity and Position in ECI ===");
fs_imu = 400; % IMU Sampling Frequency (Hz)
dt_imu = 1/fs_imu; % IMU time step
% Extract IMU accelerations in ECI
a_x_eci = IMU_Data_ECI.Accel_X_ECI;
a_y_eci = IMU_Data_ECI.Accel_Y_ECI;
a_z_eci = IMU_Data_ECI.Accel_Z_ECI;
% Integrate Acceleration to Compute Velocity (Trapezoidal Rule)
v_x_eci = cumtrapz(dt_imu, a_x_eci);
v_y_eci = cumtrapz(dt_imu, a_y_eci);
v_z_eci = cumtrapz(dt_imu, a_z_eci);
% Integrate Velocity to Compute Position
p_x_eci = cumtrapz(dt_imu, v_x_eci);
p_y_eci = cumtrapz(dt_imu, v_y_eci);
p_z_eci = cumtrapz(dt_imu, v_z_eci);
% Store computed IMU ECI Data
IMU_ECI.Position_X = p_x_eci;
IMU_ECI.Position_Y = p_y_eci;
IMU_ECI.Position_Z = p_z_eci;
IMU_ECI.Velocity_X = v_x_eci;
IMU_ECI.Velocity_Y = v_y_eci;
IMU_ECI.Velocity_Z = v_z_eci;
IMU_ECI.Accel_X = a_x_eci;
IMU_ECI.Accel_Y = a_y_eci;
IMU_ECI.Accel_Z = a_z_eci;
disp("IMU Position, Velocity, and Acceleration computed in ECI Frame.");
%% -------------------- Part B: Convert ECI Data Back to ECEF --------------------
disp("=== Converting ECI Data Back to ECEF for Verification ===");
omega_e = 7.2921150e-5; % Earth's Rotation Rate (rad/s)
t_imu = (0:length(p_x_eci)-1) * dt_imu; % Time vector for IMU
% Convert back to ECEF (rotation correction)
theta = omega_e * t_imu;
p_x_ecef = p_x_eci .* cos(theta) + p_y_eci .* sin(theta);
p_y_ecef = -p_x_eci .* sin(theta) + p_y_eci .* cos(theta);
p_z_ecef = p_z_eci; % No rotation in Z
v_x_ecef = v_x_eci .* cos(theta) + v_y_eci .* sin(theta);
v_y_ecef = -v_x_eci .* sin(theta) + v_y_eci .* cos(theta);
v_z_ecef = v_z_eci; % No rotation in Z
disp("Converted IMU ECI Data Back to ECEF Frame.");
%% -------------------- Part C: Compute GNSS Position & Velocity in ECI --------------------
disp("=== Interpolating GNSS Data for IMU Comparison ===");
fs_gnss = 10; % GNSS Sampling Frequency (Hz)
t_gnss = GNSS_Data.Posix_Time(:); % GNSS time vector
% Extract GNSS Positions in ECEF
p_x_gnss = GNSS_Data.X_ECEF_m(:);
p_y_gnss = GNSS_Data.Y_ECEF_m(:);
p_z_gnss = GNSS_Data.Z_ECEF_m(:);
% Convert GNSS ECEF Data to ECI (same method as before)
theta_gnss = omega_e * t_gnss;
p_x_eci_gnss = p_x_gnss .* cos(theta_gnss) + p_y_gnss .* sin(theta_gnss);
p_y_eci_gnss = -p_x_gnss .* sin(theta_gnss) + p_y_gnss .* cos(theta_gnss);
p_z_eci_gnss = p_z_gnss;
% Compute GNSS Velocity using Central Difference Approximation
v_x_eci_gnss = gradient(p_x_eci_gnss, t_gnss);
v_y_eci_gnss = gradient(p_y_eci_gnss, t_gnss);
v_z_eci_gnss = gradient(p_z_eci_gnss, t_gnss);
disp("GNSS Position and Velocity computed in ECI Frame.");
%% -------------------- Part D: Interpolate GNSS to IMU Rate --------------------
disp("=== Interpolating GNSS Data to IMU Sampling Rate ===");
p_x_eci_gnss_interp = interp1(t_gnss, p_x_eci_gnss, t_imu, 'linear', 'extrap');
p_y_eci_gnss_interp = interp1(t_gnss, p_y_eci_gnss, t_imu, 'linear', 'extrap');
p_z_eci_gnss_interp = interp1(t_gnss, p_z_eci_gnss, t_imu, 'linear', 'extrap');
v_x_eci_gnss_interp = interp1(t_gnss, v_x_eci_gnss, t_imu, 'linear', 'extrap');
v_y_eci_gnss_interp = interp1(t_gnss, v_y_eci_gnss, t_imu, 'linear', 'extrap');
v_z_eci_gnss_interp = interp1(t_gnss, v_z_eci_gnss, t_imu, 'linear', 'extrap');
disp("GNSS Data Interpolated to IMU Sampling Rate.");
%% -------------------- Part E: Plot Position, Velocity & Acceleration --------------------
disp("=== Plotting IMU vs GNSS Comparisons ===");
figure('Name','Position Comparison in ECI','NumberTitle','off');
subplot(3,1,1);
plot(t_imu, p_x_eci_gnss_interp, 'g-','LineWidth',1.5); hold on; grid on;
plot(t_imu, p_x_eci, 'b--','LineWidth',2);
xlabel('Time (s)'); ylabel('X (m)'); legend('GNSS','IMU');
title('X Position (ECI)'); axis tight;
subplot(3,1,2);
plot(t_imu, p_y_eci_gnss_interp, 'g-','LineWidth',1.5); hold on; grid on;
plot(t_imu, p_y_eci, 'b--','LineWidth',2);
xlabel('Time (s)'); ylabel('Y (m)'); legend('GNSS','IMU');
title('Y Position (ECI)'); axis tight;
subplot(3,1,3);
plot(t_imu, p_z_eci_gnss_interp, 'g-','LineWidth',1.5); hold on; grid on;
plot(t_imu, p_z_eci, 'b--','LineWidth',2);
xlabel('Time (s)'); ylabel('Z (m)'); legend('GNSS','IMU');
title('Z Position (ECI)'); axis tight;
figure('Name','Velocity Comparison in ECI','NumberTitle','off');
subplot(3,1,1);
plot(t_imu, v_x_eci_gnss_interp, 'g-','LineWidth',1.5); hold on; grid on;
plot(t_imu, v_x_eci, 'b--','LineWidth',2);
xlabel('Time (s)'); ylabel('X (m/s)'); legend('GNSS','IMU');
title('X Velocity (ECI)'); axis tight;
subplot(3,1,2);
plot(t_imu, v_y_eci_gnss_interp, 'g-','LineWidth',1.5); hold on; grid on;
plot(t_imu, v_y_eci, 'b--','LineWidth',2);
xlabel('Time (s)'); ylabel('Y (m/s)'); legend('GNSS','IMU');
title('Y Velocity (ECI)'); axis tight;
subplot(3,1,3);
plot(t_imu, v_z_eci_gnss_interp, 'g-','LineWidth',1.5); hold on; grid on;
plot(t_imu, v_z_eci, 'b--','LineWidth',2);
xlabel('Time (s)'); ylabel('Z (m/s)'); legend('GNSS','IMU');
title('Z Velocity (ECI)'); axis tight;
disp("=== All Plots Generated Successfully ===");
The issue with limitations of the data size so how can we resolve it ?
=== Loading Processed Data ===
=== Computing IMU Velocity and Position in ECI ===
IMU Position, Velocity, and Acceleration computed in ECI Frame.
=== Converting ECI Data Back to ECEF for Verification ===
Requested 48000x48000 (17.2GB) array exceeds maximum array size preference (15.9GB). This might
cause MATLAB to become unresponsive.
Related documentation
1 Comment
Matt J
on 19 Feb 2025
The message says your code attempted to build a 48000x48000 matrix. Is that the intended size of the result?
Answers (1)
Harald
on 20 Feb 2025
Hi,
To follow up on Matt's comment: if you intended element-wise multiplication to obtain a vector with 48000 elements, the problem is that you are mixing row and column vectors. You can transpose either the row vector theta or both column vectors p_x_eci and p_y_eci to work with vectors of matching orientation.
Best wishes,
Harald
0 Comments
See Also
Categories
Find more on Control System Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!