Add data2 impedance analysis

data2/data2_analysis.m

@@ -0,0 +1,219 @@
+clear; clc; close all;
+
+scriptDir = fileparts(mfilename('fullpath'));
+if isempty(scriptDir)
+    scriptDir = pwd;
+end
+
+load(fullfile(scriptDir, 'pp.mat'));
+
+figDir = fullfile(scriptDir, 'figures');
+resultDir = fullfile(scriptDir, 'results');
+
+if ~exist(figDir, 'dir')
+    mkdir(figDir);
+end
+if ~exist(resultDir, 'dir')
+    mkdir(resultDir);
+end
+
+sweep = (1:numel(R_peak))';
+
+features = [R_peak(:), X1_peak(:), X2_peak(:), Z_abs_peak(:), ...
+            Phase_peak(:), Y_abs_peak(:), G_peak(:), B_peak(:)];
+
+featureNames = {'R_peak','X1_peak','X2_peak','Z_abs_peak', ...
+                'Phase_peak','Y_abs_peak','G_peak','B_peak'};
+
+% Exclude pump degradation after sweep 1785
+valid = sweep < 1785;
+
+featuresValid = features(valid, :);
+sweepValid = sweep(valid);
+
+% Segment valid sweeps into 21 concentration plateaus: 0.0 to 2.0%
+nLevels = 21;
+concLevels = (0:0.1:2.0)';
+edges = round(linspace(min(sweepValid), max(sweepValid) + 1, nLevels + 1));
+
+plateauID = nan(size(sweep));
+concSweep = nan(size(sweep));
+
+for k = 1:nLevels
+    idx = sweep >= edges(k) & sweep < edges(k+1) & valid;
+    plateauID(idx) = k;
+    concSweep(idx) = concLevels(k);
+end
+
+% Figure 1: time/sweep series
+figure('Color','w');
+plot(sweep, R_peak, 'LineWidth', 1.2);
+hold on;
+xline(1785, '--r', 'Pump failure onset');
+xlabel('Sweep number');
+ylabel('R peak frequency (Hz)');
+title('R peak resonance frequency over experiment');
+grid on;
+saveas(gcf, fullfile(figDir, 'r_peak_timeseries.png'));
+
+% Plateau statistics
+meanVals = nan(nLevels, numel(featureNames));
+stdVals = nan(nLevels, numel(featureNames));
+
+for k = 1:nLevels
+    idx = plateauID == k;
+    meanVals(k, :) = mean(features(idx, :), 1, 'omitnan');
+    stdVals(k, :) = std(features(idx, :), 0, 1, 'omitnan');
+end
+
+meanTable = array2table(meanVals, 'VariableNames', featureNames);
+stdTable = array2table(stdVals, 'VariableNames', strcat(featureNames, '_std'));
+
+plateauStats = [table(concLevels, 'VariableNames', {'Concentration_percent'}), ...
+                meanTable, stdTable];
+
+writetable(plateauStats, fullfile(resultDir, 'plateau_statistics.csv'));
+
+% Calibration analysis for three representative features
+selectedIdx = [1 2 7]; % R_peak, X1_peak, G_peak
+calibResults = table();
+
+figure('Color','w');
+tiledlayout(1, 3, 'TileSpacing', 'compact');
+
+for j = 1:numel(selectedIdx)
+    col = selectedIdx(j);
+
+    x = concLevels;
+    y = meanVals(:, col);
+    yerr = stdVals(:, col);
+
+    [intercept, sensitivity, R2, xfit, yfit, yci] = simpleLinearFit(x, y);
+
+    noiseFloor = yerr(1);
+    LOD = 3 * noiseFloor / abs(sensitivity);
+
+    calibResults.Feature(j, 1) = string(featureNames{col});
+    calibResults.Sensitivity_Hz_per_percent(j, 1) = sensitivity;
+    calibResults.Intercept_Hz(j, 1) = intercept;
+    calibResults.R2(j, 1) = R2;
+    calibResults.NoiseFloor_Hz(j, 1) = noiseFloor;
+    calibResults.LOD_percent(j, 1) = LOD;
+
+    nexttile;
+    errorbar(x, y, yerr, 'o', 'LineWidth', 1.1);
+    hold on;
+
+    plot(xfit, yfit, 'r-', 'LineWidth', 1.5);
+    plot(xfit, yci, 'r--', 'LineWidth', 0.8);
+
+    xlabel('Glycerol concentration (% w/v)');
+    ylabel('Frequency (Hz)');
+    title(strrep(featureNames{col}, '_', '\_'));
+    grid on;
+end
+
+saveas(gcf, fullfile(figDir, 'calibration_curves.png'));
+writetable(calibResults, fullfile(resultDir, 'calibration_results.csv'));
+
+% PCA on plateau mean feature matrix using SVD, avoiding toolbox dependency
+X = standardizeColumns(meanVals);
+[U, S, coeff] = svd(X, 'econ');
+score = U * S;
+latent = diag(S).^2 ./ (size(X, 1) - 1);
+explained = 100 * latent ./ sum(latent);
+
+pcaTable = table((1:numel(explained))', explained, cumsum(explained), ...
+    'VariableNames', {'PC', 'Explained_percent', 'Cumulative_percent'});
+writetable(pcaTable, fullfile(resultDir, 'pca_explained_variance.csv'));
+
+figure('Color','w');
+scatter(score(:,1), score(:,2), 50, concLevels, 'filled');
+xlabel(sprintf('PC1 (%.1f%%)', explained(1)));
+ylabel(sprintf('PC2 (%.1f%%)', explained(2)));
+title('PCA of standardized impedance features');
+cb = colorbar;
+cb.Label.String = 'Glycerol concentration (% w/v)';
+grid on;
+saveas(gcf, fullfile(figDir, 'pca_result.png'));
+
+% Simple multivariate regression using first two PCs
+Xreg = score(:, 1:2);
+yreg = concLevels;
+
+[pcrIntercept, pcrBeta, pcrR2, yhat] = simpleMultipleFit(Xreg, yreg);
+rmsePCR = sqrt(mean((yreg - yhat).^2));
+
+figure('Color','w');
+plot(yreg, yhat, 'o', 'LineWidth', 1.2);
+hold on;
+plot([0 2], [0 2], 'k--');
+xlabel('Reference concentration (% w/v)');
+ylabel('Predicted concentration (% w/v)');
+title(sprintf('PCR prediction, RMSE = %.4f %%', rmsePCR));
+grid on;
+saveas(gcf, fullfile(figDir, 'pcr_prediction.png'));
+
+multivarResults = table(rmsePCR, pcrR2, ...
+    'VariableNames', {'PCR_RMSE_percent', 'PCR_R2'});
+writetable(multivarResults, fullfile(resultDir, 'multivariate_results.csv'));
+disp('Calibration results:');
+disp(calibResults);
+
+disp('PCA explained variance:');
+disp(pcaTable(1:3, :));
+
+disp('Done. Figures saved in /figures and CSV files saved in /results.');
+function [intercept, slope, R2, xfit, yfit, yci] = simpleLinearFit(x, y)
+    x = x(:);
+    y = y(:);
+    keep = isfinite(x) & isfinite(y);
+    x = x(keep);
+    y = y(keep);
+
+    A = [ones(size(x)), x];
+    beta = A \ y;
+    intercept = beta(1);
+    slope = beta(2);
+
+    yhat = A * beta;
+    residual = y - yhat;
+    sse = sum(residual.^2);
+    sst = sum((y - mean(y)).^2);
+    R2 = 1 - sse / sst;
+
+    xfit = linspace(min(x), max(x), 100)';
+    yfit = intercept + slope .* xfit;
+
+    n = numel(x);
+    s2 = sse / max(n - 2, 1);
+    xbar = mean(x);
+    sxx = sum((x - xbar).^2);
+    sePred = sqrt(s2 .* (1 + 1/n + ((xfit - xbar).^2 ./ sxx)));
+    yci = [yfit - 2 * sePred, yfit + 2 * sePred];
+end
+
+function Xs = standardizeColumns(X)
+    mu = mean(X, 1, 'omitnan');
+    sigma = std(X, 0, 1, 'omitnan');
+    sigma(sigma == 0 | ~isfinite(sigma)) = 1;
+    Xs = (X - mu) ./ sigma;
+    Xs(~isfinite(Xs)) = 0;
+end
+
+function [intercept, slopes, R2, yhat] = simpleMultipleFit(X, y)
+    y = y(:);
+    keep = all(isfinite(X), 2) & isfinite(y);
+    X = X(keep, :);
+    y = y(keep);
+
+    A = [ones(size(X, 1), 1), X];
+    beta = A \ y;
+    intercept = beta(1);
+    slopes = beta(2:end);
+    yhat = A * beta;
+
+    sse = sum((y - yhat).^2);
+    sst = sum((y - mean(y)).^2);
+    R2 = 1 - sse / sst;
+end

data2/results/calibration_results.csv

@@ -0,0 +1,4 @@
+Feature,Sensitivity_Hz_per_percent,Intercept_Hz,R2,NoiseFloor_Hz,LOD_percent
+R_peak,1534.01034070498,10009082.526501,0.977026816185807,157.792903229194,0.308588995215008
+X1_peak,2294.7818803054,10003961.9475113,0.997095488857794,292.088895834743,0.381851841791434
+G_peak,-4.0597168738683e-05,0.000781300190482523,0.328612302165231,0.000422239623051608,31.2021480440785

data2/results/multivariate_results.csv

@@ -0,0 +1,2 @@
+PCR_RMSE_percent,PCR_R2
+0.067618748601777,0.987530104102354

data2/results/pca_explained_variance.csv

@@ -0,0 +1,9 @@
+PC,Explained_percent,Cumulative_percent
+1,83.4224323054342,83.4224323054342
+2,15.1818480320896,98.6042803375238
+3,1.31783735889597,99.9221176964198
+4,0.0679975396508157,99.9901152360706
+5,0.00939210937645678,99.9995073454471
+6,0.00029271378379679,99.9998000592309
+7,0.000191112212537476,99.9999911714434
+8,8.82855657319047e-06,100

data2/results/plateau_statistics.csv

@@ -0,0 +1,22 @@
+Concentration_percent,R_peak,X1_peak,X2_peak,Z_abs_peak,Phase_peak,Y_abs_peak,G_peak,B_peak,R_peak_std,X1_peak_std,X2_peak_std,Z_abs_peak_std,Phase_peak_std,Y_abs_peak_std,G_peak_std,B_peak_std
+0,10008818.6202244,10003718.1729432,10013861.1501161,10013384.2476455,10008460.7835786,10003229.2556246,0.000736466399461239,10002806.103232,157.792903229194,292.088895834743,191.333780088756,211.412115683534,180.885982144457,351.764380149965,0.000422239623051608,418.265170070091
+0.1,10009120.4911768,10004243.8616391,10014192.9119009,10013751.3714685,10008799.745608,10003838.8003449,0.000743576080196564,10003498.9128254,48.9752860675288,80.889247359638,77.5664080432058,70.7560847655657,52.0693720271209,90.0046751774182,0.000385734141973868,114.312604993474
+0.2,10009303.2135538,10004490.2340794,10014402.6078434,10013967.8692993,10008992.4738339,10004116.1653913,0.000701618077325986,10003811.8079134,50.001056289808,80.0560516723227,80.4488521022321,72.2035759503678,53.923059953112,85.1220969865212,0.000401260598068668,98.2034284893044
+0.3,10009456.2930616,10004700.1027033,10014579.5135648,10014145.714203,10009151.8740741,10004331.510171,0.000739381973359067,10004054.4901578,45.1500543524943,70.9759918577537,82.0141830027889,76.3264874338652,46.6447924303461,69.8369619660059,0.000369611795661687,84.0295122966274
+0.4,10009619.7607434,10004872.4487497,10014838.3377494,10014367.9660531,10009318.6540919,10004525.7437826,0.00081443977138571,10004249.4119101,47.4910261824401,48.8461007122343,101.853236219637,102.45979040755,45.8074587833404,61.8671054547815,0.000297324013784061,67.1052272935271
+0.5,10009795.4986371,10005041.685585,10015113.5560904,10014652.4713274,10009494.8034464,10004689.6360326,0.000772525486270561,10004422.9404107,60.5057492727399,73.0393979988469,93.4990075515101,99.6021056416084,57.4223521371279,71.6704506141382,0.000343590032237065,76.8392567862517
+0.6,10010037.3133207,10005286.4845478,10015425.5005284,10014944.9325711,10009729.3560669,10004927.4483663,0.00078493007298431,10004647.7212692,86.121275623232,90.1356147280536,132.808038498789,125.567409809073,82.8803708187052,85.9299391959141,0.000333825706014917,89.7411960667606
+0.7,10010373.5770236,10005677.227306,10015661.0340511,10015165.2692231,10010056.1731665,10005311.3163384,0.000761066366297671,10005025.3147681,67.4515888982397,111.425175437931,84.5453420084274,54.4063166385441,63.8805420130217,110.240248540513,0.00033802390593995,112.062004506787
+0.8,10010503.7453836,10005882.6704507,10015599.5355723,10015128.1770821,10010184.9607717,10005523.8105012,0.000777836646959189,10005238.7588137,36.478289889637,61.1982716623814,56.2886987534757,49.3877100879082,34.4495087917822,55.2936796466378,0.0003000329644916,70.9250964555539
+0.9,10010627.6206675,10006064.7414002,10015568.8197771,10015113.7462554,10010307.935871,10005699.924013,0.000805370243780016,10005414.01079,45.2008709038458,73.0910229777243,54.1200006971917,52.9248254503344,45.7855392187891,74.0149233375511,0.00024182766394971,79.7159195595029
+1,10010823.3023993,10006325.0689245,10015559.2762126,10015156.7199325,10010497.155248,10005959.0668567,0.000757683518316722,10005675.7904159,45.8062657405717,69.2131270854152,61.0833489774642,45.4232013241563,44.8738523310153,74.0990894561798,0.000278846112070904,74.7614294421997
+1.1,10010941.5443888,10006525.8927346,10015535.0287654,10015130.9056204,10010613.7615047,10006154.4197024,0.000780410687902343,10005870.856033,37.809242185673,58.7484180373748,41.4028841738087,41.1020496137892,36.4680941081756,64.5533814758243,0.000215953478524817,58.0907808475414
+1.2,10011074.6059115,10006749.3686383,10015530.0304987,10015152.5172771,10010745.2107955,10006377.439563,0.000763985048223773,10006095.9646883,35.0757602380582,70.5972718740799,42.0551996288374,37.9801503510197,33.8375109408663,71.1557398003329,0.000210712460145705,68.8472853349688
+1.3,10011175.6595051,10006946.1403817,10015501.6070813,10015130.0660806,10010845.3599646,10006569.8034045,0.000735289908029619,10006293.6527588,24.3306317407354,50.0753130627641,39.7288337152158,39.780384827109,24.0920732205884,48.8003530948486,0.000220296758752495,49.3382284014013
+1.4,10011271.7995663,10007134.4507057,10015502.2390649,10015131.8768145,10010942.6904714,10006768.9793066,0.000728909554583364,10006487.6016378,27.0705059850872,61.5143702604628,38.6578022379767,38.7559517868194,27.851760719384,62.6363248111383,0.000201138470084967,68.7364795833065
+1.5,10011393.8435164,10007384.0180973,10015495.4649449,10015127.8972814,10011065.5951451,10007022.804384,0.000743737451490402,10006745.9748531,42.3682198017505,93.089553211294,42.9203421606586,36.1937196648681,42.4667607356324,90.6326169869339,0.000115603166475932,94.0751904968302
+1.6,10011514.3411225,10007624.869592,10015518.6295024,10015152.2420455,10011187.6529169,10007261.0892849,0.000706071826031988,10006988.6206598,28.7619137472294,64.078506683881,46.1235155526921,34.2003930898357,28.8693743548172,62.7953637354422,0.000157295568981629,66.8491650442687
+1.7,10011628.8744335,10007860.7178777,10015544.6399884,10015171.347603,10011304.162364,10007498.5696673,0.000704052174978828,10007235.0342037,35.2654539076451,78.6740102489567,42.76429706843,32.9140270178755,35.0201872907768,78.0942215700027,0.000109656495621124,80.5216050560525
+1.8,10011717.703354,10008080.7768344,10015506.6305957,10015125.3286157,10011393.0016038,10007718.7099455,0.000682369101861647,10007460.7702825,28.8046958168206,64.1988914115464,42.2943784247508,39.2917644444486,29.1186035311627,70.8061936746101,0.000106156210350962,70.106414842782
+1.9,10011788.269672,10008217.1388071,10015503.9638347,10015133.1671327,10011463.7667661,10007852.4732109,0.000670041508157615,10007600.3835029,21.0612034230767,47.7181498984466,42.6045715279601,35.1710378326526,21.9217002914185,47.4669925819639,0.00010422870663713,51.587308616895
+2,10011961.1960142,10008565.2452259,10015524.9740088,10015150.011244,10011653.2409125,10008243.0212937,0.000645001559024031,10007959.0499885,465.228120501886,785.338026599368,269.765327389109,199.17128420168,478.975068711591,841.756124893689,9.43522383201465e-05,727.772917587464