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.