Rocco Galati
asked on
MATLAB: Index exceeds matrix dimensions.
Hi to all,
i'm having a problem with one of my old matlab scripts.
This script worked very well on a previous version of Matlab, but now I'm running it on a different machine and on a different Matlab version (R2017b) and it gives me this error:
Please, can you try to help me?
I attached to this post all the matlab files (i renamed them with *.txt extension) and also the txt file to use to populate the vectors.
The main matlab file is 'original.m':
and the cusum.m matlab file is as the following:
moving_average.txt
periodia.txt
Kfilter_gyro_mixed.txt
Kfilter_gyro_1.txt
cusum.txt
original.txt
testxx.txt
i'm having a problem with one of my old matlab scripts.
This script worked very well on a previous version of Matlab, but now I'm running it on a different machine and on a different Matlab version (R2017b) and it gives me this error:
Index exceeds matrix dimensions.
Error in cusum (line 36)
y = x1(:,c_y) - mean(x1(:,c_y));
Error in original (line 269)
[ckc,slc]=cusum(Ic(:,2)); %% PROBLEMA
Please, can you try to help me?
I attached to this post all the matlab files (i renamed them with *.txt extension) and also the txt file to use to populate the vectors.
The main matlab file is 'original.m':
clear all, close all
addpath '/home/rocco/Documents/MATLAB/code/function'
%
% figure(1)
% set(1,'units','normalized','position',[0.142187 0.484722 0.307813 0.447222],....
% 'name','Path','numbertitle','off');
% set(gca,'FontSize',11);
% figure(2)
% set(2,'units','normalized','position',[0.457031 0.484722 0.264063 0.447917],....
% 'name','Equivalent track','numbertitle','off');
% set(gca,'FontSize',11);
%
% figure(3)
% set(3,'units','normalized','position',[0.009 0.13 0.482031 0.806944],....
% 'name','Terrain estimation','numbertitle','off');
% set(gca,'FontSize',11);
%
% figure(4)
% set(4,'units','normalized','position',[0.00351563 0.0465278 0.214453 0.893056],....
% 'name','Path','numbertitle','off');
%% Load data
% First campaign january 2015
%d=load('dati\14_01_2015\terr3.txt'); %
%d=load('dati\29_01_2015\s3.txt'); %
%d=d(1:700,:);
% second campaign September/october 2015
d=load('testxx.txt');
% Acquisition system based on ROS, this is used in the experiments
% 1: Timestamp, 2: velocit� sx (m/s), 3: velocit� dx (m/s), 4: Timestamp Encoder, 5: encoder sx, 6: encoder dx, 18: Timestamp Husky 19: dROLL (deg), 20: dPitch (deg), 21: dROLLrate (deg/s), 22: dPitchRate (deg/s), 23: dYAWRate (deg/s), 24: AccX (m/s^2), 25: AccY, 26: AccZ, 27: temperatura
% C-based implementation, not used in these experiments 1: Timestamp, 2: velocit� lineare, 3: velocit� angolare, 4: velocit� sx (m/s), 5: velocit� dx (m/s), 6: encoder sx, 7: encoder dx, 8: dROLL (deg), 9: dPitch (deg), 10: dROLLrate (deg/s), 11: dPitchRate (deg/s), 12: dYAWRate (deg/s), 13: AccX (m/s^2), 14: AccY, 15: AccZ, 16: temperatura
%% Input data
t=d(:,1);
dt=diff(d(:,1)); %sampling interval Husky
tg=d(:,18); % Timestamp gyro
dtg=diff(tg); % sampling interval gyro
roll=d(:,19); %roll angle (deg)
pitch=-d(:,20); % pitch angle (deg)
% plot(t,roll)
% hold on
% plot(t, pitch,'r')
% plot(t,'r')
% hold on
% plot(tg,'k')
vs=d(:,2); % DSP velocity left side
vd=d(:,3); % DSP velocity right side
dv=(vd-vs);
es=d(:,6); % Encoder ticks left side
ed=d(:,7); % Encoder ticks right side
% Accelerations
az=d(:,26);
az_f=moving_average(az,20);
%Include Electrical currents
t_ele=d(:,7)';
t_ele=t_ele-t_ele(1);
V_batt_raw=d(:,8)';
V_batt=sgolayfilt(V_batt_raw,4,9);
V_sx_raw=d(:,9)';
V_sx=sgolayfilt(V_sx_raw,4,41);
V_dx_raw=d(:,10)';
V_dx=sgolayfilt(V_dx_raw,4,41);
I_batt_raw=d(:,11)';
I_batt=sgolayfilt(I_batt_raw,4,41);
I_sx_raw=d(:,12)';
I_sx=sgolayfilt(I_sx_raw,1,41); I_sx=I_sx'-0.14;
I_dx_raw=d(:,13)';
I_dx=sgolayfilt(I_dx_raw,1,41); I_dx=I_dx';
t1=1; % set to 1 to discuss encoder readings
if t1==0
figure(1)
plot(vs,'r.-')
hold on
plot(vd,'b.-')
Ve=diff(es)*0.000195;
k=1; %Filtraggio encoder: ordine polinomiale
f=51; %Filtraggio encoder: frame size
%Ve=diff(es)./dt*0.000013;
W_1_f2=sgolayfilt(Ve,k,f);
W_1_f=moving_average(Ve,25);
plot(Ve,'k.-')
plot(W_1_f,'k.-')
plot(W_1_f2,'g.-')
legend('Left speed','Right speed', 'Raw Encoder', 'Moving average', 'Sgo')
ylabel('velocity (m|s)')
return
end
%% Yaw rate estimation
b=0.6; %(m), Husky width
r=-(d(:,23)); % Yaw rate from gyro (deg/sec)
t_s1=0.0667;
[x_kal measurement]=Kfilter_gyro_mixed(r,t_s1, vs, vd); %The first kinematic filter provides during straight line bias estimation, x_kal(1,:)=heading angle, x_kal(2,:)=gyro bias, and during turning motion integrates the gyro
tetag=diff(x_kal(1,:))*pi/180;
tetag=tetag'; % Husky incremental change for heading inthe sampling interval
r_f=(r-x_kal(2,:)')*pi/180;
[x_kal1]=Kfilter_gyro_1(r,t_s1); %The second kinematic filter just cleans the gyro rate of turn, x_kal1(1,:)= filtered gyro yaw rate, x_kal1(2,:)=gyro bias
% Trapezoidal integration of the gyro rate of turn, simple approach
bias=mean(r(1:40)); % estimate bias as the mean value over the first 100 samples
r=r-bias; % unbiased yaw rate
yaw_t=(r(1:end-1)+0.5*diff(r)).*(dtg*1+0.0667*0); %deg
yaw_t=[0; yaw_t]*pi/180; % yaw obtained by trapezoidal intgeration
yaw_t1=(x_kal1(1,1:end-1)'+0.5*diff(x_kal1(1,:)')).*(dtg*1+0.0667*0); %deg
yaw_t1=[0; yaw_t1]*pi/180; % yaw rate obtained by trapezoidal intgeration of the filtered gyro signal (Kfilter_gyro_1)
%% Estimate Husky pose. Here the DSP estimation is trusted more than raw encoder readings
odometry=0;
if odometry==1
% velocity integration
d_s=(vs(1:end-1)+0.5*diff(vs)).*(dt*1+0.0667*0); % m, displacemnt left side in the sampling interval
d_d=(vd(1:end-1)+0.5*diff(vd)).*(dt*1+0.0667*0); % m, displacemnt right side in the sampling interval
%% Filtering
% CF=0.3; [B,A]=butter(5,CF);
% dif_xtf=filter(B,A,dif_xt);
%% 1: Pose estimation using odmetry equations under the assumption of no slip
df=(d_s+d_d)/2; % Husky incremental displacement along Longitudinal direction
tetaf=(d_d-d_s)/b; % Husky incremental change for heading
[X0,Y0,YAW0]=odometry(0,0,-90*pi/180,df(1:end), tetaf(1:end)); % Odometry estimation
%% 2: Pose estimation using wheel encoders+gyro for heading estimation
r=-(d(:,23)); % Yaw rate from gyro (deg/sec)
t_s1=0.0667;
[x_kal measurement]=Kfilter_gyro_mixed(r,t_s1, vs, vd); %The first kinematic filter provides during straight line bias estimation, x_kal(1,:)=heading angle, x_kal(2,:)=gyro bias, and during turning motion integrates the gyro
tetag=diff(x_kal(1,:))*pi/180;
tetag=tetag'; % Husky incremental change for heading inthe sampling interval
r_f=(r-x_kal(2,:)')*pi/180;
[x_kal1]=Kfilter_gyro_1(r,t_s1); %The second kinematic filter just cleans the gyro rate of turn, x_kal1(1,:)= filtered gyro yaw rate, x_kal1(2,:)=gyro bias
% Trapezoidal integration of the gyro rate of turn, simple approach
bias=mean(r(1:80)); % estimate bias as the mean value over the first 100 samples
r=r-bias; % unbiased yaw rate
yaw_t=(r(1:end-1)+0.5*diff(r)).*(dtg*1+0.0667*0); %deg
yaw_t=[0; yaw_t]*pi/180; % yaw obtained by trapezoidal intgeration
yaw_t1=(x_kal1(1,1:end-1)'+0.5*diff(x_kal1(1,:)')).*(dtg*1+0.0667*0); %deg
yaw_t1=[0; yaw_t1]*pi/180; % yaw rate obtained by trapezoidal intgeration of the filtered gyro signal (Kfilter_gyro_1)
[Xg,Yg,YAWg]=odometry(0,0,-90*pi/180,df(1:end), tetag(1:end)); % Odometry estimation, gyro integrated by first KF
[Xgt,Ygt,YAWgt]=odometry(0,0,-90*pi/180,df(1:end), yaw_t(1:end)); % Odometry estimation, gyro integrated by trapezoidal approach
[Xgt1,Ygt1,YAWgt1]=odometry(0,0,-90*pi/180,df(1:end), yaw_t1(1:end)); % Odometry estimation, using gyro yaw rate obtained by trapezoidal intgeration of the filtered gyro signal (Kfilter_gyro_1)
%%----
%% Plotting
% fig=1;
% figure(fig)
% [ax]=plot(-Y0,X0,'k.-');
% hold on
% [axg]=plot(-Yg,Xg,'r.-');
% [axgt]=plot(-Ygt,Xgt,'b.-');
% [axgt1]=plot(-Ygt1,Xgt1,'g.-');
% hold off
% xlabel('Y(m)','FontSize',12); ylabel('X(m)','FontSize',12)
% legend('Pure odometry','Gyro odometry','Gyro odometry trapezoidal','Gyro odometry second KF')
% axis equal
% set(gca,'XDir','reverse')
%
% chg=1; % Use 1 to compare the different integration methods of gyro
% if chg==0
% figure(2)
% plot(x_kal1(1,:),'g.-')
% hold on
% plot(r,'k.-')
% hold off
%
% end
end
%--------------------------------------------------------------------
%% Equivalent track estimation through Kalman filtering, the basic equation is r*L_eq=dv
R=[1]; % Measurement noise covariance
Q=R*0.8; % Process noise covariance
P = 0.1*eye(1); %Initial estimates for error covariance
A=1.0;
%H=[mu 1];
Bu = [0];
kfinit=0;
% loop over all data
n=length(r); % solve the problem with cusum
x=zeros(n,1); %State Vector
Xp=zeros(n,1); %Projected State
Px=zeros(n,1); %Error covariance
Z=zeros(n,1);
% Parameters definition for change detection
Sv=zeros(n,1);
Zsq=zeros(n,1);
fl=zeros(n,1);
g=[];
q=0;
ni=0;
Th=0.001; % Threshold first method
Tc=0.0001; % Threshold alternative method
S1=1.0;
xp0=b; %nominal track
% Choose the yaw rate input
r=x_kal1(1,:)'*pi/180; % yaw rate filtered by the second KF
%r=(r_f);
%r=r*pi/180;
for i = 1 : n
% Kalman update
i;
H=[r(i)];
if kfinit==0
xp = xp0;
else
xp=A*x(i-1,:)' + Bu; %Projected state
Xp(i,:)=xp;
end
kfinit=1;
%-----------------------------------------------------------------------
% Change detector working in parallel with the KF
z = (dv(i) - H*xp)'; % Innovation vector, use this value to monitor the filter performance
Zc(i)=z;
g(i)=z;
cu=cumsum(g(i));
Z(i,:)=cu;
if cu>Th & i>=30
S = 2.5;
q=0;
g(1:i)=0;
Sv(i)=S;
%Q=[0.06 0.00; 0 0.0]; % It shoulkd be noted that the filter apper more
%efficient acting on S rather than Q
else
S=1; %
Sv(i)=S;
%Q=[0.03 0; 0 0.0];
end
% Alternate approach
if i>20
Ic=[(1:21)' r(i-20:i)];
[ckc,slc]=cusum(Ic(:,2)); %% PROBLEMA
[Y,I] = min(ckc(:,2));
if (ckc(end,2)-Y)>0.2
Flag(i)=1;
S1=1.1;
else
Flag(i)=0;
S1=1.0;
end
end
if abs(dv(i))>0 & abs(r(i))>10*pi/180
fl(i)=1; % Flag for track estimation step
S1=1; % Use this line to run without change detector
PP = S1*(A*P*A' + Q); % Projection of the error covariance ahead, S is a wheoghted factor, S1 for the laternate method
K = PP*H'*inv(H*PP*H'+R); % Kalman Gain
x(i,:) = (xp + K*(dv(i) - H*xp))'; % Updated state with measurement i
%x(i,1)
P = (eye(1)-K*H)*PP; % Update the error covariance
Px(i,1)= P(1,1);
pause(0.0)
else
x(i,:) = xp0;
end
end
% Save data
% save Second.txt x -ascii
% Estimate the rms error for x and y
% ux=Xgps-x(:,1);
% uy=Ygps-x(:,2);
% KF1=[x(:,1),x(:,2)];
% %save KF1.txt KF1 -ascii
% rmsux = sqrt(sum(ux.*conj(ux))/size(ux,1));
% rmsuy = sqrt(sum(uy.*conj(uy))/size(uy,1));
% rmsz = sqrt(sum(Z.*conj(Z))/size(Z,1));
%--------------------------------------------------------------------------
% Change Detection applied to all the dta set (offline approach)
% I=[(1:n)' Zc'];
% [ck,sl]=cusum(I);
% Mc=[];
% Pc=[];
% Tc=0.0008;
% for i=1:n
% [Y,I] = min(ck(10:i,2));
% if i>50 & (ck(i,2)-Y)>Tc
% Mc(i)=1;
% Pc(I+10)=1;
% Pc(I+10:n)=1;
% else
% Mc(i)=0;
% Pc(i)=0;
% end
% end
%-----------------------------------------------------------------------------
%% 4: Pose estimation using odmetry equations and equivalent track
% df=(d_s+d_d)/2; % Husky incremental displacement along Longitudinal direction
% teta_eq=(d_d-d_s)./x(1:end-1); % Husky incremental change for heading
% [X_eq,Y_eq,YAW_eq]=odometry(0,0,-90*pi/180,df(1:end), teta_eq(1:end)); % Odometry estimation
%% PLOTTING
tmin=70;
tmax=length(r_f)-70; %length(r_f); %default
%tmax=120; %Samples you want to display
ymin=0;
ymax=inf;
%--------------------------------------------------------------------------
tx=13;
fig=3;
figure(fig);
t=(1:length(x))*t_s1;
subplot(3,1,3)
plot(t(tmin:tmax),x((tmin:tmax),1),'k.-','LineWidth',1.5)
hold on
%plot(t(1:tmax),Flag(1:tmax))
% plot(dv(1:tmax),'k','LineWidth',1)
% plot(r(1:tmax),'g','LineWidth',1)
%plot(fl(1:tmax),'r','LineWidth',2); %uncomment to see the electrical currents
%plot(I_sx(1:tmax),'m','LineWidth',2)
%plot(diff(I_dx(1:tmax)),'c','LineWidth',2)
% plot(I_sx(1:tmax),'r','LineWidth',2)
% plot(I_dx(1:tmax),'b','LineWidth',2)
%plot(az_f,'m.-')
hold off
ylabel('Slip track (m)','FontSize',tx)
xlabel('Time (s)','FontSize',tx)
%legend('Equivalent track','\Delta V','Yaw rate',2);
% axis equal
h = gca
set(h,'FontSize',tx)
axis([-inf inf 0 3]);
subplot(3,1,1)
tmin=1;
tmax=10;
[ax,p1,p2] = plotyy(t(tmin:tmax),dv(tmin:tmax),t(tmin:tmax),r(tmin:tmax)*180/pi,'plot');
% plot(dv(1:tmax),'k','LineWidth',1)
% plot(r(1:tmax)*180/pi,'g','LineWidth',1)
set(p1,'Color','k', 'LineWidth',1.5), pp=0.55
set(p2,'Color',[pp pp pp], 'LineWidth',1.5)
set(ax(1),'XColor','k','YColor','k','FontSize',tx)
set(ax(2),'YColor',[pp pp pp],'FontSize',tx)
set(get(ax(1), 'Ylabel'), 'String', '\Delta V (m/s)','FontSize',tx,'Color','k');
set(get(ax(2), 'Ylabel'), 'String', 'Yaw rate (deg/s)','FontSize',tx,'Color',[pp pp pp]);
%ylabel(ax(1),'\Delta V (m/s)','Color','k') % label left y-axis
%ylabel(ax(2),'Yaw rate (deg/s)','Color',[0.7 0.7 0.7]) % label right y-axis
xlabel(ax(2),'Time (s)','FontSize',tx,'Color','k') % label x-axis
axis(ax(1),[t(tmin) t(tmax) 0 1]);
axis(ax(2),[t(tmin) t(tmax) -100 100]);
legend('Speed difference','Yaw rate')
subplot(3,1,2)
%plot(t(tmin:tmax),I_sx(tmin:tmax),'k','LineWidth',2)
plot(t(tmin:tmax),I_sx(tmin:tmax),'k','LineWidth',2)
hold on
%plot(t(tmin:tmax),I_dx(tmin:tmax),'Color',[0.6 0.6 0.6],'LineWidth',2)
plot(t(tmin:tmax),I_dx(tmin:tmax),'Color',[pp pp pp],'LineWidth',2)
ylabel('Motor current (A)','FontSize',tx)
xlabel('Time (s)','FontSize',tx)
legend('Left motor','Right motor');
% axis equal
h = gca;
set(h,'FontSize',tx)
axis([-inf inf 1 12]);
fig=2;
figure(fig);
t1=tmin; t2=tmax;
hold on
plot(t(t1:t2),roll(t1:t2),'k','LineWidth',2)
plot(t(t1:t2),pitch(t1:t2),'b','LineWidth',2)
xlabel('Time (s)','FontSize',tx); ylabel('Tilt angle (�)','FontSize',tx);
legend('Roll','Pitch');
axis([-inf inf -10 10]);
hold off
% [ax]=plot(-Y_eq,X_eq,'k.-');
% hold on
% plot(-Y0,X0,'r.-');
% % [axgt]=plot(-Ygt,Xgt,'b.-');
% % [axgt1]=plot(-Ygt1,Xgt1,'g.-');
% hold off
% xlabel('Y(m)','FontSize',12); ylabel('X(m)','FontSize',12)
% legend('Compensated odometry','Pure odometry')
% axis equal
% set(gca,'XDir','reverse')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Attitude compensation
% This part of the code aims to compensate for vehicle roll and pitch
Inl=2.9; % (A) no load current
kt=0.044; % motor constant
tau=78.7; % gear rstio
R=0.13; % wheel radius (m)
h=0.2; w=b; % cog height from the ground and vehicle width (m)
W=147*2; % Husky weigt (N)
hold on
frl=((kt*tau*2*(I_sx-Inl))/(R*W)+sind(pitch))./(cosd(pitch).*(cosd(roll)-sind(roll)*h/w)); % coefficient of resistance
frl0=((kt*tau*2*(I_sx-Inl))/(R*W)); % coefficient of resistance left side
frr=((kt*tau*2*(I_dx-Inl))/(R*W)+sind(pitch))./(cosd(pitch).*(cosd(roll)+sind(roll)*h/w)); % coefficient of resistance
frr0=((kt*tau*2*(I_dx-Inl))/(R*W)); % coefficient of resistance left side
fig=1;
figure(fig);
hold on
plot(t(t1:t2),frl(t1:t2),'k','LineWidth',2)
plot(t(t1:t2),frl0(t1:t2),'b','LineWidth',2)
plot(t(t1:t2),frr(t1:t2),'r','LineWidth',2)
plot(t(t1:t2),frr0(t1:t2),'g','LineWidth',2)
xlabel('Time (s)','FontSize',tx); ylabel('Coefficient of resistance','FontSize',tx);
legend('Left side', 'Left side (no comp)','right side', 'Right side (no comp)');
axis([-inf inf 0 2]);
hold off
%% Terrain classification
tt=1;
if tt==1
t=load('Classes.txt');
% 1: I_0
% 2: I_t (onset)
% 3: I_t (steady-state)
% 4: slip track
Inl=2.85; % I no load
kt=0.1812; %scale factor f_r
kt1=0.028; %scale factor mu_t
s=t(1:4,:);
s(3,:)=(s(3,:)-s(1,:)+Inl)*kt1; %mu_t
s(1,:)=(s(1,:)-Inl)*kt; % f_r
at=t(5:8,:);
at(3,:)=(at(3,:)-at(1,:)+Inl)*kt1; %mu_t
at(1,:)=(at(1,:)-Inl)*kt;
a=t(9:12,:);
a(3,:)=(a(3,:)-a(1,:)+Inl)*kt1; %mu_t
a(1,:)=(a(1,:)-Inl)*kt;
p=t(13:16,:);
p(3,:)=(p(3,:)-p(1,:)+Inl)*kt1; %mu_t
p(1,:)=(p(1,:)-Inl)*kt;
figure(4)
subplot(3,1,1)
plot(a(4,:), a(1,:),'sk','MarkerSize',5,'MarkerFaceColor','k')
hold on
plot(at(4,:), at(1,:),'xk','MarkerSize',10,'MarkerFaceColor','k','Linewidth',2)
plot(s(4,:), s(1,:),'ok', 'MarkerSize',5,'MarkerFaceColor','k','Linewidth',2)
plot(p(4,:), p(1,:),'^k', 'MarkerSize',5,'MarkerFaceColor','k','Linewidth',2)
% fv_asp=[1.08 2.97 8.9 0.2]; % four features are considered: Leq, I0, and I_max, dI
% fv_pav=[1.09 3.15 8.0 0.2]; % three features are considered: Leq, I0, and I_max
% fv_ter=[1.27 3.8 8.9]; % three features are considered: Leq, I0, and I_max
%
%
% plot(fv_asp(1),fv_asp(2),'r.', 'MarkerSize',6)
% hold on
%
% plot(fv_pav(1),fv_pav(2),'b.', 'MarkerSize',6)
% plot(fv_ter(1),fv_ter(2),'k.', 'MarkerSize',6)
% hold off
xlabel('Slip track (m)','FontSize',tx); ylabel('Coefficient of motion resistance','FontSize',tx)
legend('Asphalt','Dirt road', 'Sand','Plowed terrain', 2)
axis([0.6 2 -0.05 0.5])
h = gca;
set(h,'FontSize',tx)
subplot(3,1,2)
plot(a(4,:), a(3,:),'sk','MarkerSize',5,'MarkerFaceColor','k')
hold on
plot(at(4,:), at(3,:),'xk','MarkerSize',10,'MarkerFaceColor','k','Linewidth',2)
plot(s(4,:), s(3,:),'ok', 'MarkerSize',5,'MarkerFaceColor','k','Linewidth',2)
plot(p(4,:), p(3,:),'^k', 'MarkerSize',5,'MarkerFaceColor','k','Linewidth',2)
% fv_asp=[1.08 2.97 8.9 0.2]; % four features are considered: Leq, I0, and I_max, dI
% fv_pav=[1.09 3.15 8.0 0.2]; % three features are considered: Leq, I0, and I_max
% fv_ter=[1.27 3.8 8.9]; % three features are considered: Leq, I0, and I_max
%
%
% plot(fv_asp(1),fv_asp(2),'r.', 'MarkerSize',6)
% hold on
%
% plot(fv_pav(1),fv_pav(2),'b.', 'MarkerSize',6)
% plot(fv_ter(1),fv_ter(2),'k.', 'MarkerSize',6)
% hold off
xlabel('Slip track (m)','FontSize',tx); ylabel('Coefficient of lateral resistance','FontSize',tx)
legend('Asphalt','Dirt road', 'Sand','Plowed terrain', 3)
axis([0.6 2 0.1 0.3])
h = gca;
set(h,'FontSize',tx)
subplot(3,1,3)
plot(a(1,:), a(3,:),'sk','MarkerSize',5,'MarkerFaceColor','k')
hold on
plot(at(1,:), at(3,:),'xk','MarkerSize',10,'MarkerFaceColor','k','Linewidth',2)
plot(s(1,:), s(3,:),'ok', 'MarkerSize',5,'MarkerFaceColor','k','Linewidth',2)
plot(p(1,:), p(3,:),'^k', 'MarkerSize',5,'MarkerFaceColor','k','Linewidth',2)
% fv_asp=[1.08 2.97 8.9 0.2]; % four features are considered: Leq, I0, and I_max, dI
% fv_pav=[1.09 3.15 8.0 0.2]; % three features are considered: Leq, I0, and I_max
% fv_ter=[1.27 3.8 8.9]; % three features are considered: Leq, I0, and I_max
%
%
% plot(fv_asp(1),fv_asp(2),'r.', 'MarkerSize',6)
% hold on
%
% plot(fv_pav(1),fv_pav(2),'b.', 'MarkerSize',6)
% plot(fv_ter(1),fv_ter(2),'k.', 'MarkerSize',6)
% hold off
xlabel('Coefficient of motion resistance','FontSize',tx); ylabel('Coefficient of lateral resistance','FontSize',tx)
legend('Asphalt','Dirt road', 'Sand','Plowed terrain', 3)
axis([-0.02 0.5 0.1 0.3])
h = gca;
set(h,'FontSize',tx)
end
and the cusum.m matlab file is as the following:
function [ck,sl]=cusum(x1,L);
%--------------------------------------------------------------------
% cusum function cumulative sum (CUSUM) chart.
% to detect non stationarity in a series mean.
% input : - Cvenly space observations matrix [Index Value], i.e [time
% Measurement]
% - L (used to estimate the spectral density at zero frequency).
% (L need to be small compared to N).
% default is floor(0.1*N).
% - K (used to estimate the spectral density at zero frequency).
% default is 1.
% output : - [k, cumulative sum up to k].
% - False Alarm Probability.
% Reference : Koen & Lombard 1993 MNRAS 263, 287-308.
% Tested : Matlab 4.2
% By : Eran O. Ofek October 1994
% URL : http://wise-obs.tau.ac.il/~eran/matlab.html
%--------------------------------------------------------------------
%---------------------------------------------------------------------
c_x = 1;
c_y = 2;
N = length(x1(:,c_x));
if nargin==1,
L = floor(0.1.*N);
K = 1;
elseif nargin==2,
K = 1;
elseif nargin==3,
c_x = c_x;
else
error('1, 2 or 3 args only');
end
y = x1(:,c_y) - mean(x1(:,c_y));
ck = [x1(:,c_x), cumsum(y)];
% calculating the spectral density at zero frequency.
f1 = K./N;
f2 = L./N;
df = f1:(f2-f1)./(L-K):f2;
p = periodia(x1,f1,f2,df);
S0 = sum(p(:,2))./L;
% normalizing the CUSUM.
ck = ck./sqrt(N.*S0);
Dn = max(abs(ck(:,2)));
sl = 2.*exp(-2.*Dn.*Dn);
movingmean.txtmoving_average.txt
periodia.txt
Kfilter_gyro_mixed.txt
Kfilter_gyro_1.txt
cusum.txt
original.txt
testxx.txt
ASKER
Hello and thank you for your answer!
I edited my question by posting the whole cusum.m file so now the line numbers match with the line numbers in the error.
My matlab version is R2017b, but if I well remember this same code run fine on the same matlab version, but on a different machine.
Now, I tried to run the code both on Windows 10 and on Ubuntu, but I always get this error.
Hope someone can help me!
I edited my question by posting the whole cusum.m file so now the line numbers match with the line numbers in the error.
My matlab version is R2017b, but if I well remember this same code run fine on the same matlab version, but on a different machine.
Now, I tried to run the code both on Windows 10 and on Ubuntu, but I always get this error.
Hope someone can help me!
It worked on R2017b but on a different machine with same version it does not work. A little strange. A co-worker recently commented that a newer version had a different number of arguments for a function he was using. Possibly somehow on a different machine you're accessing a different function version.
In any case, you should set breakpoints at the error locations and start checking variable dimensions to get a better understanding of how far you are off.
In any case, you should set breakpoints at the error locations and start checking variable dimensions to get a better understanding of how far you are off.
ASKER
I tried to debug the script in order to check when the error occurs, but everything seems to be OK.
Anyway, I tried to run the code on my old machine few minutes ago and it doesn't work, so, maybe I'm using a different version of the script.
The matrix dimensions are OK for me, so I can't understand why I get that error.
Anyway, I tried to run the code on my old machine few minutes ago and it doesn't work, so, maybe I'm using a different version of the script.
The matrix dimensions are OK for me, so I can't understand why I get that error.
When you are in the debugger what are the sizes of the input arguments and the sizes a variables at line 36 and 37.
If no one else participates, I will try to get an old Matlab version tonight to see if your program will run on it.
In the meantime try to find your code version that you said works.
If no one else participates, I will try to get an old Matlab version tonight to see if your program will run on it.
In the meantime try to find your code version that you said works.
ASKER
On line 269 in original.m file, the statement that generates the error is:
The Ic contains:
and
and I get this when I try to print the statement on line 36:
I don't think I will able to find the old version of my code, I think I'll have to make this version work since the old one is gone :(
[ckc,slc]=cusum(Ic(:,2));
The Ic contains:
Ic =
1.0000 -0.0019
2.0000 -0.0055
3.0000 -0.0074
4.0000 0.0007
5.0000 0.0000
6.0000 -0.0025
7.0000 -0.0078
8.0000 -0.0059
9.0000 -0.0075
10.0000 -0.0027
11.0000 -0.0008
12.0000 -0.0024
13.0000 -0.0002
14.0000 -0.0030
15.0000 -0.0045
16.0000 -0.0086
17.0000 -0.0047
18.0000 -0.0018
19.0000 0.0019
20.0000 0.0040
21.0000 -0.0036
and
size(Ic)
ans =
21 2
and I get this when I try to print the statement on line 36:
x1 =
-0.0019
-0.0055
-0.0074
0.0007
0.0000
-0.0025
-0.0078
-0.0059
-0.0075
-0.0027
-0.0008
-0.0024
-0.0002
-0.0030
-0.0045
-0.0086
-0.0047
-0.0018
0.0019
0.0040
-0.0036
K>> x1(:,c_y)
Index exceeds matrix dimensions.
36 y = x1(:,c_y) - mean(x1(:,c_y));
I don't think I will able to find the old version of my code, I think I'll have to make this version work since the old one is gone :(
c_y is hardcoded to 2.
x1 is just a vector, so 2 is outside the range of the 2nd dimension.
Take a look at the size of Ic(:,2), which is the 2nd column of Ic.
x1 is just a vector, so 2 is outside the range of the 2nd dimension.
Take a look at the size of Ic(:,2), which is the 2nd column of Ic.
ASKER
If I change the instruction on line 269 to pass a vector of size 2 in this way:
then I get this error:
where perioda.m is:
I'm becoming crazy! :(
[ckc,slc]=cusum(Ic);
then I get this error:
Error using zeros
Size inputs must be scalar.
Error in periodia (line 47)
Pxw = zeros((h_f-l_f)./df,2);
Error in cusum (line 42)
p = periodia(x1,f1,f2,df);
Error in original (line 269)
[ckc,slc]=cusum(Ic);
where perioda.m is:
function [Pxw,Fm]=periodia(x,l_f,h_f,df,c_x,c_y)
% Classical power spectrum of non-evenly space time series (OBSOLETE)
% Package: timeseries
% Description: Classical power-spectrum of a non-evenly spaced time series.
% The power spectrum is normalized by the variance of the data.
% OBSOLETE: Use timeseries.period instead.
% Input : - Time series matrix, [Time, Mag], in which the first column
% is the time and the second column is the magnitude.
% - Lowest frequency (h_l).
% - Highest frequency (h_f).
% - Frequency interval (df).
% - The column of time (c_x), default is 1.
% - The column of magnitudes (c_y), default is 2.
% Output : - Periodigram matrix. normalized with the variance
% of the observations (Horne & Baliunas 1986).
% The first column is the frequency and
% the second column is the power.
% - The peaks of the periodogram sorted by the probability.
% [Frequency, Power, Period, (1-False alarm probability)].
% The probability is good only for evenly spaced data
% See Also : periodis; periodit; pdm; fitharmo; period
% Reference: Koen, C. 1990, ApJ 348, 700-702.
% Tested : Matlab 4.2
% By : Eran O. Ofek Dec 1993
% URL : http://weizmann.ac.il/home/eofek/matlab/
% Reliable: 1
%--------------------------------------------------------------------------
if nargin==4,
c_x = 1;
c_y = 2;
elseif nargin==5,
c_y = 2;
elseif nargin==6,
% do nothing
else
error('Illegal number of input arguments');
end
% set probability cutoff to zero.
pr = 0;
N0 = length(x(:,c_x));
Ni = -6.362 + 1.193.*N0 + 0.00098.*N0.*N0;
tem0 = x(:,c_y) - mean(x(:,c_y));
f_ind = l_f;
k = 1;
Pxw = zeros((h_f-l_f)./df,2);
while f_ind<h_f,
Pxw(k,1) = f_ind;
temp = abs(sum(tem0.*exp(-i*2*pi*f_ind*x(:,c_x))));
Pxw(k,2) = temp*temp/N0;
f_ind = f_ind + df;
k = k + 1;
end
% normalization of the periodogram
noise = std(x(:,c_y)).^2;
Pxw(:,2) = Pxw(:,2)./noise;
if (nargout==2),
DiffVec = diff(sign(diff([0;Pxw(:,2);0])));
IdV = find(DiffVec==-2);
Fm = [Pxw(IdV,1), Pxw(IdV,2), 1./Pxw(IdV,1), (1-exp(-Pxw(IdV,2))).^Ni];
Fm = sortrows(Fm,2);
end
I'm becoming crazy! :(
Are you saying that you solved the original error?
>> I'm becoming crazy! :(
I suggest that you choose not to become crazy. :) Instead, I suggest that you do what I did when I became stuck in matlab program. Just step through the program one line at a time, and step into functions, and when you get to the line number that will cause the error, stop and look around at all the sizes to see if that helps you figure out the error.
Another thing to get the best help (as well as help you solve your own problem) is to keep reducing the problem in a sandbox to have the minimal amount of code that causes a problem, and then ask a new question. The volunteers here like to see full minimal programs that are easy to copy and paste into their matlab, and then fix one problem for you. Sounds like you have a number of cascading problems.
>> I'm becoming crazy! :(
I suggest that you choose not to become crazy. :) Instead, I suggest that you do what I did when I became stuck in matlab program. Just step through the program one line at a time, and step into functions, and when you get to the line number that will cause the error, stop and look around at all the sizes to see if that helps you figure out the error.
Another thing to get the best help (as well as help you solve your own problem) is to keep reducing the problem in a sandbox to have the minimal amount of code that causes a problem, and then ask a new question. The volunteers here like to see full minimal programs that are easy to copy and paste into their matlab, and then fix one problem for you. Sounds like you have a number of cascading problems.
What is your level of background in Matlab, and what reference materials do you use? Is the OP code your own, and you are having trouble with the syntax? Let us know, so we can advise you better.
I got home from work too late to download all your files and see what is going on. I hope you take my advice in my previous post on how to debug Matlab (or other languages for that matter).
I got the impression that you got past your original OP question by solving the problem. It is fine to close this question and ask other question(s).
If you ask another question with a minimal program that illustrates your current question, then others who know matlab may have time during this holiday season to help.
I got home from work too late to download all your files and see what is going on. I hope you take my advice in my previous post on how to debug Matlab (or other languages for that matter).
I got the impression that you got past your original OP question by solving the problem. It is fine to close this question and ask other question(s).
If you ask another question with a minimal program that illustrates your current question, then others who know matlab may have time during this holiday season to help.
ASKER CERTIFIED SOLUTION
membership
This solution is only available to members.
To access this solution, you must be a member of Experts Exchange.
ASKER
Hwllo phoffric, sorry for my delay, but I spent all yesterday to check my code (I didn't rember anything about it since I wrote it several years ago) and now it works fine!
I checked all the vectors and their size and I found out that there were an error with a function.
Thank you for your assistance!
P.S. The code runs an algorithms which read different sensors parameters and feed them into a Kalman Filter.
I checked all the vectors and their size and I found out that there were an error with a function.
Thank you for your assistance!
P.S. The code runs an algorithms which read different sensors parameters and feed them into a Kalman Filter.
ASKER
Thank you!
Hey Rocco,
Good to hear that you solved the remaining problems.
Happy New Year!
Paul
Good to hear that you solved the remaining problems.
Happy New Year!
Paul
In your posted code the line numbers do not appear to agree the line numbers that you report . Did you post a different version of the two files , original.m and cusum.m?
For example: Error in cusum (line 36) but the cusum.m code that you posted only shows 27 lines.