diff --git a/data2/data2_analysis.m b/data2/data2_analysis.m new file mode 100644 index 0000000..ea0063c --- /dev/null +++ b/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 diff --git a/data2/figures/calibration_curves.png b/data2/figures/calibration_curves.png new file mode 100644 index 0000000..8bf8c6e Binary files /dev/null and b/data2/figures/calibration_curves.png differ diff --git a/data2/figures/pca_result.png b/data2/figures/pca_result.png new file mode 100644 index 0000000..24b9b42 Binary files /dev/null and b/data2/figures/pca_result.png differ diff --git a/data2/figures/pcr_prediction.png b/data2/figures/pcr_prediction.png new file mode 100644 index 0000000..f5fb89f Binary files /dev/null and b/data2/figures/pcr_prediction.png differ diff --git a/data2/figures/r_peak_timeseries.png b/data2/figures/r_peak_timeseries.png new file mode 100644 index 0000000..72869ab Binary files /dev/null and b/data2/figures/r_peak_timeseries.png differ diff --git a/data2/results/calibration_results.csv b/data2/results/calibration_results.csv new file mode 100644 index 0000000..3aa5394 --- /dev/null +++ b/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 diff --git a/data2/results/multivariate_results.csv b/data2/results/multivariate_results.csv new file mode 100644 index 0000000..197fc49 --- /dev/null +++ b/data2/results/multivariate_results.csv @@ -0,0 +1,2 @@ +PCR_RMSE_percent,PCR_R2 +0.067618748601777,0.987530104102354 diff --git a/data2/results/pca_explained_variance.csv b/data2/results/pca_explained_variance.csv new file mode 100644 index 0000000..d830565 --- /dev/null +++ b/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 diff --git a/data2/results/plateau_statistics.csv b/data2/results/plateau_statistics.csv new file mode 100644 index 0000000..7eef592 --- /dev/null +++ b/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