Vectorization for loops in MATLAB

Question

I'm not sure if this is possible, but my understanding of MATLAB could certainly be better.

I have some code that I want to highlight, since it causes a rather bottleneck in my program. This is part of the optimization procedure, which has many possible configurations of short-term average (STA), long-term average (LTA) and sensitivity (OnSense) for passing.

Time in vector format, FL2onSS is the main data (double Nx1), FL2onSSSTA is STA (NxSTA double), FL2onSSThresh is the Threshold value (NxLTAxOnSense double)

The idea is to calculate the red alarm matrix, which will be 4D - alarmStatexSTAxLTAxOnSense, which is used in the rest of the program.

Red = zeros(length(FL2onSS), length(STA), length(LTA), length(OnSense), 'double'); for i=1:length(STA) for j=1:length(LTA) for k=1:length(OnSense) Red(:,i,j,k) = calcRedAlarm(Time, FL2onSS, FL2onSSSTA(:,i), FL2onSSThresh(:,j,k)); end end end 

Currently, I am getting this recurring function, trying to get a little more speed from it, but obviously it would be better if all of this could be vectorized. In other words, I do not need to save the function if there is a better solution.

 function [Red] = calcRedAlarm(Time, FL2onSS, FL2onSSSTA, FL2onSSThresh) % Calculate Alarms % Alarm triggers when STA > Threshold zeroSize = length(FL2onSS); %Precompose Red = zeros(zeroSize, 1, 'double'); for i=2:zeroSize %Because of time chunks being butted up against each other, alarms can %go off when they shouldn't. To fix this, timeDiff has been %calculated to check if the last date is different to the current by 5 %seconds. If it isn't, don't generate an alarm as there is either a %validity or time gap. timeDiff = etime(Time(i,:), Time(i-1,:)); if FL2onSSSTA(i) > FL2onSSThresh(i) && FL2onSSThresh(i) ~= 0 && timeDiff == 5 %If Short Term Avg is > Threshold, Trigger Red(i) = 1; elseif FL2onSSSTA(i) < FL2onSSThresh(i) && FL2onSSThresh(i) ~= 0 && timeDiff == 5 %If Short Term Avg is < Threshold, Turn off Red(i) = 0; else %Otherwise keep current state Red(i) = Red(i-1); end end end 

The code is simple enough, so I will not explain it. If you need to find out what a particular line does, let me know.

Answer 1

Focus in brings all your data into the same form , using mostly repmat and move . Then logic is the easy part.

I needed an unpleasant trick to implement the last part (if none of the conditions are met, use the latest results). usually this logic is done using cumsum. I had to use a different matrix 2. ^ n to make sure that the values ​​that are defined (so that + 1, + 1, -1 will really give 1,1,0) - just look at the code :)

 %// define size variables for better readability N = length(Time); M = length(STA); O = length(LTA); P = length(OnSense); %// transform the main data to same dimentions (3d matrices) %// note that I flatten FL2onSSThresh to be 2D first, to make things simpler. %// anyway you don't use the fact that its 3D except traversing it. FL2onSSThresh2 = reshape(FL2onSSThresh, [N, O*P]); FL2onSSThresh3 = repmat(FL2onSSThresh2, [1, 1, M]); FL2onSSSTA3 = permute(repmat(FL2onSSSTA, [1, 1, O*P]), [1, 3, 2]); timeDiff = diff(datenum(Time))*24*60*60; timeDiff3 = repmat(timeDiff, [1, O*P, M]); %// we also remove the 1st plain from each of the matrices (the vector equiv of running i=2:zeroSize FL2onSSThresh3 = FL2onSSThresh3(2:end, :, :); FL2onSSSTA3 = FL2onSSSTA3(2:end, :, :); Red3 = zeros(N-1, O*P, M, 'double'); %// now the logic in vector form %// note the chage of && (logical operator) to & (binary operator) Red3((FL2onSSSTA3 > FL2onSSThresh3) & (FL2onSSThresh3 ~= 0) & (timeDiff3 == 5)) = 1; Red3((FL2onSSSTA3 < FL2onSSThresh3) & (FL2onSSThresh3 ~= 0) & (timeDiff3 == 5)) = -1; %// now you have a matrix with +1 where alarm should start, and -1 where it should end. %// add the 0s at the begining Red3 = [zeros(1, O*P, M); Red3]; %// reshape back to the same shape Red2 = reshape(Red3, [N, O, P, M]); Red2 = permute(Red2, [1, 4, 2, 3]); %// and now some nasty trick to convert the start/end data to 1 where alarm is on, and 0 where it is off. Weights = 2.^repmat((1:N)', [1, M, O, P]); %// ' damn SO syntax highlighting. learn MATLAB already! Red = (sign(cumsum(Weights.*Red2))+1)==2; %// and we are done. %// print sum(Red(:)!=OldRed(:)), where OldRed is Red calculated in non vector form to test this.