Link to home
Start Free TrialLog in
Avatar of Rocco Galati
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:

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 

Open in new window


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

Open in new window


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);

Open in new window

movingmean.txt
moving_average.txt
periodia.txt
Kfilter_gyro_mixed.txt
Kfilter_gyro_1.txt
cusum.txt
original.txt
testxx.txt
Avatar of phoffric
phoffric

In case someone has your two versions of Matlab , they may be able to reproduce the problem better if you specify which version works and which version gives the errors.

 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.
Avatar of Rocco Galati

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!
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.
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.
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.
On line 269 in original.m file, the statement that generates the error is:
  [ckc,slc]=cusum(Ic(:,2));

Open in new window


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

Open in new window


and

size(Ic)
ans =
    21     2

Open in new window


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));

Open in new window


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.
If I change the instruction on line 269 to pass a vector of size 2 in this way:

 [ckc,slc]=cusum(Ic);

Open in new window


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);

Open in new window


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

Open in new window


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.
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.
ASKER CERTIFIED SOLUTION
Avatar of phoffric
phoffric

Link to home
membership
This solution is only available to members.
To access this solution, you must be a member of Experts Exchange.
Start Free Trial
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.
Thank you!
Hey Rocco,
Good to hear that you solved the remaining problems.
Happy New Year!
Paul