Find automatically data of interest

Hello everyone,

I havew a question about how to improve my code.

I have a time-Temperature profile that is a bit noisy: my goal is to automatically find the heating, cooling and isotheramal portions by only looking at the first derivative of the temperature with respect to the time. Unluckily, the code is not perfect but I don't know how to improve it: can anyone help me?

Attached you can fine the test data and the code I use together with the subroutines.

close all
clear
tic

load 'conv_ageing_atm5_m60p80_01to2048s_bif_74165.txt'
conv_ageing_atm5_m60p80_01to2048s_bif_74165;
t=conv_ageing_atm5_m60p80_01to2048s_bif_74165(1:end,2);
T=conv_ageing_atm5_m60p80_01to2048s_bif_74165(1:end,3);
Cp=conv_ageing_atm5_m60p80_01to2048s_bif_74165(1:end,4);
time=conv_ageing_atm5_m60p80_01to2048s_bif_74165(1:end,5);
un=[1:1:length(T)];


ratew=(gradient(T)./gradient(time));
rate=smoothdata(ratew,'movmedian',30);
rate=nearest(rate);

flag=zeros(length(T),1);

for i=1:length(T)
   if rate(i)>0;
       flag(i)=1;
   else if rate(i)==0;
           flag(i)=0;
       else if rate(i)<0;
               flag(i)=-1;
           end
       end
   end    
end

[y]=remove_outliers(flag);

figure(1)
plot(un,flag,'o')
hold on
plot(un,y,'ro')

[Aw,Bw,Cw]=group_data(y);%Aw=heating/Bw=isother/Cw=cooling

[Vw]=organize_vector(Aw);%heating
[Uw]=organize_vector(Bw);%isotherm
[Ww]=organize_vector(Cw);%cooling


formatSpec = 'You have %d heating scans\n';
fprintf(formatSpec,length(Vw))
formatSpec2 = 'You have %d isotherms\n';
fprintf(formatSpec2,length(Uw))
formatSpec3 = 'You have %d cooling scans\n';
fprintf(formatSpec3,length(Ww))

toc

[dpb]

for i=1:length(T)
    if rate(i)>0;
        flag(i)=1;
    else if rate(i)==0;
            flag(i)=0;
    else if rate(i)<0;
            flag(i)=-1;
    end
    end
    end    
end

Use elseif command instead of nesting another if inside. The above is correctly written as

for i=1:length(T)
  if rate(i)>0;
    flag(i)=1;
  elseif rate(i)==0;
    flag(i)=0;
  elseif rate(i)<0;
    flag(i)=-1;
  end
end

Then, instead of if...elseif...end use

flag=sign(rate);
[y]=remove_outliers(flag);

Once you've converted to sign alone, there will be no outliers. If you're going to do that, probably would need to do it at least before there and most likely before the smoothing/filtering that will also have knocked any real outliers down/out.

[Jan]

Please mention, what you consider as "the code is not perfect".

[Daniele Sonaglioni]

Dear @dpb, the suggestion about sign is a real improvement to the code in term of efficiency but there are still some outliers.

[Daniele Sonaglioni]

Dear @Jan,

when I wrote "the code is not perfect" I was referring to the fact that it detects some heating and/or cooling and/or isotherm segments that do not exist.

On the data that I have attached, there are 14 heating segments, 40 isotherms and 25 cooling but the code finds more heating and isotherm scans.

By increasing the smoothing window I can solve the problem but I was wondering if there is an alternative way to solve the problem.

[dpb]

Might want to look at findchangepts in Signal Processing TB...

Answer 1

[Jan]

Some methods are "ugly":

load 'conv_ageing_atm5_m60p80_01to2048s_bif_74165.txt'

Loading data directly into the workspace creates variables dynamically. This impedes Matlab's JIT acceleration and the debugging. Prefer to load data into a specific variable:

Data = load('conv_ageing_atm5_m60p80_01to2048s_bif_74165.txt');

This line is a waste of time only:

conv_ageing_atm5_m60p80_01to2048s_bif_74165;

Use a less painful name for variables. Hiding important information in the name of variables is a programming anti-pattern:

value = Data.conv_ageing_atm5_m60p80_01to2048s_bif_74165;

Use ":" instead of "1:end", because it is faster:

% Replace: t = conv_ageing_atm5_m60p80_01to2048s_bif_74165(1:end,2); by
t = value(:, 2);

Replace the code:

flag=zeros(length(T),1);
for i=1:length(T)
   if rate(i)>0;
       flag(i)=1;
   else if rate(i)==0;
           flag(i)=0;
       else if rate(i)<0;
               flag(i)=-1;
           end
       end
   end    
end

by

flag = sign(rate);

[] is Matlab's operator for a concatenation of arrays. In

un=[1:1:length(T)];

you concatenate the vector 1:length(T) with nothing. Cleaner and faster:

un = 1:length(T);

[Daniele Sonaglioni]

Thank you for your suggestion, I will encapsulate in my code!