Thinking about it, I suppose there are two problems that I am trying to solve: First, how to extract from the data segments of when the temperature spikes and returns to baseline. Second, do a linear fit for each segment to calculate the rate of temperature decrease.
How to find rates of return to baseline after peaks, in temperature data over long timecourse?
82 views (last 30 days)
Show older comments
Hello all, new MatLab user here, find the tutorials and documentation very useful so far! I am trying to evaluate the performance of a freezer in my lab. I wish to know the rate at which it cools back down to its setpoint after a spike in temperature. I want to be able to determine these rates for events over two years
I've been able to import the data into matlab as a table with the datetime and temperature, and then use the "find local extrema" task to identify the peaks. But I am now stuck on how I can find the rates of change after these peaks in an automated fashion. I also converted the data into timetables and tried doing this in the Signal Analyzer app. Would anyone be able to point me in the right direction?
I've included a plot to provide an idea of the data I'm working with. This plot has ~250,000 data points
4 Comments
Star Strider
on 24 Nov 2025 at 19:01
The rate of cooling would likely be described as a decaying exponential (Newton's Law of Cooling), and that function would be:
nlc = @(b,t) b(1) .* exp(b(2) .*t) + b(3);
where b(1) is the original temperature excursion, b(2) is the rate you probably want to get, and b(3) is the offset. You would likely need to use the minutes or seconds function on the time variable (datetime turned into a duration) and then determine the sections of the data you want to use for each nonlinear regression. There are several funcitons that can handle that sort of regression in the Optimization Toolbox and the Statistics and Machine Learning Toolbox. The fminsearch function in core MATLAB could also work. It depends on what you have available.
Answers (2)
dpb
on 24 Nov 2025 at 20:24
Edited: dpb
on 25 Nov 2025 at 6:09
whos -file EventSelection
load EventSelection
head(EventSection)
d=dir('*.csv');
l=readlines(d.name);
l(1:20)
EventSectionPlot
t=EventSectionPlot.XData; T=EventSectionPlot.YData;
plot(t,T)
Very astute observation by @Star Strider -- with the event log, I think I would also read it and select the pairs of door open and close events; those will be the start of a time from which to find the subsequent peak after the door was opened and then can look for the time at which the temperature again is a minimum before another rise.
The present table of data doesn't have the temperature trace; that would be the cleaner one to work from rather than the plot data...
Let's find one peak to play with and see what actually looks like...
[Tmax,imax]=maxk(T,1) % just largest location to look at
Npts=500; % this was empirical to see the whole peak width
plot(t(imax-Npts:imax+Npts),T(imax-Npts:imax+Npts))
hold on
N1=0; N2=300;
plot(t(imax-N1:imax+N2),T(imax-N1:imax+N2),'r-')
nlc=@(b,t) b(1).*exp(-b(2).*t)+b(3); % note -b(2) so rate constant will be positive value
tevent=minutes(t(imax:imax+N2)-t(imax)).'; % convert to elapsed time
y=T(imax:imax+N2).';
b0=[15,0.01,-80];
f=fitnlm(tevent,y,nlc,b0)
ypred=predict(f,tevent);
plot(t(imax-N1:imax+N2),ypred,'b.')
The cycling before the event does appear as @Star Strider would suggest but the recovery after isn't so smooth as well as having the little double peak...but, the model doesn't do such a bad job of the trend although it does overpredict the peak. That may not be as signficant for the purpose of estimating the cooldown rate, however; that would probably depend upon just what is the intended or needed use.
The rate constant of b(2) = 0.0206 is in inverse minutes.
I think it would be very interesting to do a number of the cycling cycles and see how consistent this is across them and the other events. If had the event log data, could put it on the plot as well...
Note this is the largest peak in the lot there in November given the way maxk works to find N largest peaks and asked for one -- so actually is same as if had used max, anyways...
5 Comments
dpb
on 25 Nov 2025 at 19:05
I am not sure what the end objective really is, so concentrated on looking at one event in some detail over mass production.
A little more background might lead to some other ideas on processing.
For example, is the double peak in this particular incident significant on its own as a separate event? As noted, the model being a least squares fit minimizes the overall error and in this particular case, anyway, does overshoot the origin. Is the predicted -64.6 versus observed -66.9 at the peak unacceptable? If it is mandatory to go through the peak point, then one would have to switch to using a method such as fmincon that would allow for introducing constraints.
On the automatic location of endpoints for excursions, using the difference and looking for a change in sign is often very useful...if I get a chance I'll give that a go as well.
dpb
on 26 Nov 2025 at 19:00
Edited: dpb
on 27 Nov 2025 at 15:02
To explore alternate ways to find end of segment, beginning of next let's expand our view and split the overall and the change into two...
load EventSelection
t=EventSection.DateTime; T=EventSection.RTD; % some shorter names
[Tmax,imax]=maxk(T,1); % just largest location to look at
Npts=500; % this was empirical to see the whole peak width
N1=50; N2=300;
hAx=subplot(3,1,1);
plot(t(imax-N1:imax+Npts),T(imax-N1:imax+Npts))
hold on
plot(t(imax:imax+N2),T(imax:imax+N2),'r-')
% now look at the temperature gradient...
hAx(2)=subplot(3,1,2);
dT=round([nan;(diff(T(imax:end)))],1);
plot(t(imax:end),dT,'x-')
xlim(hAx(2),xlim(hAx(1))), ylim(0.5*[0 1])
[U,~,iU]=unique(dT(dT>0)); U
hAx(3)=subplot(3,1,3);
dT=2*(dT>0)-1; % turn into bi-state
plot(t(imax:end),dT,'x-')
xlim(hAx(3),xlim(hAx(1)))
Now the change is very obvious; applying findpeaks here and walking back would lead to the point just before the next transient, even for the normal cycling if those are also of importance.
Alternatively, as the third plot shows, turning the difference into a bi-level flag then when
ixPlus=(diff(dT)==2);
since the change in state from -1 to +1 is now identically 2.
Star Strider
on 25 Nov 2025 at 2:53
Edited: Torsten
on 26 Nov 2025 at 15:54
This calculates all the cooling rates that my findpeaks call can identify, then plots a few of them at the end.
I did not examine its results closely, however it appears to work. (It took a few hours and a fair amount of experimentation to get it to run as I want it to.) It identifies and estimates the parameters of 345 peaks, and produces a summary table of their times and rates at the end.
Try something like this --
LD = load('EventSelection.mat')
Events1 = LD.EventSection
Events2 = readtable('Meyer Lab -80_170564_EventLog_11_18_2025_01_08_45.csv', VariableNamingRule='preserve')
[pks,plocs] = findpeaks(Events1.(2), MinPeakProminence=0.5);
for k = 1:numel(plocs)-1
[minv,idx] = min(Events1.(2)(plocs(k) : plocs(k+1)));
idxrng{k,:} = plocs(k):(plocs(k)+idx);
end
figure
plot(Events1.(1), Events1.(2))
hold on
plot(Events1.(1)(plocs), pks, '^r')
hold off
grid
Time = minutes(Events1.(1) - Events1.(1)(1));
nlc = @(b,t) b(1) .* exp(b(2) .*t) + b(3);
% locs = [locs; size(Events1,1)];
%
% QQQ = numel(idxrng)
B = NaN(3,numel(idxrng));
rn = NaN(size(B(1,:)));
for k = 1:numel(idxrng)
%k
% idxrng = locs(k) : locs(k+1);
% QQ = Events1(idxrng,:)
if isempty(idxrng{k}) | (min(Events1.(2)(idxrng{k})) == max(Events1.(2)(idxrng{k})))
continue
end
[tmin,tmax] = bounds(Time(idxrng{k}));
timek{k} = Time(idxrng{k})-tmin;
B0 = [max(Events1.(2)(idxrng{k})); -log(min(Events1.(2)(idxrng{k}))/max(Events1.(2)(idxrng{k}))) / (tmax-tmin); min(Events1.(2)(idxrng{k}))];
[B(:,k),rn(k)] = lsqcurvefit(nlc, B0, timek{k}, Events1.(2)(idxrng{k}),[],[],[],[],[],[],[],optimset('Display','none'));
% idxrngc{k} = idxrng;
tminmax(k,:) = [tmin tmax];
tminmaxidx(k,:) = [idxrng{k}(1) idxrng{k}(end)];
end
disp(B)
idxv = floor(linspace(1, size(B,2), 16))
figure
tiledlayout(4,4)
for k = 1:numel(idxv)
nexttile
plot(Time(idxrng{k}), Events1.(2)(idxrng{k}))
hold on
plot(Time(idxrng{k}), nlc(B(:,k), timek{k}), '-r')
hold off
grid
title("Curve #"+idxv(k),"Rate = "+B(2,idxv(k)))
end
sgtitle("Selected Sample")
Cooling_Rates = B(2,:)
Cooling_Rates_Table = table(Events1.(1)(tminmaxidx(:,1)), Events1.(1)(tminmaxidx(:,2)), B(2,:).', VariableNames=["Start Time","End Time","Rate"])
.
2 Comments
Star Strider
on 25 Nov 2025 at 18:03
My pleasure!
I will be glad to help with any necessary modifications.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)