Correcting effects of Humidity on sensors

Hi All

I have gas sensor, that gets effected by hummdity that needs to be corrected. So was hoping to see if we can correct this ?

How can i run my code on the this support forum with my data file, so it can be run?

Answer 1

[Star Strider]

What sort of correction do you want to do to your data?

Do you also have the humidity data?

Are there any published ways to correct the readings for humidity? If so, please share them.

To run your code with your data here, first upload the data file, using the 'paperclip' icon in the top toolbar (just to the right of the Σ). Click on the 'insert a line of doce' icon in the top toolbar (farthest left icon in the CODE section, or ALT+ENTER) to create a code line, then type or copy-paste your code in it. To run it, press the green arrrow in the top toolbar.

x = linspace(0, 2*pi);
y = sin(x) .* cos(x);

figure
plot(x, y)
grid

Your code should run here essentially the same way it runs on your computer, including reading the file.

,

[Dharmesh]

Thank you for your reply.

As humidity increases and decreases — or due to the rate of change — it introduces transients into my analogue signal. This behaviour is known to occur in electrochemical sensors. I would like to explore whether this effect can be corrected, possibly using a quadratic or adaptive algorithm that could later be implemented or ported into other software code.

If you run my code with my dataset, you’ll notice a clear correlation between humidity and certain variations in the signal. The signal is also affected by temperature, but that influence is much easier to handle, so my current focus is on correcting the humidity-related effects.

%% AAN-803-05 corrections: kT for NO, nT for NO2 & OX (O3)
% Dharmesh – edit the folder if needed:
folder  = 'G:\AVR_Project\MATLAB Projects\Air Sensor\Temp_test';
matFile = fullfile(folder, 'SensorLog_ALL.mat');
% --- Load the combined table ---
S = load(matFile);
fn = fieldnames(S);
T  = S.(fn{find(structfun(@istable, S), 1)});
% --- Time & environment ---
if ~isdatetime(T.Timestamp)
    T.Timestamp = datetime(T.Timestamp,'InputFormat','yyyy-MM-dd HH:mm:ss');
end
time  = T.Timestamp;
TempC = T.("Temperature (°C)");
Hum   = T.("Humidity (%)");
% ================== CALIBRATION CONSTANTS ==================
% ---------- ELECTRONICS OFFSETS (hardware bias; set if measured) ----------
ELEC.NO_WE   = 0.292;  ELEC.NO_AUX   = 0.259;
ELEC.OX_WE   = 0.229;  ELEC.OX_AUX   = 0.227;   % OX uses O3 columns
ELEC.NO2_WE  = 0.238;  ELEC.NO2_AUX  = 0.225;
% ---------- SENSOR ZERO BASELINES (clean air, include electronics if measured that way) ----------
ZERO.NO_WE   = 0.3158;  ZERO.NO_AUX   = 0.3017;   % NO-B4
ZERO.OX_WE   = 0.2366;  ZERO.OX_AUX   = 0.2278;   % OX-A431 (O3)
ZERO.NO2_WE  = 0.2396;  ZERO.NO2_AUX  = 0.2253;   % NO2-B43F
% ===========================================================
% --- Ensure expected columns exist (adjust here if names differ) ---
need = ["NO WE","NO AUX","O3 WE","O3 AUX","NO2 WE","NO2 AUX"];
have = string(T.Properties.VariableNames);
if ~all(ismember(need, have))
    error('Missing columns. Expected: %s\nHave: %s', strjoin(need,', '), strjoin(have,', '));
end
% --- Initialise structs to avoid dot-indexing errors ---
WE = struct(); AE = struct(); WEo = struct(); AEo = struct();
% --- Convert zeros to sensor-only baselines (remove electronics bias) ---
WEo.NO  = ZERO.NO_WE  - ELEC.NO_WE;   AEo.NO  = ZERO.NO_AUX  - ELEC.NO_AUX;
WEo.OX  = ZERO.OX_WE  - ELEC.OX_WE;   AEo.OX  = ZERO.OX_AUX  - ELEC.OX_AUX;
WEo.NO2 = ZERO.NO2_WE - ELEC.NO2_WE;  AEo.NO2 = ZERO.NO2_AUX - ELEC.NO2_AUX;
% --- Raw -> electronics-corrected live readings ---
WE.NO   = T.("NO WE")   - ELEC.NO_WE;    AE.NO   = T.("NO AUX")  - ELEC.NO_AUX;
WE.OX   = T.("O3 WE")   - ELEC.OX_WE;    AE.OX   = T.("O3 AUX")  - ELEC.OX_AUX;
WE.NO2  = T.("NO2 WE")  - ELEC.NO2_WE;   AE.NO2  = T.("NO2 AUX") - ELEC.NO2_AUX;
% ---------------- Temperature compensation tables ----------------
T_pts   = [-30 -20 -10  0  10  20  30  40  50];
% NO (kT, ratio model)
kT_NO   = interp1(T_pts, [1.8 1.8 1.4 1.1 1.1 1.0 0.9 0.9 0.8], TempC, 'linear','extrap');
% NO2 & OX (nT, simple model) – replace OX row with your own if you have it
nT_NO2  = interp1(T_pts, [1.3 1.3 1.3 1.3 1.0 0.6 0.4 0.2 -1.5], TempC, 'linear','extrap');
nT_OX   = interp1(T_pts, [1.3 1.3 1.3 1.3 1.0 0.6 0.4 0.2 -1.5], TempC, 'linear','extrap');
% ---------------- Apply algorithms ----------------
% NO (Algorithm 2 / ratio): WEc = (WEu - WEe) - kT*(WEo/AEo)*(AEu - AEe)
WEc_NO   = (WE.NO) - kT_NO .* (ZERO.NO_WE ./ZERO.NO_AUX ) .* (AE.NO );
% NO2 (Algorithm 1 / nT):   WEc = (WEu - WEe) - nT*(AEu - AEe)
%WEc_NO2  = (WE.NO2 - WEo.NO2) - nT_NO2 .* (AE.NO2 - AEo.NO2);
WEc_NO2  = (WE.NO2) - nT_NO2 .* (AE.NO2);
% OX (Algorithm 1 / nT):    WEc = (WEu - WEe) - nT*(AEu - AEe)
WEc_OX   = (WE.OX ) - nT_OX  .* (AE.OX );
% ---------------- Plot sets ----------------
makeFig("NO (k_T ratio)",  time, TempC, Hum, T.("NO WE") ,  T.("NO AUX"),  WEc_NO);
makeFig("NO2 (n_T simple)",time, TempC, Hum, T.("NO2 WE") , T.("NO2 AUX"), WEc_NO2);
makeFig("OX/O3 (n_T simple)",time,TempC, Hum, T.("O3 WE") ,  T.("O3 AUX"),  WEc_OX);
%% --------------- Plot helper ---------------
function makeFig(name,t,tc,hum,WE,AUX,WEc)
f = figure('Name',name,'Color','w');
tl = tiledlayout(f,2,1,'TileSpacing','compact','Padding','compact');
% Top: Temperature & Humidity
ax1 = nexttile(tl);
yyaxis left,  plot(t,tc,'LineWidth',1.2), ylabel('Temperature (°C)')
yyaxis right, plot(t,hum,'LineWidth',1.2), ylabel('Humidity (%)')
grid on, title([name,' — Temp & Humidity'])
% Bottom: WE/AUX vs Corrected
ax2 = nexttile(tl);
yyaxis left
plot(t,WE,'LineWidth',1.2); hold on
plot(t,AUX,'LineWidth',1.2); ylabel('Raw Sensor (V)')
yyaxis right
plot(t,WEc,'LineWidth',1.5); ylabel('Corrected WE_c (V)')
xlabel('Time'), grid on
title([name,' — WE/AUX (left), WE_c (right)'])
legend({'WE','AUX','WE Corrected*'},'Location','best')
linkaxes([ax1 ax2],'x');
end

[Star Strider]

My pleasure!

You provided the .m file for your code, however not the data ('SensorLog_ALL.mat'). (I checked -- the data are not included., at least that I could see.)

It would help to have it.

[Dharmesh]

Sorry its attached now

%% AAN-803-05 corrections: kT for NO, nT for NO2 & OX (O3)
% Dharmesh – edit the folder if needed:
%folder  = 'G:\AVR_Project\MATLAB Projects\Air Sensor\Temp_test';
matFile = 'SensorLog_ALL.mat';
% --- Load the combined table ---
S = load(matFile);
fn = fieldnames(S);
T  = S.(fn{find(structfun(@istable, S), 1)});
% --- Time & environment ---
if ~isdatetime(T.Timestamp)
    T.Timestamp = datetime(T.Timestamp,'InputFormat','yyyy-MM-dd HH:mm:ss');
end
time  = T.Timestamp;
TempC = T.("Temperature (°C)");
Hum   = T.("Humidity (%)");
% ================== CALIBRATION CONSTANTS ==================
% ---------- ELECTRONICS OFFSETS (hardware bias; set if measured) ----------
ELEC.NO_WE   = 0.292;  ELEC.NO_AUX   = 0.259;
ELEC.OX_WE   = 0.229;  ELEC.OX_AUX   = 0.227;   % OX uses O3 columns
ELEC.NO2_WE  = 0.238;  ELEC.NO2_AUX  = 0.225;
% ---------- SENSOR ZERO BASELINES (clean air, include electronics if measured that way) ----------
ZERO.NO_WE   = 0.3158;  ZERO.NO_AUX   = 0.3017;   % NO-B4
ZERO.OX_WE   = 0.2366;  ZERO.OX_AUX   = 0.2278;   % OX-A431 (O3)
ZERO.NO2_WE  = 0.2396;  ZERO.NO2_AUX  = 0.2253;   % NO2-B43F
% ===========================================================
% --- Ensure expected columns exist (adjust here if names differ) ---
need = ["NO WE","NO AUX","O3 WE","O3 AUX","NO2 WE","NO2 AUX"];
have = string(T.Properties.VariableNames);
if ~all(ismember(need, have))
    error('Missing columns. Expected: %s\nHave: %s', strjoin(need,', '), strjoin(have,', '));
end
% --- Initialise structs to avoid dot-indexing errors ---
WE = struct(); AE = struct(); WEo = struct(); AEo = struct();
% --- Convert zeros to sensor-only baselines (remove electronics bias) ---
WEo.NO  = ZERO.NO_WE  - ELEC.NO_WE;   AEo.NO  = ZERO.NO_AUX  - ELEC.NO_AUX;
WEo.OX  = ZERO.OX_WE  - ELEC.OX_WE;   AEo.OX  = ZERO.OX_AUX  - ELEC.OX_AUX;
WEo.NO2 = ZERO.NO2_WE - ELEC.NO2_WE;  AEo.NO2 = ZERO.NO2_AUX - ELEC.NO2_AUX;
% --- Raw -> electronics-corrected live readings ---
WE.NO   = T.("NO WE")   - ELEC.NO_WE;    AE.NO   = T.("NO AUX")  - ELEC.NO_AUX;
WE.OX   = T.("O3 WE")   - ELEC.OX_WE;    AE.OX   = T.("O3 AUX")  - ELEC.OX_AUX;
WE.NO2  = T.("NO2 WE")  - ELEC.NO2_WE;   AE.NO2  = T.("NO2 AUX") - ELEC.NO2_AUX;
% ---------------- Temperature compensation tables ----------------
T_pts   = [-30 -20 -10  0  10  20  30  40  50];
% NO (kT, ratio model)
kT_NO   = interp1(T_pts, [1.8 1.8 1.4 1.1 1.1 1.0 0.9 0.9 0.8], TempC, 'linear','extrap');
% NO2 & OX (nT, simple model) – replace OX row with your own if you have it
nT_NO2  = interp1(T_pts, [1.3 1.3 1.3 1.3 1.0 0.6 0.4 0.2 -1.5], TempC, 'linear','extrap');
nT_OX   = interp1(T_pts, [1.3 1.3 1.3 1.3 1.0 0.6 0.4 0.2 -1.5], TempC, 'linear','extrap');
% ---------------- Apply algorithms ----------------
% NO (Algorithm 2 / ratio): WEc = (WEu - WEe) - kT*(WEo/AEo)*(AEu - AEe)
WEc_NO   = (WE.NO) - kT_NO .* (ZERO.NO_WE ./ZERO.NO_AUX ) .* (AE.NO );
% NO2 (Algorithm 1 / nT):   WEc = (WEu - WEe) - nT*(AEu - AEe)
%WEc_NO2  = (WE.NO2 - WEo.NO2) - nT_NO2 .* (AE.NO2 - AEo.NO2);
WEc_NO2  = (WE.NO2) - nT_NO2 .* (AE.NO2);
% OX (Algorithm 1 / nT):    WEc = (WEu - WEe) - nT*(AEu - AEe)
WEc_OX   = (WE.OX ) - nT_OX  .* (AE.OX );
% ---------------- Plot sets ----------------
makeFig("NO (k_T ratio)",  time, TempC, Hum, T.("NO WE") ,  T.("NO AUX"),  WEc_NO);

makeFig("NO2 (n_T simple)",time, TempC, Hum, T.("NO2 WE") , T.("NO2 AUX"), WEc_NO2);

●

makeFig("OX/O3 (n_T simple)",time,TempC, Hum, T.("O3 WE") ,  T.("O3 AUX"),  WEc_OX);

%% --------------- Plot helper ---------------
function makeFig(name,t,tc,hum,WE,AUX,WEc)
f = figure('Name',name,'Color','w');
tl = tiledlayout(f,2,1,'TileSpacing','compact','Padding','compact');
% Top: Temperature & Humidity
ax1 = nexttile(tl);
yyaxis left,  plot(t,tc,'LineWidth',1.2), ylabel('Temperature (°C)')
yyaxis right, plot(t,hum,'LineWidth',1.2), ylabel('Humidity (%)')
grid on, title([name,' — Temp & Humidity'])
% Bottom: WE/AUX vs Corrected
ax2 = nexttile(tl);
yyaxis left
plot(t,WE,'LineWidth',1.2); hold on
plot(t,AUX,'LineWidth',1.2); ylabel('Raw Sensor (V)')
yyaxis right
plot(t,WEc,'LineWidth',1.5); ylabel('Corrected WE_c (V)')
xlabel('Time'), grid on
title([name,' — WE/AUX (left), WE_c (right)'])
legend({'WE','AUX','WE Corrected*'},'Location','best')
linkaxes([ax1 ax2],'x');
end

[Star Strider]

Torsten -- Thank you!

@Dharmesh -- I am having a bit of trouble understnading this.

What are the various variables? I do not completely understand 'WE" and 'AUX'.

Do you have the actual NOx values? (Is this a calibration test?)

I do not understand the technology of the sensor, and this may be necessary to devise a way of correcting its readings. Do you have a published way of correcting the sensor values for the humidity? Is there a specific way that humidity interferes with the sensor (does it 'poison' it, physically block it, or something else)?

EDIT -- (24 Oct 2025 at 00:45)

See if the sensors dexcribed in HUMIDITY AND TEMPERATURE CORRECTION FACTORS FOR NOX EMISSIONS FROM DIESEL ENGINES matches your sensor. If so, we can probably use one of these correction approaches. Those appear to be relatively straightforward, however it will be necessary to state the units of your data, and make any necessary unit conversions to work with these equations. Specifically, Equation (4) would be easy to code.

.

[Dharmesh]

@Star Strider

Sorry. , the sensor outputs two signals: WE and AUX. The difference between them is that WE includes the gas concentration, while AUX does not. The purpose of this is so that AUX can be used to identify whether any noise or environmental factors have affected the sensor.

However, both signals are influenced by temperature and humidity. Temperature is relatively simple to correct you ensure that the gas concentration is zero (or very close to zero) and then create a baseline by measuring at various temperature steps.

In this case, we can assume that the gas concentration is close to zero. At this stage, we do not need to calibrate the concentration.

The signal we need to focus on is the corrected WE signal, shown in red. Primarily, you can see humidity transients, which simply affect the current and therefore the final voltage.

There is some information on these sensors in the application note, please see the last page:

aan_110_app-note_v0.pdf

Sensor data sheet for NO2. All other sensors would in principle be similar

alphasense_no2-b43f_datasheet_en_4.pdf

[Star Strider]

I am still having problems understanding this.

I thought this would be more straightforward than it turned out to be. (I have a chemistry background, however this is far ffrom my areas of expertise.)

I looked up and implemented Equation (3) from the paper I cited in my previous Comment. (There actually is not nuch information on this, at least that I can find.) I wrote a simple function for the formula that appeared to be most applicable, however the absolute humidity units (this formula gives g/m^3, while the equations in the paper require g/kg, so that may require an additional calculation) I then did an example calculation using observed data from the first two subplots for the temperatire, relative humidity, and NO concentration to calculate a 'corrected' NO value. Those are at the end of the code block.

This is the best that I can do. I provided a sample calculation and plots of although I did not apply the correction to the data.

%% AAN-803-05 corrections: kT for NO, nT for NO2 & OX (O3)
% Dharmesh – edit the folder if needed:
%folder  = 'G:\AVR_Project\MATLAB Projects\Air Sensor\Temp_test';
matFile = 'SensorLog_ALL.mat';
% --- Load the combined table ---
S = load(matFile);
fn = fieldnames(S);
T  = S.(fn{find(structfun(@istable, S), 1)});
% --- Time & environment ---
if ~isdatetime(T.Timestamp)
    T.Timestamp = datetime(T.Timestamp,'InputFormat','yyyy-MM-dd HH:mm:ss');
end
time  = T.Timestamp;
TempC = T.("Temperature (°C)");
Hum   = T.("Humidity (%)");
% ================== CALIBRATION CONSTANTS ==================
% ---------- ELECTRONICS OFFSETS (hardware bias; set if measured) ----------
ELEC.NO_WE   = 0.292;  ELEC.NO_AUX   = 0.259;
ELEC.OX_WE   = 0.229;  ELEC.OX_AUX   = 0.227;   % OX uses O3 columns
ELEC.NO2_WE  = 0.238;  ELEC.NO2_AUX  = 0.225;
% ---------- SENSOR ZERO BASELINES (clean air, include electronics if measured that way) ----------
ZERO.NO_WE   = 0.3158;  ZERO.NO_AUX   = 0.3017;   % NO-B4
ZERO.OX_WE   = 0.2366;  ZERO.OX_AUX   = 0.2278;   % OX-A431 (O3)
ZERO.NO2_WE  = 0.2396;  ZERO.NO2_AUX  = 0.2253;   % NO2-B43F
% ===========================================================
% --- Ensure expected columns exist (adjust here if names differ) ---
need = ["NO WE","NO AUX","O3 WE","O3 AUX","NO2 WE","NO2 AUX"];
have = string(T.Properties.VariableNames);
if ~all(ismember(need, have))
    error('Missing columns. Expected: %s\nHave: %s', strjoin(need,', '), strjoin(have,', '));
end
% --- Initialise structs to avoid dot-indexing errors ---
WE = struct(); AE = struct(); WEo = struct(); AEo = struct();
% --- Convert zeros to sensor-only baselines (remove electronics bias) ---
WEo.NO  = ZERO.NO_WE  - ELEC.NO_WE;   AEo.NO  = ZERO.NO_AUX  - ELEC.NO_AUX;
WEo.OX  = ZERO.OX_WE  - ELEC.OX_WE;   AEo.OX  = ZERO.OX_AUX  - ELEC.OX_AUX;
WEo.NO2 = ZERO.NO2_WE - ELEC.NO2_WE;  AEo.NO2 = ZERO.NO2_AUX - ELEC.NO2_AUX;
% --- Raw -> electronics-corrected live readings ---
WE.NO   = T.("NO WE")   - ELEC.NO_WE;    AE.NO   = T.("NO AUX")  - ELEC.NO_AUX;
WE.OX   = T.("O3 WE")   - ELEC.OX_WE;    AE.OX   = T.("O3 AUX")  - ELEC.OX_AUX;
WE.NO2  = T.("NO2 WE")  - ELEC.NO2_WE;   AE.NO2  = T.("NO2 AUX") - ELEC.NO2_AUX;
% ---------------- Temperature compensation tables ----------------
T_pts   = [-30 -20 -10  0  10  20  30  40  50];
% NO (kT, ratio model)
kT_NO   = interp1(T_pts, [1.8 1.8 1.4 1.1 1.1 1.0 0.9 0.9 0.8], TempC, 'linear','extrap');
% NO2 & OX (nT, simple model) – replace OX row with your own if you have it
nT_NO2  = interp1(T_pts, [1.3 1.3 1.3 1.3 1.0 0.6 0.4 0.2 -1.5], TempC, 'linear','extrap');
nT_OX   = interp1(T_pts, [1.3 1.3 1.3 1.3 1.0 0.6 0.4 0.2 -1.5], TempC, 'linear','extrap');
% ---------------- Apply algorithms ----------------
% NO (Algorithm 2 / ratio): WEc = (WEu - WEe) - kT*(WEo/AEo)*(AEu - AEe)
WEc_NO   = (WE.NO) - kT_NO .* (ZERO.NO_WE ./ZERO.NO_AUX ) .* (AE.NO );
% NO2 (Algorithm 1 / nT):   WEc = (WEu - WEe) - nT*(AEu - AEe)
%WEc_NO2  = (WE.NO2 - WEo.NO2) - nT_NO2 .* (AE.NO2 - AEo.NO2);
WEc_NO2  = (WE.NO2) - nT_NO2 .* (AE.NO2);
% OX (Algorithm 1 / nT):    WEc = (WEu - WEe) - nT*(AEu - AEe)
WEc_OX   = (WE.OX ) - nT_OX  .* (AE.OX );
% ---------------- Plot sets ----------------
makeFig("NO (k_T ratio)",  time, TempC, Hum, T.("NO WE") ,  T.("NO AUX"),  WEc_NO);

makeFig("NO2 (n_T simple)",time, TempC, Hum, T.("NO2 WE") , T.("NO2 AUX"), WEc_NO2);

makeFig("OX/O3 (n_T simple)",time,TempC, Hum, T.("O3 WE") ,  T.("O3 AUX"),  WEc_OX);

●

%% --------------- Plot helper ---------------
function makeFig(name,t,tc,hum,WE,AUX,WEc)
f = figure('Name',name,'Color','w');
tl = tiledlayout(f,2,1,'TileSpacing','compact','Padding','compact');
% Top: Temperature & Humidity
ax1 = nexttile(tl);
yyaxis left,  plot(t,tc,'LineWidth',1.2), ylabel('Temperature (°C)')
yyaxis right, plot(t,hum,'LineWidth',1.2), ylabel('Humidity (%)')
grid on, title([name,' — Temp & Humidity'])
% Bottom: WE/AUX vs Corrected
ax2 = nexttile(tl);
yyaxis left
plot(t,WE,'LineWidth',1.2); hold on
plot(t,AUX,'LineWidth',1.2); ylabel('Raw Sensor (V)')
yyaxis right
plot(t,WEc,'LineWidth',1.5); ylabel('Corrected WE_c (V)')
xlabel('Time'), grid on
title([name,' — WE/AUX (left), WE_c (right)'])
legend({'WE','AUX','WE Corrected*'},'Location','best')
linkaxes([ax1 ax2],'x');
end


KNOx = Fritz(40,80)                     % Calculated Correction Factor

KNOx = 0.5042

NOcorr = [0.45; 0.35]*KNOx              % Corrected Values

NOcorr = 2×1

0.2269 0.1765

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

KNOxv = Fritz(TempC, Hum);              % Correctio9ns Vector Correspoinding TO Top Plot In Each Figure

figure
plot(time, KNOxv)
grid
xlabel('Time')
ylabel('KNO_x')
title('KNO_x Correction Factor')

figure
plot3(TempC, Hum, KNOxv)
grid
xlabel('T (ºC)')
ylabel('Relative Humidity (%)')
zlabel('KNO_x')
title('KNO_x Correction Factor')

function KNOx = Fritz(T,RH)
% Source KNOx: https://wiki.unece.org/download/attachments/51972522/EPPR-21-17%20Humidity%20and%20Temperature%20correction%20factors%20for%20NOx%20from%20diesel.pdf?api=v2
% T = Ambient Temperature ºC
% RH = Relative Humidity %

% Source Habs: https://www.reddit.com/r/embedded/comments/f4iino/whats_the_right_way_to_convert_relative_humidity/
Habs = 216.7*(RH/100 .* 6.112 .* exp((17.62.*T)./(243.12+T)))./(273.15+T);  %Units: g/m^3

KNOx = 1+(0.00446*(T-25) - 0.018708*(Habs-10.71)); % Source KNOx Eqn (3)

end

.

[Dharmesh]

Thanks, I’ll add this to my code, but I don’t think it’s the relative humidity — it’s more the transients caused by the rate of humidity rise and fall that create the issue.

[Star Strider]

My pleasure!

This is the only correction that I can find that makes any sense. I noticed in the application notes that humidity changes cause transients with exponential decaay characteristics back to baseline. I have no way of modeling that, because I do not understand the sensor technoloby thoroughly enough. It would be necessary to do that experiment with your sensor and then create the appropriate regression using the input-output characteristics. Apparently only the transients are significant, and the sensor seems to correct itself to the ambient humidity otherwise.

If this is not what you want, then I will delete my answer.

[Dharmesh]

@Star Strider

Yes, what you’ve observed is correct as humidity changes, there is an exponential delay in returning to the baseline. Humidity transients are more noticeable in NO₂ and OX sensors, while NO sensors do not appear to be affected as much. I will ask the manufacturer why this occurs.

Here’s a small idea please advise if it’s feasible.

From simple observation, I can estimate what the baseline signal should be, assuming the humidity remains stable. At this stage, it might not be 100% accurate, but it could, in principle, generate a signal representing temperatures from 10°C to 45°C. This could be implemented as a simple lookup table.

We could then use this generated signal as a reference for what we expect the signal to be, and compare it with the actual WE signal using a model ,perhaps with the Regression Learner app.

If this approach works, we need to consider how to include humidity in the model. I don’t think relative humidity (RH) alone would be sufficient; instead, we may need to look at the rate of change (± humidity variation) over a certain duration. We can test different time windows for example, 1 min, 5 min, and 10 min . Some times humdity could be increasing slowly throught out the day.

Another important consideration is that the model should not rely on an absolute reference to the baseline. If a sensor has a different baseline, the correction should still work. Therefore, the correction factor should be expressed as a ± percentage relative to its baseline.

[Star Strider]

The 'lookup table' coould be implemented with one of the interpolation functions in core MATLAB. That part would not be difficult.

I have also been thinking about this with respect to the humidity transients, since that seems to be the most important problem.

Looking at the 'Humidity Transients' plot in 'App_note_v0', one option would be to re-create that experiment with numidity 'steps' and record the output. I am not certain what the sampling frequency would be, however considering that the note mentions that the numidity transients 'decay in about 10 minutes', sampling at 1 Hz might be enough, although more data is always better, so use a higher sampling frequency if system memory permits. Then use the System Identification Toolbox to derive the system structure (most likely a state space representation) choosing the appropriate order and using the compare function to determine the best system order. With that information, you can use the recorded humidity as an input, calculate the resulting output (simulate the state space representation with the recorded humidity signal), and then subtract that output signal from the observed result. That is the only way I can think of to correct for the transients. I doubt that any other approach would work.

The 'Humidity Transients' plot does not report either the input signal (probably some sort of step-wise humidity change) or the time, so using that plot as it exists would not be appropriate. (Without knowing the precise input signal and its shape, and the associated times, that plot is quantitatively useless.) It would be necessary to duplicate that experiment and mathematically identify the system in order to model the sensor response to humidity.

I can help with the system identification process once you have those data.

.

[Dharmesh]

@Star Strider I’ve collected new and improved data, please see below. I used an environmental chamber, which allowed me to stabilise the humidity much better than before.

Before I present my code and data for us to work on the humidity model, I’d like to clean up the data a little so we can present a reference signal that represents the expected behaviour.

The areas that need cleaning are as follows:

1. Start and End: These sections correspond to when I started and stopped the test, during which the machine was powering up or shutting down. I can simply remove these timestamps from the data file.
2. Temperature Steps: Each time the temperature increases, it causes a step change, and the humidity shifts slightly, creating small transients in the signal. After each temperature step, I wait for 10 minutes before moving to the next step. Therefore, I’d like to discard the step period and the first 5 minutes following each step.

Is there a function that can do this automatically instead of manually removing the data?

[Star Strider]

I would not eliminate or edit anything. The transients in the step changes are important. If you know the inpuits (specifically the temperature changes and whatever else was the input signal, such as humidity changes), and have that recorrd, that is all that would be necessary. The input signal needs to be provided as part of the record. I would include the start and end, unless the instrumentation was just coming online and was not stable in those regions.

I am not certain how the humidity changes with temperature, although that might be automatically included in the estimated system. It is linear, (or at least linear in the parameters) in the region-of-interest, that would be enough.

I am still not certain what you are doing, however the data appear to be good.

The parameter estimation will model the system put to it in the data. If you are simply interested in humidity transients, you need to increase, decresas, (or both) the humidity, record those changes in the input signal, and record the system output. Temperature variations will model temperature. Everything else is inference.

[Dharmesh]

@Star Strider I’ve uploaded my new data. The irregular readings at the start occurred while the chamber was stabilising itself, and those at the end appeared when I set the machine to reduce its temperature before opening the door.

To explain what I’m trying to achieve:

The sensor’s chemistry is affected by both temperature and humidity. When measuring gas concentration, it’s necessary to subtract the sensor’s baseline, which itself varies with temperature. My goal is to keep the humidity level as stable as possible (though this might be somewhat unrealistic) and vary the temperature from 0 °C to 50 °C. Assuming there is no gas present in the environment, this should produce a temperature-dependent baseline for the sensor.

Across the 0 °C–50 °C temperature range, the effect of humidity on the baseline appears relatively minor, except during temperature step changes. This happens because, as temperature rises, air can hold more moisture, which increases relative humidity. However, the chamber compensates for this by activating its fans, which gradually restore the environment to a stable condition, and the sensor readings also stabilise.

Therefore, when constructing the baseline, I’m mainly interested in using the data recorded when the sensor is in its most stable state. So the reason for asking, how we can remove data at certain points in time.

Once this baseline is established,which represents the sensor’s output in the absence of gas and becomes our reference signal for any temperature window,we can then study the effects of humidity. This can be done using the original data or by setting the environmental chamber to introduce various humidity transients (both rises and falls) to observe how the sensor responds and how long it takes to stabilise after each transient.

So Step 1: Keep Humdity Static, and produces a function for Temperture Change.

Step 2: Keep Temperture Static, and produces a function for Humdity Change.

[Star Strider]

I am still having significant problems with your data, in part because I still do not ocmpletely understand what you are doing.

I calculated the absolute hiumidity from the tempearature and relative humidity, according to one of the functions I discovered for that, then attempted a system identification using that and one input. (Apparently, absolute humidity is necessary for the calculations you want to do.) They are delineated in my modified version of your code as such at the end of your original code. (I ran this in MATLAB Online, that being the reson for the websave call.)

Signal processing techniques, including syystem identification, require regularly-sampled data (and assume that the data is refularly-sampled), so I calculated the mean sampling interval, copied and converted your 'T' table to a timetable, and then resampled it using the retime funciton. (These are all core MATLAB functions.) I then used System Identification Toolbox functions iddata and ssest (estimates the system as a state space realisation) to see if I could reproduce the data. The fit of a order system (using compare to visualise it) is moderately good, although it obviously needs more inputs or a different output. (I do not understand your system well enough to decide what the inputs and outputs should be.)

The important part of this exercise is that it sets your data up for further experiments with system identification. The System Identification Toolbox can estimate multi-input, multi-output (MIMO) systems as well as the single-input, single-output (SISO) system I modeled here.

In short, we need to know what the inputs and outputs should be in order to model this correctly. You understand what those are. At present, I do not.

Try this --

close all

matFile = websave('SensorLog_ALLS1.mat','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1842649/SensorLog_ALLS1.mat');

%% AAN-803-05 corrections: kT for NO, nT for NO2 & OX (O3)
% Dharmesh – edit the folder if needed:
%folder  = 'G:\AVR_Project\MATLAB Projects\Air Sensor\Temp_test';
% matFile = 'SensorLog_ALL.mat';
% --- Load the combined table ---
S = load(matFile);
fn = fieldnames(S);
T  = S.(fn{find(structfun(@istable, S), 1)});
% --- Time & environment ---
if ~isdatetime(T.Timestamp)
    T.Timestamp = datetime(T.Timestamp,'InputFormat','yyyy-MM-dd HH:mm:ss');
end
time  = T.Timestamp;
TempC = T.("Temperature (°C)");
Hum   = T.("Humidity (%)");
% ================== CALIBRATION CONSTANTS ==================
% ---------- ELECTRONICS OFFSETS (hardware bias; set if measured) ----------
ELEC.NO_WE   = 0.292;  ELEC.NO_AUX   = 0.259;
ELEC.OX_WE   = 0.229;  ELEC.OX_AUX   = 0.227;   % OX uses O3 columns
ELEC.NO2_WE  = 0.238;  ELEC.NO2_AUX  = 0.225;
% ---------- SENSOR ZERO BASELINES (clean air, include electronics if measured that way) ----------
ZERO.NO_WE   = 0.3158;  ZERO.NO_AUX   = 0.3017;   % NO-B4
ZERO.OX_WE   = 0.2366;  ZERO.OX_AUX   = 0.2278;   % OX-A431 (O3)
ZERO.NO2_WE  = 0.2396;  ZERO.NO2_AUX  = 0.2253;   % NO2-B43F
% ===========================================================
% --- Ensure expected columns exist (adjust here if names differ) ---
need = ["NO WE","NO AUX","O3 WE","O3 AUX","NO2 WE","NO2 AUX"];
have = string(T.Properties.VariableNames);
if ~all(ismember(need, have))
    error('Missing columns. Expected: %s\nHave: %s', strjoin(need,', '), strjoin(have,', '));
end
% --- Initialise structs to avoid dot-indexing errors ---
WE = struct(); AE = struct(); WEo = struct(); AEo = struct();
% --- Convert zeros to sensor-only baselines (remove electronics bias) ---
WEo.NO  = ZERO.NO_WE  - ELEC.NO_WE;   AEo.NO  = ZERO.NO_AUX  - ELEC.NO_AUX;
WEo.OX  = ZERO.OX_WE  - ELEC.OX_WE;   AEo.OX  = ZERO.OX_AUX  - ELEC.OX_AUX;
WEo.NO2 = ZERO.NO2_WE - ELEC.NO2_WE;  AEo.NO2 = ZERO.NO2_AUX - ELEC.NO2_AUX;
% --- Raw -> electronics-corrected live readings ---
WE.NO   = T.("NO WE")   - ELEC.NO_WE;    AE.NO   = T.("NO AUX")  - ELEC.NO_AUX;
WE.OX   = T.("O3 WE")   - ELEC.OX_WE;    AE.OX   = T.("O3 AUX")  - ELEC.OX_AUX;
WE.NO2  = T.("NO2 WE")  - ELEC.NO2_WE;   AE.NO2  = T.("NO2 AUX") - ELEC.NO2_AUX;
% ---------------- Temperature compensation tables ----------------
T_pts   = [-30 -20 -10  0  10  20  30  40  50];
% NO (kT, ratio model)
kT_NO   = interp1(T_pts, [1.8 1.8 1.4 1.1 1.1 1.0 0.9 0.9 0.8], TempC, 'linear','extrap');
% NO2 & OX (nT, simple model) – replace OX row with your own if you have it
nT_NO2  = interp1(T_pts, [1.3 1.3 1.3 1.3 1.0 0.6 0.4 0.2 -1.5], TempC, 'linear','extrap');
nT_OX   = interp1(T_pts, [1.3 1.3 1.3 1.3 1.0 0.6 0.4 0.2 -1.5], TempC, 'linear','extrap');
% ---------------- Apply algorithms ----------------
% NO (Algorithm 2 / ratio): WEc = (WEu - WEe) - kT*(WEo/AEo)*(AEu - AEe)
WEc_NO   = (WE.NO) - kT_NO .* (ZERO.NO_WE ./ZERO.NO_AUX ) .* (AE.NO );
% NO2 (Algorithm 1 / nT):   WEc = (WEu - WEe) - nT*(AEu - AEe)
%WEc_NO2  = (WE.NO2 - WEo.NO2) - nT_NO2 .* (AE.NO2 - AEo.NO2);
WEc_NO2  = (WE.NO2) - nT_NO2 .* (AE.NO2);
% OX (Algorithm 1 / nT):    WEc = (WEu - WEe) - nT*(AEu - AEe)
WEc_OX   = (WE.OX ) - nT_OX  .* (AE.OX );
% ---------------- Plot sets ----------------
makeFig("NO (k_T ratio)",  time, TempC, Hum, T.("NO WE") ,  T.("NO AUX"),  WEc_NO);

makeFig("NO2 (n_T simple)",time, TempC, Hum, T.("NO2 WE") , T.("NO2 AUX"), WEc_NO2);

makeFig("OX/O3 (n_T simple)",time,TempC, Hum, T.("O3 WE") ,  T.("O3 AUX"),  WEc_OX);

% ===== System Identification Start =====

Ts = mean(diff(T.Timestamp));
Ts.Format = 'hh:mm:ss.SSS'

Ts = duration
   00:00:05.513

Tsd = std(diff(T.Timestamp));
Tsd.Format = 'hh:mm:ss.SSS'

Tsd = duration
   00:00:01.229

Tss = seconds(Ts)

Tss = 5.5140

Fs = 1/Tss                                          % Sampling Frequency (Hz)

Fs = 0.1814

% Source Habs: https://www.reddit.com/r/embedded/comments/f4iino/whats_the_right_way_to_convert_relative_humidity/
Habs = @(RH,T) 216.7*(RH/100 .* 6.112 .* exp((17.62.*T)./(243.12+T)))./(273.15+T);  %Units: g/m^3

HumAbs = Habs(T.('Humidity (%)'), T.('Temperature (°C)'));

Tr1 = T;
Tr1 = addvars(Tr1, HumAbs, 'Before','NO WE', NewVariableNames='Absolute Humidity');
TTr1 = table2timetable(Tr1(:,1:end-1));
fprintf('\nResampled Data --\n')

Resampled Data --

TTr2 = retime(TTr1, 'regular', 'linear', 'TimeStep',seconds(Tss))

TTr2 = 2683×9 timetable
         Timestamp         Temperature (°C)    Humidity (%)    Absolute Humidity     NO WE     NO AUX      O3 WE     O3 AUX     NO2 WE     NO2 AUX
    ___________________    ________________    ____________    _________________    _______    _______    _______    _______    _______    _______

    2025-10-30 11:04:49         9.6125            55.839             5.112          0.25413      0.248      0.238    0.23587    0.26687    0.23513
    2025-10-30 11:04:55         9.5463            55.651            5.0736          0.25523      0.248      0.238    0.23477    0.26577    0.23623
    2025-10-30 11:05:00         9.4835             55.46            5.0361          0.25567      0.248      0.238    0.23433    0.26533      0.237
    2025-10-30 11:05:06          9.411            55.197            4.9891          0.25543      0.248      0.238      0.235      0.266      0.237
    2025-10-30 11:05:11         9.3385             54.93            4.9421            0.256      0.248    0.23854      0.235    0.26546      0.237
    2025-10-30 11:05:17         9.2987            54.785            4.9166            0.256      0.248      0.239      0.235      0.265      0.237
    2025-10-30 11:05:22         9.2583            54.655            4.8923          0.25708    0.24854      0.239      0.235    0.26554    0.23754
    2025-10-30 11:05:28         9.2142            54.516            4.8662            0.258      0.249      0.239      0.235      0.266      0.238
    2025-10-30 11:05:33         9.1776            54.373            4.8421          0.25875    0.24975    0.23975      0.235    0.26525      0.238
    2025-10-30 11:05:39         9.1445            54.229            4.8191            0.259       0.25       0.24    0.23585      0.265    0.23715
    2025-10-30 11:05:44         9.1019            54.086            4.7933          0.25995       0.25       0.24      0.236      0.265      0.237
    2025-10-30 11:05:50         9.0683            53.957            4.7717          0.26006       0.25    0.24006    0.23506      0.265    0.23706
    2025-10-30 11:05:55          9.036            53.913            4.7579          0.26113    0.25013    0.24113      0.236      0.265      0.238
    2025-10-30 11:06:01         9.0075            53.926            4.7504            0.262    0.25106    0.24194      0.236    0.26494      0.238
    2025-10-30 11:06:06         8.9634            53.819            4.7277            0.262      0.252    0.24067      0.236    0.26384      0.238
    2025-10-30 11:06:12          8.914            53.484            4.6835          0.26173    0.25173     0.2382      0.236    0.26273      0.238

% return

figure(Name='Absolute Humidity (Calculated)')
plot(T.Timestamp, HumAbs)
grid
xlabel('Time')
ylabel('Absolute Humidity')

SysOrd = 4;

y = TTr2.('NO WE');
u = TTr2.('Absolute Humidity');
Ts = Tss

Ts = 5.5140

HumAbsData = iddata(y,u,Ts)

HumAbsData = 
Time domain data set with 2683 samples.
Sample time: 5.51399 seconds           
                                       
Outputs      Unit (if specified)       
   y1                                  
                                       
Inputs       Unit (if specified)       
   u1                                  
                                       

SysHumAbs = ssest(HumAbsData,SysOrd)

SysHumAbs =
  Continuous-time identified state-space model:
      dx/dt = A x(t) + B u(t) + K e(t)
       y(t) = C x(t) + D u(t) + e(t)
 
  A = 
               x1          x2          x3          x4
   x1  -0.0007405  -0.0009573   0.0009337   -0.002375
   x2    0.006086     0.01964    -0.01224   -0.002013
   x3    0.001077      0.0654     -0.0151     -0.2565
   x4   -0.003161    -0.04787     0.03169    -0.06633
 
  B = 
               u1
   x1   3.622e-05
   x2  -0.0004882
   x3  -0.0009286
   x4    0.001012
 
  C = 
              x1         x2         x3         x4
   y1      9.039     0.0142  -0.001957   0.001669
 
  D = 
       u1
   y1   0
 
  K = 
            y1
   x1   0.0233
   x2    1.111
   x3    1.068
   x4  0.08223
 
Parameterization:
   FREE form (all coefficients in A, B, C free).
   Feedthrough: none
   Disturbance component: estimate
   Number of free coefficients: 28
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                
Estimated using SSEST on time domain data "HumAbsData".
Fit to estimation data: 99.39% (prediction focus)      
FPE: 3.285e-07, MSE: 3.246e-07                         
 

figure
compare(HumAbsData, SysHumAbs)

% ===== System Identification End =====

%% --------------- Plot helper ---------------
function makeFig(name,t,tc,hum,WE,AUX,WEc)
f = figure('Name',name,'Color','w');
tl = tiledlayout(f,2,1,'TileSpacing','compact','Padding','compact');
% Top: Temperature & Humidity
ax1 = nexttile(tl);
yyaxis left,  plot(t,tc,'LineWidth',1.2), ylabel('Temperature (°C)')
yyaxis right, plot(t,hum,'LineWidth',1.2), ylabel('Humidity (%)')
grid on, title([name,' — Temp & Humidity'])
% Bottom: WE/AUX vs Corrected
ax2 = nexttile(tl);
yyaxis left
plot(t,WE,'LineWidth',1.2); hold on
plot(t,AUX,'LineWidth',1.2); ylabel('Raw Sensor (V)')
yyaxis right
plot(t,WEc,'LineWidth',1.5); ylabel('Corrected WE_c (V)')
xlabel('Time'), grid on
title([name,' — WE/AUX (left), WE_c (right)'])
legend({'WE','AUX','WE Corrected*'},'Location','best')
linkaxes([ax1 ax2],'x');
end

% ========================================================================

.

[Dharmesh]

@Star Strider

What I am trying to achieve is the following:

1. Correct the signal for the effect of temperature.
2. Correct the signal for the effect of humidity changes/transients.

Although both factors are related, they need to be addressed individually. Therefore, I tried to keep the humidity stable while varying the temperature so that I can calculate the coefficients for each temperature and create a corrected baseline.

In theory, this would be a two-step process according to the manufacturer’s documentation for the sensor modules.

The code I have so far performs the temperature correction (step 1). The next step is to generate data with varying humidity to better understand its effect. Please see the scatter plots — these produce the baseline I would expect as the temperature changes.

The calculated trend line is now my theoretical reference signal. The next step will be to collect data with varying humidity — both rising and falling — to observe its effect and determine how long the sensor takes to stabilise. I believe this data will help us accurately model the impact of humidity.

For this data collection/test, we will keep the temperature static at 20 °C while varying the humidity. For example, the sensor at 20 °C may produce a value of 230 mV. We need to gain a better understanding of how this value changes with humidity.

I will be conducting some additional tests later this week, after which I will share the results. Based on those, I believe we will need to develop a model to describe this behaviour.

matFile = fullfile(folder, 'SensorLog_ALLS2.mat');

% --- Load the combined table ---
S = load(matFile);
fn = fieldnames(S);
T  = S.(fn{find(structfun(@istable, S), 1)});

% --- Time & environment ---
if ~isdatetime(T.Timestamp)
    T.Timestamp = datetime(T.Timestamp,'InputFormat','yyyy-MM-dd HH:mm:ss');
end
time  = T.Timestamp;
TempC = T.("Temperature (°C)");
Hum   = T.("Humidity (%)");

% ================== CALIBRATION CONSTANTS ==================
% ---------- ELECTRONICS OFFSETS (hardware bias; set if measured) ----------
ELEC.NO_WE   = 0.292;  ELEC.NO_AUX   = 0.259;
ELEC.OX_WE   = 0.229;  ELEC.OX_AUX   = 0.227;   % OX uses O3 columns
ELEC.NO2_WE  = 0.238;  ELEC.NO2_AUX  = 0.225;

% ---------- SENSOR ZERO BASELINES (clean air, include electronics if measured that way) ----------
ZERO.NO_WE   = 0.3158;  ZERO.NO_AUX   = 0.3017;   % NO-B4
ZERO.OX_WE   = 0.2366;  ZERO.OX_AUX   = 0.2278;   % OX-A431 (O3)
ZERO.NO2_WE  = 0.2396;  ZERO.NO2_AUX  = 0.2253;   % NO2-B43F
% ===========================================================

% --- Ensure expected columns exist (adjust here if names differ) ---
need = ["NO WE","NO AUX","O3 WE","O3 AUX","NO2 WE","NO2 AUX"];
have = string(T.Properties.VariableNames);
if ~all(ismember(need, have))
    error('Missing columns. Expected: %s\nHave: %s', strjoin(need,', '), strjoin(have,', '));
end

% --- Initialise structs to avoid dot-indexing errors ---
WE = struct(); AE = struct(); WEo = struct(); AEo = struct();

% --- Convert zeros to sensor-only baselines (remove electronics bias) ---
WEo.NO  = ZERO.NO_WE  - ELEC.NO_WE;   AEo.NO  = ZERO.NO_AUX  - ELEC.NO_AUX;
WEo.OX  = ZERO.OX_WE  - ELEC.OX_WE;   AEo.OX  = ZERO.OX_AUX  - ELEC.OX_AUX;
WEo.NO2 = ZERO.NO2_WE - ELEC.NO2_WE;  AEo.NO2 = ZERO.NO2_AUX - ELEC.NO2_AUX;

% --- Raw -> electronics-corrected live readings ---
WE.NO   = T.("NO WE")   - ELEC.NO_WE;    AE.NO   = T.("NO AUX")  - ELEC.NO_AUX;
WE.OX   = T.("O3 WE")   - ELEC.OX_WE;    AE.OX   = T.("O3 AUX")  - ELEC.OX_AUX;
WE.NO2  = T.("NO2 WE")  - ELEC.NO2_WE;   AE.NO2  = T.("NO2 AUX") - ELEC.NO2_AUX;

% ---------------- Temperature compensation tables ----------------
T_pts   = [-30 -20 -10  0  10  20  30  40  50];

% NO (kT, ratio model)
kT_NO   = interp1(T_pts, [1.8 1.8 1.4 1.1 1.1 1.0 0.9 0.9 0.8], TempC, 'linear','extrap');

% NO2 & OX (nT, simple model) – replace OX row with your own if you have it
nT_NO2  = interp1(T_pts, [1.3 1.3 1.3 1.3 1.0 0.6 0.4 0.2 -1.5], TempC, 'linear','extrap');
nT_OX   = interp1(T_pts, [1.3 1.3 1.3 1.3 1.0 0.6 0.4 0.2 -1.5], TempC, 'linear','extrap');

% ---------------- Apply algorithms ----------------
% NO (Algorithm 2 / ratio): WEc = (WEu - WEe) - kT*(WEo/AEo)*(AEu - AEe)
WEc_NO   = (WE.NO) - kT_NO .* (ZERO.NO_WE ./ZERO.NO_AUX ) .* (AE.NO );
NO_NET = T.("NO WE") -T.("NO AUX") ; 

% NO2 (Algorithm 1 / nT):   WEc = (WEu - WEe) - nT*(AEu - AEe)
%WEc_NO2  = (WE.NO2 - WEo.NO2) - nT_NO2 .* (AE.NO2 - AEo.NO2);
WEc_NO2  = (WE.NO2) - nT_NO2 .* (AE.NO2);
NO2_NET = T.("NO2 WE") - T.("NO2 AUX");

% OX (Algorithm 1 / nT):    WEc = (WEu - WEe) - nT*(AEu - AEe)
WEc_OX   = (WE.OX ) - nT_OX  .* (AE.OX );
OX_NET = T.("O3 WE") - T.("O3 AUX");

% ---------------- Plot sets ----------------
makeFig("NO (k_T ratio)",  time, TempC, Hum, T.("NO WE") ,  T.("NO AUX"),  WEc_NO,NO_NET);
makeFig("NO2 (n_T simple)",time, TempC, Hum, T.("NO2 WE") , T.("NO2 AUX"), WEc_NO2,NO2_NET);
makeFig("OX/O3 (n_T simple)",time,TempC, Hum, T.("O3 WE") ,  T.("O3 AUX"), WEc_OX,OX_NET);

%% --------------- Plot helper ---------------
function makeFig(name,t,tc,hum,WE,AUX,WEc,NET)
    f = figure('Name',name,'Color','w');
    tl = tiledlayout(f,2,1,'TileSpacing','compact','Padding','compact');

    % Top: Temperature & Humidity
    ax1 = nexttile(tl);
    yyaxis left,  plot(t,tc,'LineWidth',1.2), ylabel('Temperature (°C)')
    yyaxis right, plot(t,hum,'LineWidth',1.2), ylabel('Humidity (%)')
    grid on, title([name,' — Temp & Humidity'])

    % Bottom: WE/AUX vs Corrected
    ax2 = nexttile(tl);
    yyaxis left
    plot(t,WE,'LineWidth',1.2); hold on
    plot(t,AUX,'LineWidth',1.2); ylabel('Raw Sensor (V)')
    yyaxis right
    plot(t,WEc,'LineWidth',1.5); ylabel('Corrected WE_c (V)'); hold on
    %plot(t,NET,'LineWidth',1.5); 
    plot(t, NET, 'LineWidth', 1.5, 'Color', 'g', 'LineStyle', '-');
    xlabel('Time'), grid on
    title([name,' — WE/AUX (left), WE_c (right)'])
    legend({'WE','AUX','WE Corrected*(Alpha Model)','Net Corrected'},'Location','best')
    linkaxes([ax1 ax2],'x');

    %f2 = figure('Name',[name ' — NET vs Temperature'],'Color','w');
    
       % --- Time window ---
    t0 = datetime(2025,10,30,11,41,0);
    t1 = datetime(2025,10,30,15,05,0);
    
    % --- Filter (and drop NaNs/Infs) ---
    mask = (t >= t0) & (t <= t1) & isfinite(tc) & isfinite(NET);
    tc_f  = tc(mask);
    NET_f = NET(mask);
    WE_f = WE(mask);
    
% --- Plot scatter ---
figure('Color','w');
%scatter(tc_f, NET_f, 5, 'filled'); grid on; hold on
scatter(tc_f, WE_f, 5, 'filled'); grid on; hold on
xlabel('Temperature (°C)');
ylabel('NET (V)');
title("NET vs Temperature — " + string(name));

% --- Polynomial fit (choose degree) ---
deg = 3;                                        % << change degree here
deg = min([deg, numel(tc_f)-1, 6]);             % guard: enough points & avoid wild fits

% Fit and trend line
p    = polyfit(tc_f(:), WE_f(:), deg);         % coefficients, highest power first
xfit = linspace(min(tc_f), max(tc_f), 400);
yfit = polyval(p, xfit);
plot(xfit, yfit, 'LineWidth', 1.6);

% R^2
yhat  = polyval(p, tc_f(:));
SSres = sum((WE_f(:) - yhat).^2);
SStot = sum((WE_f(:) - mean(WE_f(:))).^2);
R2 = 1 - SSres/SStot;

% Legend: show full equation for deg <= 3, otherwise compact label
if deg == 1
    eqStr = sprintf('y = %.3g x %+ .3g  (R^2 = %.3f)', p(1), p(2), R2);
elseif deg == 2
    a=p(1); b=p(2); c=p(3);
    eqStr = sprintf('y = %.3g x^2 %+ .3g x %+ .3g  (R^2 = %.3f)', a,b,c,R2);
elseif deg == 3
    a=p(1); b=p(2); c=p(3); d=p(4);
    eqStr = sprintf('y = %.3g x^3 %+ .3g x^2 %+ .3g x %+ .3g  (R^2 = %.3f)', a,b,c,d,R2);
else
    eqStr = sprintf('Poly%d fit (R^2 = %.3f)', deg, R2);
end
legend('Data', eqStr, 'Location','best');




end

[Star Strider]

I have been thinking about this for a few hours.

You need to capture the transients, and a third-order polynomial fit is not going to do it. I was able to get a reasonably decent fit of 'WE_f' using 'TempC' and 'Absolute Humidity' as inputs. (The polyfit and polyval functions can do multiple linear regression with two or more inputs, however it is inconvenient. There are other Toolboxes that offer this capability, however linear regression with even a high-order polynomial are not the best options here.)

I included a sample of that estimation at the end (as previously). It is difficult for me to go through your code and then regenerate your variables using my resampled data (required for this), so I leave that to you for the time being. You can use the System Identification Toolbox here, assuming that you do not have it licensed and installed on your own computer.

Also, since 'SensorLog_ALLS2.mat' was not provided, I used 'SensorLog_ALLS1.mat' that was.

Try this --

% matFile = fullfile(folder, 'SensorLog_ALLS2.mat');
dir

.                    ..                   SensorLog_ALLS1.mat  

matFile = 'SensorLog_ALLS1.mat';

% --- Load the combined table ---
S = load(matFile);
fn = fieldnames(S);
T  = S.(fn{find(structfun(@istable, S), 1)});

% --- Time & environment ---
if ~isdatetime(T.Timestamp)
    T.Timestamp = datetime(T.Timestamp,'InputFormat','yyyy-MM-dd HH:mm:ss');
end
time  = T.Timestamp;
TempC = T.("Temperature (°C)");
Hum   = T.("Humidity (%)");

% ================== CALIBRATION CONSTANTS ==================
% ---------- ELECTRONICS OFFSETS (hardware bias; set if measured) ----------
ELEC.NO_WE   = 0.292;  ELEC.NO_AUX   = 0.259;
ELEC.OX_WE   = 0.229;  ELEC.OX_AUX   = 0.227;   % OX uses O3 columns
ELEC.NO2_WE  = 0.238;  ELEC.NO2_AUX  = 0.225;

% ---------- SENSOR ZERO BASELINES (clean air, include electronics if measured that way) ----------
ZERO.NO_WE   = 0.3158;  ZERO.NO_AUX   = 0.3017;   % NO-B4
ZERO.OX_WE   = 0.2366;  ZERO.OX_AUX   = 0.2278;   % OX-A431 (O3)
ZERO.NO2_WE  = 0.2396;  ZERO.NO2_AUX  = 0.2253;   % NO2-B43F
% ===========================================================

% --- Ensure expected columns exist (adjust here if names differ) ---
need = ["NO WE","NO AUX","O3 WE","O3 AUX","NO2 WE","NO2 AUX"];
have = string(T.Properties.VariableNames);
if ~all(ismember(need, have))
    error('Missing columns. Expected: %s\nHave: %s', strjoin(need,', '), strjoin(have,', '));
end

% --- Initialise structs to avoid dot-indexing errors ---
WE = struct(); AE = struct(); WEo = struct(); AEo = struct();

% --- Convert zeros to sensor-only baselines (remove electronics bias) ---
WEo.NO  = ZERO.NO_WE  - ELEC.NO_WE;   AEo.NO  = ZERO.NO_AUX  - ELEC.NO_AUX;
WEo.OX  = ZERO.OX_WE  - ELEC.OX_WE;   AEo.OX  = ZERO.OX_AUX  - ELEC.OX_AUX;
WEo.NO2 = ZERO.NO2_WE - ELEC.NO2_WE;  AEo.NO2 = ZERO.NO2_AUX - ELEC.NO2_AUX;

% --- Raw -> electronics-corrected live readings ---
WE.NO   = T.("NO WE")   - ELEC.NO_WE;    AE.NO   = T.("NO AUX")  - ELEC.NO_AUX;
WE.OX   = T.("O3 WE")   - ELEC.OX_WE;    AE.OX   = T.("O3 AUX")  - ELEC.OX_AUX;
WE.NO2  = T.("NO2 WE")  - ELEC.NO2_WE;   AE.NO2  = T.("NO2 AUX") - ELEC.NO2_AUX;

% ---------------- Temperature compensation tables ----------------
T_pts   = [-30 -20 -10  0  10  20  30  40  50];

% NO (kT, ratio model)
kT_NO   = interp1(T_pts, [1.8 1.8 1.4 1.1 1.1 1.0 0.9 0.9 0.8], TempC, 'linear','extrap');

% NO2 & OX (nT, simple model) – replace OX row with your own if you have it
nT_NO2  = interp1(T_pts, [1.3 1.3 1.3 1.3 1.0 0.6 0.4 0.2 -1.5], TempC, 'linear','extrap');
nT_OX   = interp1(T_pts, [1.3 1.3 1.3 1.3 1.0 0.6 0.4 0.2 -1.5], TempC, 'linear','extrap');

% ---------------- Apply algorithms ----------------
% NO (Algorithm 2 / ratio): WEc = (WEu - WEe) - kT*(WEo/AEo)*(AEu - AEe)
WEc_NO   = (WE.NO) - kT_NO .* (ZERO.NO_WE ./ZERO.NO_AUX ) .* (AE.NO );
NO_NET = T.("NO WE") -T.("NO AUX") ; 

% NO2 (Algorithm 1 / nT):   WEc = (WEu - WEe) - nT*(AEu - AEe)
%WEc_NO2  = (WE.NO2 - WEo.NO2) - nT_NO2 .* (AE.NO2 - AEo.NO2);
WEc_NO2  = (WE.NO2) - nT_NO2 .* (AE.NO2);
NO2_NET = T.("NO2 WE") - T.("NO2 AUX");

% OX (Algorithm 1 / nT):    WEc = (WEu - WEe) - nT*(AEu - AEe)
WEc_OX   = (WE.OX ) - nT_OX  .* (AE.OX );
OX_NET = T.("O3 WE") - T.("O3 AUX");

% ---------------- Plot sets ----------------
makeFig("NO (k_T ratio)",  time, TempC, Hum, T.("NO WE") ,  T.("NO AUX"),  WEc_NO,NO_NET);

makeFig("NO2 (n_T simple)",time, TempC, Hum, T.("NO2 WE") , T.("NO2 AUX"), WEc_NO2,NO2_NET);

makeFig("OX/O3 (n_T simple)",time,TempC, Hum, T.("O3 WE") ,  T.("O3 AUX"), WEc_OX,OX_NET);

%% --------------- Plot helper ---------------
function [tc_f,NET_f,WE_f] = makeFig(name,t,tc,hum,WE,AUX,WEc,NET)
    f = figure('Name',name,'Color','w');
    tl = tiledlayout(f,2,1,'TileSpacing','compact','Padding','compact');

    % Top: Temperature & Humidity
    ax1 = nexttile(tl);
    yyaxis left,  plot(t,tc,'LineWidth',1.2), ylabel('Temperature (°C)')
    yyaxis right, plot(t,hum,'LineWidth',1.2), ylabel('Humidity (%)')
    grid on, title([name,' — Temp & Humidity'])

    % Bottom: WE/AUX vs Corrected
    ax2 = nexttile(tl);
    yyaxis left
    plot(t,WE,'LineWidth',1.2); hold on
    plot(t,AUX,'LineWidth',1.2); ylabel('Raw Sensor (V)')
    yyaxis right
    plot(t,WEc,'LineWidth',1.5); ylabel('Corrected WE_c (V)'); hold on
    %plot(t,NET,'LineWidth',1.5); 
    plot(t, NET, 'LineWidth', 1.5, 'Color', 'g', 'LineStyle', '-');
    xlabel('Time'), grid on
    title([name,' — WE/AUX (left), WE_c (right)'])
    legend({'WE','AUX','WE Corrected*(Alpha Model)','Net Corrected'},'Location','best')
    linkaxes([ax1 ax2],'x');

    %f2 = figure('Name',[name ' — NET vs Temperature'],'Color','w');
    
       % --- Time window ---
    t0 = datetime(2025,10,30,11,41,0);
    t1 = datetime(2025,10,30,15,05,0);
    
    % --- Filter (and drop NaNs/Infs) ---
    mask = (t >= t0) & (t <= t1) & isfinite(tc) & isfinite(NET);
    tc_f  = tc(mask);
    NET_f = NET(mask);
    WE_f = WE(mask);
    
% --- Plot scatter ---
figure('Color','w');
%scatter(tc_f, NET_f, 5, 'filled'); grid on; hold on
scatter(tc_f, WE_f, 5, 'filled'); grid on; hold on
xlabel('Temperature (°C)');
ylabel('NET (V)');
title("NET vs Temperature — " + string(name));

% --- Polynomial fit (choose degree) ---
deg = 3;                                        % << change degree here
deg = min([deg, numel(tc_f)-1, 6]);             % guard: enough points & avoid wild fits

% Fit and trend line
p    = polyfit(tc_f(:), WE_f(:), deg);         % coefficients, highest power first
xfit = linspace(min(tc_f), max(tc_f), 400);
yfit = polyval(p, xfit);
plot(xfit, yfit, 'LineWidth', 1.6);

% R^2
yhat  = polyval(p, tc_f(:));
SSres = sum((WE_f(:) - yhat).^2);
SStot = sum((WE_f(:) - mean(WE_f(:))).^2);
R2 = 1 - SSres/SStot;

% Legend: show full equation for deg <= 3, otherwise compact label
if deg == 1
    eqStr = sprintf('y = %.3g x %+ .3g  (R^2 = %.3f)', p(1), p(2), R2);
elseif deg == 2
    a=p(1); b=p(2); c=p(3);
    eqStr = sprintf('y = %.3g x^2 %+ .3g x %+ .3g  (R^2 = %.3f)', a,b,c,R2);
elseif deg == 3
    a=p(1); b=p(2); c=p(3); d=p(4);
    eqStr = sprintf('y = %.3g x^3 %+ .3g x^2 %+ .3g x %+ .3g  (R^2 = %.3f)', a,b,c,d,R2);
else
    eqStr = sprintf('Poly%d fit (R^2 = %.3f)', deg, R2);
end
legend('Data', eqStr, 'Location','best');




end


% ===== System Identification Start =====
Ts = mean(diff(T.Timestamp));
Ts.Format = 'hh:mm:ss.SSS'

Ts = duration
   00:00:05.513

Tsd = std(diff(T.Timestamp));
Tsd.Format = 'hh:mm:ss.SSS'

Tsd = duration
   00:00:01.229

Tss = seconds(Ts)

Tss = 5.5140

Fs = 1/Tss                                          % Sampling Frequency (Hz)

Fs = 0.1814

% Source Habs: https://www.reddit.com/r/embedded/comments/f4iino/whats_the_right_way_to_convert_relative_humidity/
Habs = @(RH,T) 216.7*(RH/100 .* 6.112 .* exp((17.62.*T)./(243.12+T)))./(273.15+T);  %Units: g/m^3
HumAbs = Habs(T.('Humidity (%)'), T.('Temperature (°C)'));

Tr1 = T;
Tr1 = addvars(Tr1, HumAbs, 'Before','NO WE', NewVariableNames='Absolute Humidity');
TTr1 = table2timetable(Tr1(:,1:end-1));
fprintf('\nResampled Data --\n')

Resampled Data --

TTr2 = retime(TTr1, 'regular', 'linear', 'TimeStep',seconds(Tss))

TTr2 = 2683×9 timetable
         Timestamp         Temperature (°C)    Humidity (%)    Absolute Humidity     NO WE     NO AUX      O3 WE     O3 AUX     NO2 WE     NO2 AUX
    ___________________    ________________    ____________    _________________    _______    _______    _______    _______    _______    _______

    2025-10-30 11:04:49         9.6125            55.839             5.112          0.25413      0.248      0.238    0.23587    0.26687    0.23513
    2025-10-30 11:04:55         9.5463            55.651            5.0736          0.25523      0.248      0.238    0.23477    0.26577    0.23623
    2025-10-30 11:05:00         9.4835             55.46            5.0361          0.25567      0.248      0.238    0.23433    0.26533      0.237
    2025-10-30 11:05:06          9.411            55.197            4.9891          0.25543      0.248      0.238      0.235      0.266      0.237
    2025-10-30 11:05:11         9.3385             54.93            4.9421            0.256      0.248    0.23854      0.235    0.26546      0.237
    2025-10-30 11:05:17         9.2987            54.785            4.9166            0.256      0.248      0.239      0.235      0.265      0.237
    2025-10-30 11:05:22         9.2583            54.655            4.8923          0.25708    0.24854      0.239      0.235    0.26554    0.23754
    2025-10-30 11:05:28         9.2142            54.516            4.8662            0.258      0.249      0.239      0.235      0.266      0.238
    2025-10-30 11:05:33         9.1776            54.373            4.8421          0.25875    0.24975    0.23975      0.235    0.26525      0.238
    2025-10-30 11:05:39         9.1445            54.229            4.8191            0.259       0.25       0.24    0.23585      0.265    0.23715
    2025-10-30 11:05:44         9.1019            54.086            4.7933          0.25995       0.25       0.24      0.236      0.265      0.237
    2025-10-30 11:05:50         9.0683            53.957            4.7717          0.26006       0.25    0.24006    0.23506      0.265    0.23706
    2025-10-30 11:05:55          9.036            53.913            4.7579          0.26113    0.25013    0.24113      0.236      0.265      0.238
    2025-10-30 11:06:01         9.0075            53.926            4.7504            0.262    0.25106    0.24194      0.236    0.26494      0.238
    2025-10-30 11:06:06         8.9634            53.819            4.7277            0.262      0.252    0.24067      0.236    0.26384      0.238
    2025-10-30 11:06:12          8.914            53.484            4.6835          0.26173    0.25173     0.2382      0.236    0.26273      0.238

% y = TTr2.('NO WE');
% u = TTr2.('Absolute Humidity');

t0 = datetime(2025,10,30,11,41,0);
t1 = datetime(2025,10,30,15,05,0);

time  = TTr2.Timestamp;
TempC = TTr2.("Temperature (°C)");
Hum   = TTr2.("Humidity (%)");
AHum = TTr2.('Absolute Humidity');

NO_NET = TTr2.("NO WE") -TTr2.("NO AUX") ; 
NO2_NET = TTr2.("NO2 WE") - TTr2.("NO2 AUX");
OX_NET = TTr2.("O3 WE") - TTr2.("O3 AUX");

NET = NO_NET;
mask = (time >= t0) & (time <= t1) & isfinite(TTr2{:,2}) & isfinite(NET);
tc_f  = TempC(mask);
NET_f = NO_NET(mask);
WE_f = TTr2.("NO WE")(mask);
A_Hum = AHum(mask);

% u = [tc_f(:) NET_f];
u = [tc_f A_Hum];                           % System Input
y = WE_f(:);                                % System Output
Ts = Tss                                    % Sampling Time (Sampling Interval)

Ts = 5.5140

tc_WE_Data = iddata(y,u,Ts)

tc_WE_Data = 
Time domain data set with 2220 samples.
Sample time: 5.51399 seconds           
                                       
Outputs      Unit (if specified)       
   y1                                  
                                       
Inputs       Unit (if specified)       
   u1                                  
   u2                                  
                                       

SysOrd = 8;

tc_WE_Sys = ssest(tc_WE_Data,SysOrd)

tc_WE_Sys =
  Continuous-time identified state-space model:
      dx/dt = A x(t) + B u(t) + K e(t)
       y(t) = C x(t) + D u(t) + e(t)
 
  A = 
               x1          x2          x3          x4          x5          x6          x7          x8
   x1  -3.843e-05   0.0001523   0.0004896   0.0004103   8.043e-05    -0.00102   0.0001702    0.001481
   x2   -0.001303    0.002956       -0.13    -0.04106     0.01189     0.06043    -0.02728    -0.09427
   x3  -0.0003734      0.1073    -0.04462    -0.05033    -0.02847      0.1011   -0.005904     -0.1543
   x4  -0.0008889     0.02689    -0.01625    -0.02561    -0.03728     0.08115    -0.01322      -0.131
   x5  -0.0008022    0.006755  -0.0008481     0.01861     0.02526      0.2683    -0.04223     -0.2239
   x6   0.0006724   -0.006607     0.04552     0.03402      -0.186     -0.1156      -0.105      0.1796
   x7  -4.361e-05   -0.006482    -0.00212    0.001634    -0.02713     0.08044     0.03744       0.385
   x8   -0.002009     0.01199    -0.05217    -0.03259     0.09367      0.1406     -0.3246     -0.2297
 
  B = 
               u1          u2
   x1  -0.0001616   0.0001545
   x2     0.01481    -0.01094
   x3      0.0365    -0.01926
   x4     0.01656    -0.01341
   x5     0.04532    -0.01363
   x6    -0.06324     0.02734
   x7    -0.05442     0.04341
   x8     0.05403   -0.006948
 
  C = 
               x1          x2          x3          x4          x5          x6          x7          x8
   y1        5.17  -0.0006356   -0.004146   -0.001723    0.001816    0.001116   -0.001702   -0.001249
 
  D = 
       u1  u2
   y1   0   0
 
  K = 
            y1
   x1  0.04208
   x2  -0.4957
   x3   -2.569
   x4   -1.758
   x5  -0.5251
   x6     2.11
   x7   0.2163
   x8   -1.996
 
Parameterization:
   FREE form (all coefficients in A, B, C free).
   Feedthrough: none
   Disturbance component: estimate
   Number of free coefficients: 96
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                
Estimated using SSEST on time domain data "tc_WE_Data".
Fit to estimation data: 99.46% (prediction focus)      
FPE: 2.448e-07, MSE: 2.362e-07                         
 

figure
compare(tc_WE_Data, tc_WE_Sys)

% ===== System Identification End =====

Experiment with other orders as well. I could not get a good fit with a lower-order estimation, however 'good fit' is subjective.

.

[Dharmesh]

@Star Strider

Thank you. Just to confirm, did you mean that I need to provide some data that includes transient variations? Yes, the polynomial is simply the temperature correction for the baseline, prior to applying the humidity correction. I will capture some data and share my findings shortly.

[Star Strider]

My pleasure!

'Just to confirm, did you mean that I need to provide some data that includes transient variations?'

Yes.

Ideally, the input to the system identification functions needs to be the exact input data in 'u'. You then estimate the system by supplying the measured output in 'y'.

For example, I performed a system identification of a rat leg muscle (specifically rectus femoris) preparation a few decades ago. The input was a random stepwise variation of frequency and voltage (within limits) and the output was the displacement of a cantilever beam displacement transducer the muscle was tied to (and kept slightly stretched when relaxed). The input to the system identification function was the recorded stimulus signal in one channel, and the output was the simultaneously recorded deflection in another channel. I then did the system identification by using the stimulus signal as 'u' and the deflection as 'y'. I was able to get a decent fit to the data with a moderatey low-order state space realisation.

This is the sort of experiment you need to do. The input needs to be recorded (unlike my experiment, you will probably not be able to measure t, so use what you want to provide in terms of temperature and humidity steps as the input), and then measure the output. (I am not certain how the NOx signal fits with this, so it may be necessary to not record any NOx signals while you are varying the temperature and humidity.) That will be the system you want to estimate. You can then use the resulting model with the actual temperature and humidity data to produce a signal you can use to correct the observed output, probably by subtracting it from the observed output.

Again, I am not sufficiently familiar with your instrumentation (either theory or realisation) to specify an experimental design, so I must leave that to you.

EDIT -- Corrected typographical errors.

[Dharmesh]

@Star Strider

I have attached the new data. During this test, the temperature was held relatively stable, varying only between 18–20 °C, while the relative humidity fluctuated between 50 % and 80 %. You can see that changes in humidity have a noticeable effect on the WE signal of the sensor.

I’ve also attached my updated code below. I made some changes to handle the differences in timestamps between the datasets used for temperature compensation and humidity compensation. I was getting a few timestamp-related errors, so I added a variable temperature_test = false; to adjust the plots accordingly.

The expected WE signal outputs for each sensor under stable temperature and humidity conditions are approximately:

• NO: 322 mV
• NO₂: 228 mV
• OX: 218 mV

These reference values were derived from previous experiments using the data in SensorLog_ALLS2.mat, where we generated scatter plots and trend lines. Using the trend line at around 20 °C, I obtained the values above.

In the most recent dataset (at about 18–20 °C), we are achieving the following results, which are reasonably close to expectations. Any slight anomalies are likely due to environmental factors beyond our control:

• NO: 295 mV
• NO₂: 233 mV
• OX: 223 mV

These readings were taken during the period when humidity was relatively stable, between 14:00 and 14:24.

Ideally lets first focus on the NO₂ data.

%matFile = fullfile(folder, 'SensorLog_ALLS2.mat');

matFile = fullfile(folder, 'HumditydataS2.mat');

temperature_test = false; % set to true for the original dataset

% --- Load the combined table ---

S = load(matFile);

fn = fieldnames(S);

T = S.(fn{find(structfun(@istable, S), 1)});

% --- Time & environment ---

if ~isdatetime(T.Timestamp)

T.Timestamp = datetime(T.Timestamp,'InputFormat','yyyy-MM-dd HH:mm:ss');

end

time = T.Timestamp;

TempC = T.("Temperature (°C)");

Hum = T.("Humidity (%)");

% ================== CALIBRATION CONSTANTS ==================

% ---------- ELECTRONICS OFFSETS (hardware bias; set if measured) ----------

ELEC.NO_WE = 0.292; ELEC.NO_AUX = 0.259;

ELEC.OX_WE = 0.229; ELEC.OX_AUX = 0.227; % OX uses O3 columns

ELEC.NO2_WE = 0.238; ELEC.NO2_AUX = 0.225;

% ---------- SENSOR ZERO BASELINES (clean air, include electronics if measured that way) ----------

ZERO.NO_WE = 0.3158; ZERO.NO_AUX = 0.3017; % NO-B4

ZERO.OX_WE = 0.2366; ZERO.OX_AUX = 0.2278; % OX-A431 (O3)

ZERO.NO2_WE = 0.2396; ZERO.NO2_AUX = 0.2253; % NO2-B43F

% ===========================================================

% --- Ensure expected columns exist (adjust here if names differ) ---

need = ["NO WE","NO AUX","O3 WE","O3 AUX","NO2 WE","NO2 AUX"];

have = string(T.Properties.VariableNames);

if ~all(ismember(need, have))

error('Missing columns. Expected: %s\nHave: %s', strjoin(need,', '), strjoin(have,', '));

end

% --- Initialise structs to avoid dot-indexing errors ---

WE = struct(); AE = struct(); WEo = struct(); AEo = struct();

% --- Convert zeros to sensor-only baselines (remove electronics bias) ---

WEo.NO = ZERO.NO_WE - ELEC.NO_WE; AEo.NO = ZERO.NO_AUX - ELEC.NO_AUX;

WEo.OX = ZERO.OX_WE - ELEC.OX_WE; AEo.OX = ZERO.OX_AUX - ELEC.OX_AUX;

WEo.NO2 = ZERO.NO2_WE - ELEC.NO2_WE; AEo.NO2 = ZERO.NO2_AUX - ELEC.NO2_AUX;

% --- Raw -> electronics-corrected live readings ---

WE.NO = T.("NO WE") - ELEC.NO_WE; AE.NO = T.("NO AUX") - ELEC.NO_AUX;

WE.OX = T.("O3 WE") - ELEC.OX_WE; AE.OX = T.("O3 AUX") - ELEC.OX_AUX;

WE.NO2 = T.("NO2 WE") - ELEC.NO2_WE; AE.NO2 = T.("NO2 AUX") - ELEC.NO2_AUX;

% ---------------- Temperature compensation tables ----------------

T_pts = [-30 -20 -10 0 10 20 30 40 50];

% NO (kT, ratio model)

kT_NO = interp1(T_pts, [1.8 1.8 1.4 1.1 1.1 1.0 0.9 0.9 0.8], TempC, 'linear','extrap');

% NO2 & OX (nT, simple model) – replace OX row with your own if you have it

nT_NO2 = interp1(T_pts, [1.3 1.3 1.3 1.3 1.0 0.6 0.4 0.2 -1.5], TempC, 'linear','extrap');

nT_OX = interp1(T_pts, [1.3 1.3 1.3 1.3 1.0 0.6 0.4 0.2 -1.5], TempC, 'linear','extrap');

% ---------------- Apply algorithms ----------------

% NO (Algorithm 2 / ratio): WEc = (WEu - WEe) - kT*(WEo/AEo)*(AEu - AEe)

WEc_NO = (WE.NO) - kT_NO .* (ZERO.NO_WE ./ZERO.NO_AUX ) .* (AE.NO );

NO_NET = T.("NO WE") -T.("NO AUX") ;

% NO2 (Algorithm 1 / nT): WEc = (WEu - WEe) - nT*(AEu - AEe)

%WEc_NO2 = (WE.NO2 - WEo.NO2) - nT_NO2 .* (AE.NO2 - AEo.NO2);

WEc_NO2 = (WE.NO2) - nT_NO2 .* (AE.NO2);

NO2_NET = T.("NO2 WE") - T.("NO2 AUX");

% OX (Algorithm 1 / nT): WEc = (WEu - WEe) - nT*(AEu - AEe)

WEc_OX = (WE.OX ) - nT_OX .* (AE.OX );

OX_NET = T.("O3 WE") - T.("O3 AUX");

% ---------------- Plot sets ----------------

makeFig("NO (k_T ratio)", time, TempC, Hum, T.("NO WE") , T.("NO AUX"), WEc_NO,NO_NET,temperature_test);

makeFig("NO2 (n_T simple)",time, TempC, Hum, T.("NO2 WE") , T.("NO2 AUX"), WEc_NO2,NO2_NET,temperature_test);

makeFig("OX/O3 (n_T simple)",time,TempC, Hum, T.("O3 WE") , T.("O3 AUX"), WEc_OX,OX_NET,temperature_test);

%% --------------- Plot helper ---------------

function makeFig(name,t,tc,hum,WE,AUX,WEc,NET,temperature_test)

% if ~exist('temperature_test','var')

% temperature_test = false; % default

% end

f = figure('Name',name,'Color','w');

tl = tiledlayout(f,2,1,'TileSpacing','compact','Padding','compact');

% Top: Temperature & Humidity

ax1 = nexttile(tl);

yyaxis left, plot(t,tc,'LineWidth',1.2), ylabel('Temperature (°C)')

yyaxis right, plot(t,hum,'LineWidth',1.2), ylabel('Humidity (%)')

grid on, title([name,' — Temp & Humidity'])

% Bottom: WE/AUX vs Corrected

ax2 = nexttile(tl);

yyaxis left

if temperature_test

plot(t,WE,'LineWidth',1.2); hold on

plot(t,AUX,'LineWidth',1.2); ylabel('Raw Sensor (V)')

else

plot(t, WE, 'b', 'LineWidth', 1.2); hold on % blue line for WE

plot(t, AUX, 'r', 'LineWidth', 1.2); % red line for AUX

ylabel('Raw Sensor (V)');

end;

if temperature_test

yyaxis right

plot(t,WEc,'LineWidth',1.5); ylabel('Corrected WE_c (V)'); hold on

%plot(t,NET,'LineWidth',1.5);

plot(t, NET, 'LineWidth', 1.5, 'Color', 'g', 'LineStyle', '-');

end

xlabel('Time'), grid on

title([name,' — WE/AUX (left), WE_c (right)'])

if temperature_test

legend({'WE','AUX','WE Corrected*(Alpha Model)','Net Corrected'},'Location','best')

else

legend({'WE','AUX'},'Location','best')

end

linkaxes([ax1 ax2],'x');

%f2 = figure('Name',[name ' — NET vs Temperature'],'Color','w');

% Use fixed window for the original temperature test; otherwise use full span

if temperature_test

% --- Fixed window for original dataset ---

t0 = datetime(2025,10,30,11,41,0);

t1 = datetime(2025,10,30,15,05,0);

else

% --- Full-range window for current data (with a small padding) ---

tt = t(~isnat(t)); % remove NaT, if any

if isempty(tt)

warning('No valid timestamps in t; skipping.');

return

end

t0 = min(tt) - minutes(1);

t1 = max(tt) + minutes(1);

end

% --- Build one mask only ---

inwin = (t >= t0) & (t <= t1);

good = isfinite(tc) & isfinite(NET) & isfinite(WE); % include WE since you use WE_f too

mask = inwin & good;

% --- Filtered vectors ---

tc_f = tc(mask);

NET_f = NET(mask);

WE_f = WE(mask);

if temperature_test

% --- Plot scatter ---

figure('Color','w');

%scatter(tc_f, NET_f, 5, 'filled'); grid on; hold on

scatter(tc_f, WE_f, 5, 'filled'); grid on; hold on

xlabel('Temperature (°C)');

ylabel('NET (V)');

title("NET vs Temperature — " + string(name));

% --- Polynomial fit (choose degree) ---

deg = 3; % << change degree here

deg = min([deg, numel(tc_f)-1, 6]); % guard: enough points & avoid wild fits

% Fit and trend line

p = polyfit(tc_f(:), WE_f(:), deg); % coefficients, highest power first

xfit = linspace(min(tc_f), max(tc_f), 400);

yfit = polyval(p, xfit);

plot(xfit, yfit, 'LineWidth', 1.6);

% R^2

yhat = polyval(p, tc_f(:));

SSres = sum((WE_f(:) - yhat).^2);

SStot = sum((WE_f(:) - mean(WE_f(:))).^2);

R2 = 1 - SSres/SStot;

% Legend: show full equation for deg <= 3, otherwise compact label

if deg == 1

eqStr = sprintf('y = %.3g x %+ .3g (R^2 = %.3f)', p(1), p(2), R2);

elseif deg == 2

a=p(1); b=p(2); c=p(3);

eqStr = sprintf('y = %.3g x^2 %+ .3g x %+ .3g (R^2 = %.3f)', a,b,c,R2);

elseif deg == 3

a=p(1); b=p(2); c=p(3); d=p(4);

eqStr = sprintf('y = %.3g x^3 %+ .3g x^2 %+ .3g x %+ .3g (R^2 = %.3f)', a,b,c,d,R2);

else

eqStr = sprintf('Poly%d fit (R^2 = %.3f)', deg, R2);

end

legend('Data', eqStr, 'Location','best');

end

end

[Star Strider]

I have copied this over to MATLAB Online for reference and further exploration. (I changed 'makeFig' so that it returns the variables of interest in elements of cell arrays to make this easier for me.)

I re-read this several times in the last few days, however I still do not understand it as well as I would like. The way I am now understanding this is that your instrumentation outputs voltages for each sensor, and you want to correct the NO and NOx voltages for the temperature and humidity signals, and that this correction has to allow for transients that can be significant and usually return to the previous baseline in about 10 minutes. I assume that you are testing the NO sensor with a known NO concentration, and that you want to understand how to model the temperature and humidity (absolute humidity is what the literature that I found seems to recommend) to correct the NO or NOx or OX signals.

The signals from about 14:00 to 14:24 show an interesting decaying exponential transient in the NO2 and OX responses to the step-wise temperature and humidity fluctuations. If temperature and humidity are the inputs, absolute humidity can be calculated. When we model this using the System Identification Toolbox, what input and output variables should the signals be?

Right now, I would model it with temperature and humidity as the inputs (absolute humidity is not linearly related to them), and NO or NOx (respectively) as the output. I am not at all certain that this is correct, or would give the desired information.

We need inputs and known outputs in order to correctly model the system. The System Identification Toolbox is necessary since there are transients involved, and linear or nonlinear regression will not be capable of accurately describing the transients. Currently, I do not understand how the various signals combine in order to produce an accurate NO, NOx or OX output.

I will continue with this once I understand what we want to model.

files = dir('*.mat');
folder = files.folder;

%matFile = fullfile(folder, 'SensorLog_ALLS2.mat');
matFile = fullfile(folder, 'HumditydataS2.mat');
temperature_test = false;   % set to true for the original dataset
% --- Load the combined table ---
S = load(matFile);
fn = fieldnames(S);
T  = S.(fn{find(structfun(@istable, S), 1)});
% --- Time & environment ---
if ~isdatetime(T.Timestamp)
    T.Timestamp = datetime(T.Timestamp,'InputFormat','yyyy-MM-dd HH:mm:ss');
end
time  = T.Timestamp;
TempC = T.("Temperature (°C)");
Hum   = T.("Humidity (%)");
% ================== CALIBRATION CONSTANTS ==================
% ---------- ELECTRONICS OFFSETS (hardware bias; set if measured) ----------
ELEC.NO_WE   = 0.292;  ELEC.NO_AUX   = 0.259;
ELEC.OX_WE   = 0.229;  ELEC.OX_AUX   = 0.227;   % OX uses O3 columns
ELEC.NO2_WE  = 0.238;  ELEC.NO2_AUX  = 0.225;
% ---------- SENSOR ZERO BASELINES (clean air, include electronics if measured that way) ----------
ZERO.NO_WE   = 0.3158;  ZERO.NO_AUX   = 0.3017;   % NO-B4
ZERO.OX_WE   = 0.2366;  ZERO.OX_AUX   = 0.2278;   % OX-A431 (O3)
ZERO.NO2_WE  = 0.2396;  ZERO.NO2_AUX  = 0.2253;   % NO2-B43F
% ===========================================================
% --- Ensure expected columns exist (adjust here if names differ) ---
need = ["NO WE","NO AUX","O3 WE","O3 AUX","NO2 WE","NO2 AUX"];
have = string(T.Properties.VariableNames);
if ~all(ismember(need, have))
    error('Missing columns. Expected: %s\nHave: %s', strjoin(need,', '), strjoin(have,', '));
end
% --- Initialise structs to avoid dot-indexing errors ---
WE = struct(); AE = struct(); WEo = struct(); AEo = struct();
% --- Convert zeros to sensor-only baselines (remove electronics bias) ---
WEo.NO  = ZERO.NO_WE  - ELEC.NO_WE;   AEo.NO  = ZERO.NO_AUX  - ELEC.NO_AUX;
WEo.OX  = ZERO.OX_WE  - ELEC.OX_WE;   AEo.OX  = ZERO.OX_AUX  - ELEC.OX_AUX;
WEo.NO2 = ZERO.NO2_WE - ELEC.NO2_WE;  AEo.NO2 = ZERO.NO2_AUX - ELEC.NO2_AUX;
% --- Raw -> electronics-corrected live readings ---
WE.NO   = T.("NO WE")   - ELEC.NO_WE;    AE.NO   = T.("NO AUX")  - ELEC.NO_AUX;
WE.OX   = T.("O3 WE")   - ELEC.OX_WE;    AE.OX   = T.("O3 AUX")  - ELEC.OX_AUX;
WE.NO2  = T.("NO2 WE")  - ELEC.NO2_WE;   AE.NO2  = T.("NO2 AUX") - ELEC.NO2_AUX;
% ---------------- Temperature compensation tables ----------------
T_pts   = [-30 -20 -10  0  10  20  30  40  50];
% NO (kT, ratio model)
kT_NO   = interp1(T_pts, [1.8 1.8 1.4 1.1 1.1 1.0 0.9 0.9 0.8], TempC, 'linear','extrap');
% NO2 & OX (nT, simple model) – replace OX row with your own if you have it
nT_NO2  = interp1(T_pts, [1.3 1.3 1.3 1.3 1.0 0.6 0.4 0.2 -1.5], TempC, 'linear','extrap');
nT_OX   = interp1(T_pts, [1.3 1.3 1.3 1.3 1.0 0.6 0.4 0.2 -1.5], TempC, 'linear','extrap');
% ---------------- Apply algorithms ----------------
% NO (Algorithm 2 / ratio): WEc = (WEu - WEe) - kT*(WEo/AEo)*(AEu - AEe)
WEc_NO   = (WE.NO) - kT_NO .* (ZERO.NO_WE ./ZERO.NO_AUX ) .* (AE.NO );
NO_NET = T.("NO WE") -T.("NO AUX") ; 
% NO2 (Algorithm 1 / nT):   WEc = (WEu - WEe) - nT*(AEu - AEe)
%WEc_NO2  = (WE.NO2 - WEo.NO2) - nT_NO2 .* (AE.NO2 - AEo.NO2);
WEc_NO2  = (WE.NO2) - nT_NO2 .* (AE.NO2);
NO2_NET = T.("NO2 WE") - T.("NO2 AUX");
% OX (Algorithm 1 / nT):    WEc = (WEu - WEe) - nT*(AEu - AEe)
WEc_OX   = (WE.OX ) - nT_OX  .* (AE.OX );
OX_NET = T.("O3 WE") - T.("O3 AUX");
% ---------------- Plot sets ----------------
[tc_f{1},NET_f{1},WE_f{1}] = makeFig("NO (k_T ratio)",  time, TempC, Hum, T.("NO WE") ,  T.("NO AUX"),  WEc_NO,NO_NET,temperature_test);

[tc_f{2},NET_f{2},WE_f{2}] = makeFig("NO2 (n_T simple)",time, TempC, Hum, T.("NO2 WE") , T.("NO2 AUX"), WEc_NO2,NO2_NET,temperature_test);

[tc_f{3},NET_f{3},WE_f{3}] = makeFig("OX/O3 (n_T simple)",time,TempC, Hum, T.("O3 WE") ,  T.("O3 AUX"), WEc_OX,OX_NET,temperature_test);

%% --------------- Plot helper ---------------
function [tc_f,NET_f,WE_f] = makeFig(name,t,tc,hum,WE,AUX,WEc,NET,temperature_test)
   % if ~exist('temperature_test','var')
   %              temperature_test = false;   % default
   % end
    f = figure('Name',name,'Color','w');
    tl = tiledlayout(f,2,1,'TileSpacing','compact','Padding','compact');
    % Top: Temperature & Humidity
    ax1 = nexttile(tl);
    yyaxis left,  plot(t,tc,'LineWidth',1.2), ylabel('Temperature (°C)')
    yyaxis right, plot(t,hum,'LineWidth',1.2), ylabel('Humidity (%)')
    grid on, title([name,' — Temp & Humidity'])
    % Bottom: WE/AUX vs Corrected
    ax2 = nexttile(tl);
    yyaxis left
    if temperature_test
    plot(t,WE,'LineWidth',1.2); hold on
    plot(t,AUX,'LineWidth',1.2); ylabel('Raw Sensor (V)')
    else
     plot(t, WE, 'b', 'LineWidth', 1.2); hold on      % blue line for WE
     plot(t, AUX, 'r', 'LineWidth', 1.2);             % red line for AUX
     ylabel('Raw Sensor (V)');
    end;
    
    if temperature_test
        yyaxis right
        plot(t,WEc,'LineWidth',1.5); ylabel('Corrected WE_c (V)'); hold on
        %plot(t,NET,'LineWidth',1.5); 
        plot(t, NET, 'LineWidth', 1.5, 'Color', 'g', 'LineStyle', '-');
    end
    xlabel('Time'), grid on
    title([name,' — WE/AUX (left), WE_c (right)'])
    
    if temperature_test
        legend({'WE','AUX','WE Corrected*(Alpha Model)','Net Corrected'},'Location','best')
    else
        legend({'WE','AUX'},'Location','best')
    end
    linkaxes([ax1 ax2],'x');
    %f2 = figure('Name',[name ' — NET vs Temperature'],'Color','w');
            % Use fixed window for the original temperature test; otherwise use full span
         
               
            if temperature_test 
                % --- Fixed window for original dataset ---
                t0 = datetime(2025,10,30,11,41,0);
                t1 = datetime(2025,10,30,15,05,0);
            else
                % --- Full-range window for current data (with a small padding) ---
                tt = t(~isnat(t));                 % remove NaT, if any
                if isempty(tt)
                    warning('No valid timestamps in t; skipping.');
                    return
                end
                t0 = min(tt) - minutes(1);
                t1 = max(tt) + minutes(1);
            end
            
            % --- Build one mask only ---
            inwin = (t >= t0) & (t <= t1);
            good  = isfinite(tc) & isfinite(NET) & isfinite(WE);   % include WE since you use WE_f too
            mask  = inwin & good;
            
            % --- Filtered vectors ---
            tc_f  = tc(mask);
            NET_f = NET(mask);
            WE_f  = WE(mask);
            
        if temperature_test 
                % --- Plot scatter ---
                figure('Color','w');
                %scatter(tc_f, NET_f, 5, 'filled'); grid on; hold on
                scatter(tc_f, WE_f, 5, 'filled'); grid on; hold on
                xlabel('Temperature (°C)');
                ylabel('NET (V)');
                title("NET vs Temperature — " + string(name));
                
                % --- Polynomial fit (choose degree) ---
                deg = 3;                                        % << change degree here
                deg = min([deg, numel(tc_f)-1, 6]);             % guard: enough points & avoid wild fits
                
                % Fit and trend line
                p    = polyfit(tc_f(:), WE_f(:), deg);         % coefficients, highest power first
                xfit = linspace(min(tc_f), max(tc_f), 400);
                yfit = polyval(p, xfit);
                plot(xfit, yfit, 'LineWidth', 1.6);
                
                % R^2
                yhat  = polyval(p, tc_f(:));
                SSres = sum((WE_f(:) - yhat).^2);
                SStot = sum((WE_f(:) - mean(WE_f(:))).^2);
                R2 = 1 - SSres/SStot;
                
                % Legend: show full equation for deg <= 3, otherwise compact label
                if deg == 1
                    eqStr = sprintf('y = %.3g x %+ .3g  (R^2 = %.3f)', p(1), p(2), R2);
                elseif deg == 2
                    a=p(1); b=p(2); c=p(3);
                    eqStr = sprintf('y = %.3g x^2 %+ .3g x %+ .3g  (R^2 = %.3f)', a,b,c,R2);
                elseif deg == 3
                    a=p(1); b=p(2); c=p(3); d=p(4);
                    eqStr = sprintf('y = %.3g x^3 %+ .3g x^2 %+ .3g x %+ .3g  (R^2 = %.3f)', a,b,c,d,R2);
                else
                    eqStr = sprintf('Poly%d fit (R^2 = %.3f)', deg, R2);
                end
                legend('Data', eqStr, 'Location','best');
        end
end

EDIT -- (11 Nov 2025 at 3:11)

I envision something like this --

To start with. compare with a known concentration of the analyte of interest (here, NOx, can be any). There can be more inputs and outputs, this is a schematic of the simplest system. All of these are vectors as functions of time analysed at the same times.

Signals Known Values

---------------------

T(ºC) --> | |

| |

Habs --> | System | --> NOx

| |

NOx --> | |

---------------------

The 'Signals' are from the instrumentation. 'Habs' is the absolute humidity.

The 'System' is then estimated.

With that information, the analysis with the estimated 'System' becomes --

Signals Correct Values

---------------------

T(ºC) --> | |

| |

Habs --> | System | --> NOx

| |

NOx --> | |

---------------------

This is the sort of approach I believe we need to do.

.

[Dharmesh]

Hi,

Thank you for your help.

Regarding your understanding — yes, each sensor outputs a voltage that is affected by both temperature and humidity. There are transients caused independently by each of these factors, but I don’t necessarily know if the stabilisation time is 10 minutes, less, or more, as this depends on the characteristics of the sensor.

At the moment, in relation to gas concentration, we are testing in zero gas, which means we are simply establishing a corrected baseline. The baseline signifies the voltage at each temperature step when the gas concentration is zero. When we subtract this baseline from the actual signal (with gas present) and divide by the manufacturer’s sensitivity factor, we obtain the actual gas concentration. So at this stage, we are focusing only on the baseline to better understand how humidity and temperature affect the voltage output.

The period between 14:00 and 14:24 is the most stable, and I believe the decay you are observing indicates that the sensor is stabilising gradually compared to other periods, as the humidity is not fluctuating much.

Regarding inputs and outputs — the inputs would be temperature and humidity, and the output would be the WE signal. However, for the model, we would need to teach it what we expect so it can learn the compensation relationship.

In the case of NO₂, for example, the inputs might be temperature, humidity, and WE, and the output could be a signal around 233 mV — not exactly flat but with minor variations. I might be mistaken, but perhaps expressing this as a percentage could help, so that when there are peaks caused by gas concentration, they are not interpreted as transients.

I believe you are from MATLAB Support? Do you have the option for an online call?

[Star Strider]

My pleasure!

At this point, I actually have something useful. You will need to have access to the System Identification Toolbox for this. I cannot simulate it any other way, even after having extracted the parameters (transfer function coefficients in this instance). The reason of course for using this approach is the dynamic nature of the signals and the system itself. I did all of it in MATLAB Online.

This is your code, with my addition being everything in the 'System Identification & Parameter Estimation' section. I initially began with a state=space identification, since that is usually preferable in my experience. It was not here. So I switched to transfer function identification, and got a decent result. The input to it is '[Temperature Humidity WE]' and the output was simply 'WE'. I then identified it with the appropriate parameter estimation to get the identified system. To get the Temperature-Humidity offset signal, I used the same identified system, however setting the 'WE' signal to a vector of zeros. I then subtracted that from the 'WE' signal in every situation to get the corrected 'WE' signal. Those results are plotted here.

I believe that I returned the correct signalss from 'makeFig' however I cannot guarantee that. I added outputs to itto get the results you originally plotted, then resampled them to a reasonable time interval, and then used that for the system identification dat. With rare exceptions, all such signal processing, system identification, and parameter estimation routines require regularly sampled signals. The relevant functions assume that the signaols are consistently samples, and will produce erroneous results if the signals are not consistently sampled. That is the reason for the table2timetable calls and the resulting restime calls.

The 'NumMtx' and 'DenMtx' are the numerator and denominator transfer function coefficients, in the event you want to experiment with them, perhaps using the transfer functions as filters for the 'WE' signal for each variable (Temperature and Humidity). I experimented with this, however I could not get that approach to work, the reason that I implement all of this in the System Identificatioon Toolbox. Tthe Control System Toolbox could also work.)

The 'NO' result is the best, so my approach works, at least for it. For the others, the result is less convincing, leading me to believe that in those, the Temperature and Humidity may be of lesser importance.

It is not easy for me to understand your code, however I did my best to understand it and provide a reasonable result.

Try this --

clear all
close all

matFile = websave('HumditydataS2.mat','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1842936/HumditydataS2.mat');

% Source Habs: https://www.reddit.com/r/embedded/comments/f4iino/whats_the_right_way_to_convert_relative_humidity/
Habs = @(RH,T) 216.7*(RH/100 .* 6.112 .* exp((17.62.*T)./(243.12+T)))./(273.15+T);  %Units: g/m^3
% HumAbs = Habs(T.('Humidity (%)'), T.('Temperature (°C)'));

% files = dir('*.mat');
% folder = files.folder;
% 
% %matFile = fullfile(folder, 'SensorLog_ALLS2.mat');
% matFile = fullfile(folder, 'HumditydataS2.mat');
temperature_test = false;   % set to true for the original dataset
% --- Load the combined table ---
S = load(matFile);
fn = fieldnames(S);
T  = S.(fn{find(structfun(@istable, S), 1)});
VN = T.Properties.VariableNames;

% --- Time & environment ---
if ~isdatetime(T.Timestamp)
    T.Timestamp = datetime(T.Timestamp,'InputFormat','yyyy-MM-dd HH:mm:ss');
end
time  = T.Timestamp;
TempC = T.("Temperature (°C)");
Hum   = T.("Humidity (%)");
HumAbs = Habs(T.('Humidity (%)'), T.('Temperature (°C)'));      % Calculate Approximate Absolute Humidity
% ================== CALIBRATION CONSTANTS ==================
% ---------- ELECTRONICS OFFSETS (hardware bias; set if measured) ----------
ELEC.NO_WE   = 0.292;  ELEC.NO_AUX   = 0.259;
ELEC.OX_WE   = 0.229;  ELEC.OX_AUX   = 0.227;   % OX uses O3 columns
ELEC.NO2_WE  = 0.238;  ELEC.NO2_AUX  = 0.225;
% ---------- SENSOR ZERO BASELINES (clean air, include electronics if measured that way) ----------
ZERO.NO_WE   = 0.3158;  ZERO.NO_AUX   = 0.3017;   % NO-B4
ZERO.OX_WE   = 0.2366;  ZERO.OX_AUX   = 0.2278;   % OX-A431 (O3)
ZERO.NO2_WE  = 0.2396;  ZERO.NO2_AUX  = 0.2253;   % NO2-B43F
% ===========================================================
% --- Ensure expected columns exist (adjust here if names differ) ---
need = ["NO WE","NO AUX","O3 WE","O3 AUX","NO2 WE","NO2 AUX"];
have = string(T.Properties.VariableNames);
if ~all(ismember(need, have))
    error('Missing columns. Expected: %s\nHave: %s', strjoin(need,', '), strjoin(have,', '));
end
% --- Initialise structs to avoid dot-indexing errors ---
WE = struct(); AE = struct(); WEo = struct(); AEo = struct();
% --- Convert zeros to sensor-only baselines (remove electronics bias) ---
WEo.NO  = ZERO.NO_WE  - ELEC.NO_WE;   AEo.NO  = ZERO.NO_AUX  - ELEC.NO_AUX;
WEo.OX  = ZERO.OX_WE  - ELEC.OX_WE;   AEo.OX  = ZERO.OX_AUX  - ELEC.OX_AUX;
WEo.NO2 = ZERO.NO2_WE - ELEC.NO2_WE;  AEo.NO2 = ZERO.NO2_AUX - ELEC.NO2_AUX;
% --- Raw -> electronics-corrected live readings ---
WE.NO   = T.("NO WE")   - ELEC.NO_WE;    AE.NO   = T.("NO AUX")  - ELEC.NO_AUX;
WE.OX   = T.("O3 WE")   - ELEC.OX_WE;    AE.OX   = T.("O3 AUX")  - ELEC.OX_AUX;
WE.NO2  = T.("NO2 WE")  - ELEC.NO2_WE;   AE.NO2  = T.("NO2 AUX") - ELEC.NO2_AUX;
% ---------------- Temperature compensation tables ----------------
T_pts   = [-30 -20 -10  0  10  20  30  40  50];
% NO (kT, ratio model)
kT_NO   = interp1(T_pts, [1.8 1.8 1.4 1.1 1.1 1.0 0.9 0.9 0.8], TempC, 'linear','extrap');
% NO2 & OX (nT, simple model) – replace OX row with your own if you have it
nT_NO2  = interp1(T_pts, [1.3 1.3 1.3 1.3 1.0 0.6 0.4 0.2 -1.5], TempC, 'linear','extrap');
nT_OX   = interp1(T_pts, [1.3 1.3 1.3 1.3 1.0 0.6 0.4 0.2 -1.5], TempC, 'linear','extrap');
% ---------------- Apply algorithms ----------------
% NO (Algorithm 2 / ratio): WEc = (WEu - WEe) - kT*(WEo/AEo)*(AEu - AEe)
WEc_NO   = (WE.NO) - kT_NO .* (ZERO.NO_WE ./ZERO.NO_AUX ) .* (AE.NO );
NO_NET = T.("NO WE") -T.("NO AUX") ; 
% NO2 (Algorithm 1 / nT):   WEc = (WEu - WEe) - nT*(AEu - AEe)
%WEc_NO2  = (WE.NO2 - WEo.NO2) - nT_NO2 .* (AE.NO2 - AEo.NO2);
WEc_NO2  = (WE.NO2) - nT_NO2 .* (AE.NO2);
NO2_NET = T.("NO2 WE") - T.("NO2 AUX");
% OX (Algorithm 1 / nT):    WEc = (WEu - WEe) - nT*(AEu - AEe)
WEc_OX   = (WE.OX ) - nT_OX  .* (AE.OX );
OX_NET = T.("O3 WE") - T.("O3 AUX");
% ---------------- Plot sets ----------------
[tq{1},tc_fq{1},NET_fq{1},WE_fq{1},Humq{1},nameq{1}] = makeFig("NO (k_T ratio)",  time, TempC, Hum, T.("NO WE") ,  T.("NO AUX"),  WEc_NO,NO_NET,temperature_test);

[tq{2},tc_fq{2},NET_fq{2},WE_fq{2},Humq{2},nameq{2}] = makeFig("NO2 (n_T simple)",time, TempC, Hum, T.("NO2 WE") , T.("NO2 AUX"), WEc_NO2,NO2_NET,temperature_test);

[tq{3},tc_fq{3},NET_fq{3},WE_fq{3},Humq{3},nameq{3}] = makeFig("OX/O3 (n_T simple)",time,TempC, Hum, T.("O3 WE") ,  T.("O3 AUX"), WEc_OX,OX_NET,temperature_test);

% ========== System Identification & Parameter Estimation ==========

fprintf('\nStart time:  %s\n\n',datetime('now'))

Start time:  17-Nov-2025 20:49:55

tic

SysOrdSS = [6; 6; 6];
SysOrdTF = [12; 14; 10];

for k = 1:numel(tq)
    fprintf('\n====================\nk = %d\nName: %s --\n\n',k,nameq{k})
    Tss = mean(diff(tq{k}));
    Tss.Format = "hh:mm:ss.SSS"
    Tsv(k) = floor(seconds(Tss));
    Fs = 1/Tsv(k)
    Tss = seconds(Tsv(k))
    DataTable{k} = table(tq{k},tc_fq{k},NET_fq{k},WE_fq{k},Humq{k}, VariableNames=["Timestamp","Temperature (°C)","NET","WE","Humidity (%)"]);
    HumAbs = Habs(Humq{k},tc_fq{k});
    Tr1 = DataTable{k};
    Tr1 = addvars(Tr1, HumAbs, NewVariableNames='Absolute Humidity');
    % TTr1 = table2timetable(Tr1(:,1:end-1));
    TTr1 = table2timetable(Tr1);
    fprintf('\nResampled Data --\n')
    TTr2 = retime(TTr1, 'regular', 'linear', TimeStep=Tss)
    
    y = TTr2.WE;
    % u = [TTr2.("Temperature (°C)") TTr2.("Absolute Humidity")];
    u = [TTr2.("Temperature (°C)"), TTr2.("Humidity (%)"), TTr2.WE];

    SIData = iddata(y,u,Tsv(k));
    % Sys{k} = ssest(SIData,SysOrdSS(k),Ts=Tsv(k));
    Sys{k} = tfest(SIData,SysOrdTF(k), Ts=Tsv(k));

    % fprintf('\n====================\nk = %d\nName: %s --\n\n',k,nameq{k})

    IDSystem = Sys{k}

    Parameters = IDSystem.Report.Parameters.ParVector;
    NPar = numel(Parameters)
    % PQ = [(1:NPar).' Parameters]
    % Szu = size(u)
    NPmtx = reshape(Parameters, [], size(u,2));
    % SzNP = size(NPmtx)
    % NPmtx2 = [(1:size(NPmtx,1)).' NPmtx]

    NumMtx = NPmtx(2:SysOrdTF(k),:)
    DenMtx = NPmtx(SysOrdTF(k)+1:size(NPmtx,1)-1,:)

    figure
    compare(SIData, Sys{k})
    hc = get(gca);
    xv = hc.Children(1).Children.XData;

    % for k2 = 1:2
    %     ofstvct = [-1 ones(1, size(DenMtx,1)-1)];
    %     partsys = tf(NumMtx(:,k2).', [1 (ofstvct .* DenMtx(:,k2).')], Ts=Tsv(k))
    %     out(:,k2) = lsim(partsys, u(:,k2).', xv(:).');
    % end

% QQ = out

    u2 = u;
    u2(:,3) = zeros(size(u,1),1);
    out = lsim(Sys{k}, u2, xv);

    figure
    plot(xv, out, LineWidth=2)
    grid
    ylim([-1 1]*0.05)
    xlabel('Time')
    ylabel('Temperature (°C) + Humidity (%) Offset Signal')
    title(nameq{k})

    figure
    plot(xv, TTr2.WE-out)
    grid
    xlabel('Time')
    ylabel('WE Signal - (Temperature (°C) + Humidity (%) Offset Signal)')
    title(nameq{k})
    % return

end

====================
k = 1
Name: NO (k_T ratio) --

Tss = duration
   00:00:05.685

Fs = 0.2000

Tss = duration
   5 sec

Resampled Data --

TTr2 = 3223×5 timetable
         Timestamp         Temperature (°C)      NET        WE      Humidity (%)    Absolute Humidity
    ___________________    ________________    _______    ______    ____________    _________________

    2025-11-07 11:24:20         20.214         -0.0094    0.2882       74.248            12.963      
    2025-11-07 11:24:25         20.104         -0.0084    0.2852       74.128             12.86      
    2025-11-07 11:24:30         19.994         -0.0086     0.281       74.026             12.76      
    2025-11-07 11:24:35         19.896         -0.0078    0.2766       73.972            12.678      
    2025-11-07 11:24:40           19.8         -0.0076     0.272       73.942            12.601      
    2025-11-07 11:24:45         19.652         -0.0086    0.2682       73.954            12.495      
    2025-11-07 11:24:50          19.54          -0.009     0.266           74            12.421      
    2025-11-07 11:24:55         19.472          -0.009    0.2646       74.054             12.38      
    2025-11-07 11:25:00         19.402         -0.0098    0.2632       74.126            12.341      
    2025-11-07 11:25:05         19.336         -0.0114    0.2608       74.212            12.308      
    2025-11-07 11:25:10         19.276          -0.012     0.259       74.296            12.278      
    2025-11-07 11:25:15          19.22         -0.0144     0.257       74.382            12.252      
    2025-11-07 11:25:20         19.174         -0.0176    0.2532       74.466            12.233      
    2025-11-07 11:25:25         19.134         -0.0166     0.252       74.556            12.219      
    2025-11-07 11:25:30         19.094          -0.016    0.2516        74.65            12.205      
    2025-11-07 11:25:35         19.058         -0.0164     0.251        74.75            12.196      

IDSystem =
 
  From input "u1" to output "y1":
                                                                                                                                                             
   0.0108 z^-1 - 0.01825 z^-2 + 0.006978 z^-3 + 0.01401 z^-4 - 0.01089 z^-5 - 0.002783 z^-6 - 0.009682 z^-7 + 0.01425 z^-8 + 0.0008552 z^-9 - 0.02089 z^-10  
                                                                                                                                                             
                                                                                                                                               + 0.0156 z^-11
                                                                                                                                                             
  -------------------------------------------------------------------------------------------------------------------------------------------------------------
                                                                                                                                                               
  1 + 0.082 z^-1 - 0.3421 z^-2 + 0.3774 z^-3 + 0.1827 z^-4 + 0.09624 z^-5 - 0.7252 z^-6 - 0.253 z^-7 + 0.4788 z^-8 - 0.4241 z^-9 - 0.1422 z^-10 - 0.1586 z^-11 
                                                                                                                                                               
                                                                                                                                              - 0.1139 z^-12   
                                                                                                                                                               
 
  From input "u2" to output "y1":
                                                                                                                                                             
   -9.514e-05 z^-1 + 0.0002647 z^-2 - 0.0001822 z^-3 - 0.0001455 z^-4 + 0.0002082 z^-5 - 5.979e-06 z^-6 - 3.486e-05 z^-7 + 7.935e-05 z^-8 - 0.0001754 z^-9   
                                                                                                                                                             
                                                                                                                          + 7.717e-05 z^-10 + 9.644e-06 z^-11
                                                                                                                                                             
  -------------------------------------------------------------------------------------------------------------------------------------------------------------
                                                                                                                                                               
  1 - 0.6568 z^-1 - 0.577 z^-2 + 0.3794 z^-3 - 0.4464 z^-4 - 0.2022 z^-5 - 0.0663 z^-6 + 0.1396 z^-7 + 0.5433 z^-8 - 0.3393 z^-9 + 0.4507 z^-10 + 0.6686 z^-11 
                                                                                                                                                               
                                                                                                                                              - 0.8935 z^-12   
                                                                                                                                                               
 
  From input "u3" to output "y1":
                                                                                                                                                           
   1.119 z^-1 - 0.5246 z^-2 + 0.5787 z^-3 - 0.307 z^-4 + 0.7303 z^-5 - 0.4016 z^-6 + 0.6047 z^-7 - 0.4135 z^-8 + 0.8077 z^-9 - 0.5584 z^-10 + 0.04216 z^-11
                                                                                                                                                           
  ----------------------------------------------------------------------------------------------------------------------------------------------------------
                                                                                                                                                            
  1 - 0.2666 z^-1 + 0.3026 z^-2 - 0.09455 z^-3 + 0.5163 z^-4 - 0.1507 z^-5 + 0.3424 z^-6 - 0.1697 z^-7 + 0.5862 z^-8 - 0.2027 z^-9 - 0.05469 z^-10          
                                                                                                                                                            
                                                                                                                              + 0.08059 z^-11 - 0.1727 z^-12
                                                                                                                                                            
 
Sample time: 5 seconds
Discrete-time identified transfer function.

Parameterization:
   Number of poles: [12 12 12]   Number of zeros: [11 11 11]
   Number of free coefficients: 69
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                            
Estimated using TFEST on time domain data "SIData".
Fit to estimation data: 91.21%                     
FPE: 4.954e-07, MSE: 4.641e-07                     
 

NPar = 75

NumMtx = 11×3

0.0108 -0.0001 1.1192 -0.0183 0.0003 -0.5246 0.0070 -0.0002 0.5787 0.0140 -0.0001 -0.3070 -0.0109 0.0002 0.7303 -0.0028 -0.0000 -0.4016 -0.0097 -0.0000 0.6047 0.0143 0.0001 -0.4135 0.0009 -0.0002 0.8077 -0.0209 0.0001 -0.5584 0.0156 0.0000 0.0422

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

●

DenMtx = 12×3

0.0820 -0.6568 -0.2666 -0.3421 -0.5770 0.3026 0.3774 0.3794 -0.0946 0.1827 -0.4464 0.5163 0.0962 -0.2022 -0.1507 -0.7252 -0.0663 0.3424 -0.2530 0.1396 -0.1697 0.4788 0.5433 0.5862 -0.4241 -0.3393 -0.2027 -0.1422 0.4507 -0.0547 -0.1586 0.6686 0.0806 -0.1139 -0.8935 -0.1727

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

====================
k = 2
Name: NO2 (n_T simple) --

Tss = duration
   00:00:05.685

Fs = 0.2000

Tss = duration
   5 sec

Resampled Data --

TTr2 = 3223×5 timetable
         Timestamp         Temperature (°C)      NET        WE      Humidity (%)    Absolute Humidity
    ___________________    ________________    _______    ______    ____________    _________________

    2025-11-07 11:24:20         20.214         -0.0254    0.2086       74.248            12.963      
    2025-11-07 11:24:25         20.104         -0.0244    0.2096       74.128             12.86      
    2025-11-07 11:24:30         19.994         -0.0252      0.21       74.026             12.76      
    2025-11-07 11:24:35         19.896         -0.0254    0.2106       73.972            12.678      
    2025-11-07 11:24:40           19.8         -0.0244    0.2122       73.942            12.601      
    2025-11-07 11:24:45         19.652         -0.0228    0.2148       73.954            12.495      
    2025-11-07 11:24:50          19.54          -0.021     0.217           74            12.421      
    2025-11-07 11:24:55         19.472           -0.02     0.218       74.054             12.38      
    2025-11-07 11:25:00         19.402           -0.02    0.2184       74.126            12.341      
    2025-11-07 11:25:05         19.336         -0.0196    0.2198       74.212            12.308      
    2025-11-07 11:25:10         19.276         -0.0186    0.2214       74.296            12.278      
    2025-11-07 11:25:15          19.22         -0.0172    0.2228       74.382            12.252      
    2025-11-07 11:25:20         19.174         -0.0164    0.2244       74.466            12.233      
    2025-11-07 11:25:25         19.134         -0.0166    0.2254       74.556            12.219      
    2025-11-07 11:25:30         19.094         -0.0152    0.2264        74.65            12.205      
    2025-11-07 11:25:35         19.058         -0.0136    0.2278        74.75            12.196      

IDSystem =
 
  From input "u1" to output "y1":
                                                                                                                                                             
   -0.0005966 z^-1 + 0.002684 z^-2 - 0.002937 z^-3 + 0.001092 z^-4 - 0.001974 z^-5 + 0.002478 z^-6 - 0.0009937 z^-7 + 0.001086 z^-8 - 0.002645 z^-9          
                                                                                                                                                             
                                                                                       + 0.002546 z^-10 - 0.0008166 z^-11 + 0.0007472 z^-12 - 0.0006697 z^-13
                                                                                                                                                             
  -------------------------------------------------------------------------------------------------------------------------------------------------------------
                                                                                                                                                               
  1 - 1.444 z^-1 + 0.05416 z^-2 + 0.7511 z^-3 - 0.1106 z^-4 - 0.5981 z^-5 + 0.4543 z^-6 - 0.6751 z^-7 + 1.145 z^-8 - 0.2812 z^-9 - 0.8826 z^-10 + 0.7122 z^-11 
                                                                                                                                                               
                                                                                                                + 0.2981 z^-12 - 0.6721 z^-13 + 0.2593 z^-14   
                                                                                                                                                               
 
  From input "u2" to output "y1":
                                                                                                                                                            
  0.001145 z^-1 - 0.001986 z^-2 + 0.0002689 z^-3 + 0.001286 z^-4 - 0.001076 z^-5 + 0.0004332 z^-6 - 0.0002389 z^-7 - 9.016e-05 z^-8 + 0.0007838 z^          
                                                                                                                                                            
                                                                                    -9 - 0.001072 z^-10 + 2.319e-05 z^-11 + 0.001357 z^-12 - 0.0008357 z^-13
                                                                                                                                                            
  ----------------------------------------------------------------------------------------------------------------------------------------------------------
                                                                                                                                                            
  1 - 0.7576 z^-1 - 0.5825 z^-2 + 0.8611 z^-3 - 0.2566 z^-4 - 0.06305 z^-5 - 0.1491 z^-6 - 0.2316 z^-7 + 0.4486 z^-8 - 0.3881 z^-9 - 0.4974 z^-10           
                                                                                                                                                            
                                                                                               + 0.7163 z^-11 - 0.002167 z^-12 - 0.1959 z^-13 + 0.1362 z^-14
                                                                                                                                                            
 
  From input "u3" to output "y1":
                                                                                                                                                            
  1.179 z^-1 + 0.0958 z^-2 - 0.003798 z^-3 - 0.1792 z^-4 + 0.4969 z^-5 + 0.1418 z^-6 - 0.03119 z^-7 + 0.4185 z^-8 - 0.264 z^-9 + 0.162 z^-10 - 0.1246 z^-11 
                                                                                                                                                            
                                                                                                                               + 0.159 z^-12 + 0.01494 z^-13
                                                                                                                                                            
  ----------------------------------------------------------------------------------------------------------------------------------------------------------
                                                                                                                                                            
  1 + 0.2766 z^-1 - 0.02588 z^-2 - 0.1642 z^-3 + 0.3952 z^-4 + 0.2188 z^-5 + 0.007276 z^-6 + 0.3659 z^-7 - 0.1889 z^-8 + 0.1354 z^-9 - 0.02257 z^-10        
                                                                                                                                                            
                                                                                               + 0.1082 z^-11 - 0.05823 z^-12 - 0.00718 z^-13 + 0.0211 z^-14
                                                                                                                                                            
 
Sample time: 5 seconds
Discrete-time identified transfer function.

Parameterization:
   Number of poles: [14 14 14]   Number of zeros: [13 13 13]
   Number of free coefficients: 81
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                            
Estimated using TFEST on time domain data "SIData".
Fit to estimation data: 27.21%                     
FPE: 2.117e-05, MSE: 2.013e-05                     
 

NPar = 87

NumMtx = 13×3

-0.0006 0.0011 1.1786 0.0027 -0.0020 0.0958 -0.0029 0.0003 -0.0038 0.0011 0.0013 -0.1792 -0.0020 -0.0011 0.4969 0.0025 0.0004 0.1418 -0.0010 -0.0002 -0.0312 0.0011 -0.0001 0.4185 -0.0026 0.0008 -0.2640 0.0025 -0.0011 0.1620 -0.0008 0.0000 -0.1246 0.0007 0.0014 0.1590 -0.0007 -0.0008 0.0149

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

DenMtx = 14×3

-1.4443 -0.7576 0.2766 0.0542 -0.5825 -0.0259 0.7511 0.8611 -0.1642 -0.1106 -0.2566 0.3952 -0.5981 -0.0631 0.2188 0.4543 -0.1491 0.0073 -0.6751 -0.2316 0.3659 1.1447 0.4486 -0.1889 -0.2812 -0.3881 0.1354 -0.8826 -0.4974 -0.0226 0.7122 0.7163 0.1082 0.2981 -0.0022 -0.0582 -0.6721 -0.1959 -0.0072 0.2593 0.1362 0.0211

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

====================
k = 3
Name: OX/O3 (n_T simple) --

Tss = duration
   00:00:05.685

Fs = 0.2000

Tss = duration
   5 sec

Resampled Data --

TTr2 = 3223×5 timetable
         Timestamp         Temperature (°C)      NET        WE      Humidity (%)    Absolute Humidity
    ___________________    ________________    _______    ______    ____________    _________________

    2025-11-07 11:24:20         20.214         -0.0244     0.207       74.248            12.963      
    2025-11-07 11:24:25         20.104         -0.0234     0.207       74.128             12.86      
    2025-11-07 11:24:30         19.994          -0.023    0.2064       74.026             12.76      
    2025-11-07 11:24:35         19.896         -0.0224    0.2066       73.972            12.678      
    2025-11-07 11:24:40           19.8          -0.022     0.207       73.942            12.601      
    2025-11-07 11:24:45         19.652         -0.0214    0.2076       73.954            12.495      
    2025-11-07 11:24:50          19.54           -0.02    0.2085           74            12.421      
    2025-11-07 11:24:55         19.472         -0.0186    0.2094       74.054             12.38      
    2025-11-07 11:25:00         19.402         -0.0176    0.2104       74.126            12.341      
    2025-11-07 11:25:05         19.336          -0.017    0.2114       74.212            12.308      
    2025-11-07 11:25:10         19.276          -0.017     0.212       74.296            12.278      
    2025-11-07 11:25:15          19.22         -0.0166    0.2124       74.382            12.252      
    2025-11-07 11:25:20         19.174         -0.0156    0.2134       74.466            12.233      
    2025-11-07 11:25:25         19.134          -0.015     0.214       74.556            12.219      
    2025-11-07 11:25:30         19.094         -0.0146    0.2144        74.65            12.205      
    2025-11-07 11:25:35         19.058          -0.014    0.2154        74.75            12.196      

IDSystem =
 
  From input "u1" to output "y1":
  -0.0001083 z^-1 + 8.177e-05 z^-2 - 8.476e-05 z^-3 + 0.0001396 z^-4 - 5.401e-05 z^-5 + 4.415e-05 z^-6 - 3.235e-05 z^-7 - 1.972e-05 z^-8 + 3.465e-05 z^-9
  -------------------------------------------------------------------------------------------------------------------------------------------------------
    1 - 0.7071 z^-1 - 0.3312 z^-2 - 0.001249 z^-3 - 0.04412 z^-4 + 0.4732 z^-5 - 0.04704 z^-6 + 0.3481 z^-7 - 0.7103 z^-8 - 0.3272 z^-9 + 0.5477 z^-10
 
  From input "u2" to output "y1":
  0.0007567 z^-1 - 0.001473 z^-2 + 0.001213 z^-3 - 0.00145 z^-4 + 0.001413 z^-5 - 0.0008355 z^-6 + 0.0008114 z^-7 - 0.0005198 z^-8 + 8.41e-05 z^-9
  ------------------------------------------------------------------------------------------------------------------------------------------------
  1 - 0.7872 z^-1 + 0.1703 z^-2 - 0.6647 z^-3 + 0.2157 z^-4 + 0.1862 z^-5 + 0.1994 z^-6 - 0.1049 z^-7 - 0.1325 z^-8 - 0.05429 z^-9 + 0.09688 z^-10
 
  From input "u3" to output "y1":
             1.314 z^-1 - 2.064 z^-2 + 0.5294 z^-3 + 0.4738 z^-4 + 0.8795 z^-5 - 1.669 z^-6 + 1.05 z^-7 - 0.3542 z^-8 + 0.01494 z^-9
  ---------------------------------------------------------------------------------------------------------------------------------------------
  1 - 1.309 z^-1 - 0.03812 z^-2 + 0.5617 z^-3 + 0.725 z^-4 - 1.177 z^-5 + 0.485 z^-6 + 0.0509 z^-7 - 0.1108 z^-8 - 0.03958 z^-9 + 0.02717 z^-10
 
Sample time: 5 seconds
Discrete-time identified transfer function.

Parameterization:
   Number of poles: [10 10 10]   Number of zeros: [9 9 9]
   Number of free coefficients: 57
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                            
Estimated using TFEST on time domain data "SIData".
Fit to estimation data: 38.14%                     
FPE: 1.936e-05, MSE: 1.869e-05                     
 

NPar = 63

NumMtx = 9×3

-0.0001 0.0008 1.3145 0.0001 -0.0015 -2.0643 -0.0001 0.0012 0.5294 0.0001 -0.0014 0.4738 -0.0001 0.0014 0.8795 0.0000 -0.0008 -1.6685 -0.0000 0.0008 1.0499 -0.0000 -0.0005 -0.3542 0.0000 0.0001 0.0149

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

DenMtx = 10×3

-0.7071 -0.7872 -1.3094 -0.3312 0.1703 -0.0381 -0.0012 -0.6647 0.5617 -0.0441 0.2157 0.7250 0.4732 0.1862 -1.1772 -0.0470 0.1994 0.4850 0.3481 -0.1049 0.0509 -0.7103 -0.1325 -0.1108 -0.3272 -0.0543 -0.0396 0.5477 0.0969 0.0272

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fprintf('\n\nFinish time:  %s\n\n',datetime('now'))

Finish time:  17-Nov-2025 20:50:16

toc

Elapsed time is 21.123490 seconds.

fprintf('\n\nFinished!\n==========\n')

Finished!
==========

%% --------------- Plot helper ---------------
function [t,tc_f,NET_f,WE_f,hum,name] = makeFig(name,t,tc,hum,WE,AUX,WEc,NET,temperature_test)
   % if ~exist('temperature_test','var')
   %              temperature_test = false;   % default
   % end
    f = figure('Name',name,'Color','w');
    tl = tiledlayout(f,2,1,'TileSpacing','compact','Padding','compact');
    % Top: Temperature & Humidity
    ax1 = nexttile(tl);
    yyaxis left,  plot(t,tc,'LineWidth',1.2), ylabel('Temperature (°C)')
    yyaxis right, plot(t,hum,'LineWidth',1.2), ylabel('Humidity (%)')
    grid on, title([name,' — Temp & Humidity'])
    % Bottom: WE/AUX vs Corrected
    ax2 = nexttile(tl);
    yyaxis left
    if temperature_test
    plot(t,WE,'LineWidth',1.2); hold on
    plot(t,AUX,'LineWidth',1.2); ylabel('Raw Sensor (V)')
    else
     plot(t, WE, 'b', 'LineWidth', 1.2); hold on      % blue line for WE
     plot(t, AUX, 'r', 'LineWidth', 1.2);             % red line for AUX
     ylabel('Raw Sensor (V)');
    end;
    
    if temperature_test
        yyaxis right
        plot(t,WEc,'LineWidth',1.5); ylabel('Corrected WE_c (V)'); hold on
        %plot(t,NET,'LineWidth',1.5); 
        plot(t, NET, 'LineWidth', 1.5, 'Color', 'g', 'LineStyle', '-');
    end
    xlabel('Time'), grid on
    title([name,' — WE/AUX (left), WE_c (right)'])
    
    if temperature_test
        legend({'WE','AUX','WE Corrected*(Alpha Model)','Net Corrected'},'Location','best')
    else
        legend({'WE','AUX'},'Location','best')
    end
    linkaxes([ax1 ax2],'x');
    %f2 = figure('Name',[name ' — NET vs Temperature'],'Color','w');
            % Use fixed window for the original temperature test; otherwise use full span
         
               
            if temperature_test 
                % --- Fixed window for original dataset ---
                t0 = datetime(2025,10,30,11,41,0);
                t1 = datetime(2025,10,30,15,05,0);
            else
                % --- Full-range window for current data (with a small padding) ---
                tt = t(~isnat(t));                 % remove NaT, if any
                if isempty(tt)
                    warning('No valid timestamps in t; skipping.');
                    return
                end
                t0 = min(tt) - minutes(1);
                t1 = max(tt) + minutes(1);
            end
            
            % --- Build one mask only ---
            inwin = (t >= t0) & (t <= t1);
            good  = isfinite(tc) & isfinite(NET) & isfinite(WE);   % include WE since you use WE_f too
            mask  = inwin & good;
            
            % --- Filtered vectors ---
            tc_f  = tc(mask);
            NET_f = NET(mask);
            WE_f  = WE(mask);
            
        if temperature_test 
                % --- Plot scatter ---
                figure('Color','w');
                %scatter(tc_f, NET_f, 5, 'filled'); grid on; hold on
                scatter(tc_f, WE_f, 5, 'filled'); grid on; hold on
                xlabel('Temperature (°C)');
                ylabel('NET (V)');
                title("NET vs Temperature — " + string(name));
                
                % --- Polynomial fit (choose degree) ---
                deg = 3;                                        % << change degree here
                deg = min([deg, numel(tc_f)-1, 6]);             % guard: enough points & avoid wild fits
                
                % Fit and trend line
                p    = polyfit(tc_f(:), WE_f(:), deg);         % coefficients, highest power first
                xfit = linspace(min(tc_f), max(tc_f), 400);
                yfit = polyval(p, xfit);
                plot(xfit, yfit, 'LineWidth', 1.6);
                
                % R^2
                yhat  = polyval(p, tc_f(:));
                SSres = sum((WE_f(:) - yhat).^2);
                SStot = sum((WE_f(:) - mean(WE_f(:))).^2);
                R2 = 1 - SSres/SStot;
                
                % Legend: show full equation for deg <= 3, otherwise compact label
                if deg == 1
                    eqStr = sprintf('y = %.3g x %+ .3g  (R^2 = %.3f)', p(1), p(2), R2);
                elseif deg == 2
                    a=p(1); b=p(2); c=p(3);
                    eqStr = sprintf('y = %.3g x^2 %+ .3g x %+ .3g  (R^2 = %.3f)', a,b,c,R2);
                elseif deg == 3
                    a=p(1); b=p(2); c=p(3); d=p(4);
                    eqStr = sprintf('y = %.3g x^3 %+ .3g x^2 %+ .3g x %+ .3g  (R^2 = %.3f)', a,b,c,d,R2);
                else
                    eqStr = sprintf('Poly%d fit (R^2 = %.3f)', deg, R2);
                end
                legend('Data', eqStr, 'Location','best');
        end
end




% =========================================================================

This is the best I can do with your data to provide the result you want.

.