/usr/share/psychtoolbox-3/PsychCal/CalibrateFitLinMod.m is in psychtoolbox-3-common 3.0.11.20140816.dfsg1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 | function cal = CalibrateFitLinMod(cal)
% cal = CalibrateFitLinMod(cal)
%
% Fit the linear model to spectral calibration data.
%
% 3/26/02 dhb Pulled out of CalibrateMonDrvr.
% 3/27/02 dhb Add case of nPrimaryBases == 0.
% 2/15/10 dhb Fix so that first basis vector is good approximation to max
% input primary spectrum.
% dhb Normalize basis vectors so that their max power matches that
% of first component.
% 4/30/10 dhb Execute yoked fit if yokedGamma flag is set.
% 5/25/10 dhb, ar Change yoked field names to match
% 5/26/10 dhb, ar Still fussing with names.
% 5/28/10 dhb, ar Pull out yoked fitting from here -- too confusing.
% 5/27/12 dhb Handle case where there are more measurements than wavelength samples.
Pmon = zeros(cal.describe.S(3),cal.nDevices*cal.nPrimaryBases);
mGammaRaw = zeros(cal.describe.nMeas,cal.nDevices*cal.nPrimaryBases);
monSVs = zeros(min([cal.describe.nMeas cal.describe.S(3)]),cal.nDevices);
for i = 1:cal.nDevices
tempMon = reshape(cal.rawdata.mon(:,i),cal.describe.S(3),cal.describe.nMeas);
monSVs(:,i) = svd(tempMon);
% Build a linear model
if (cal.nPrimaryBases ~= 0)
% Get full linear model
[monB,monW] = FindLinMod(tempMon,cal.nPrimaryBases);
% Express max measurement within the full linear model.
% This is the first basis function.
tempB = monB*monW(:,cal.describe.nMeas);
maxPow = max(abs(tempB));
% Get residuals with respect to first component
residMon = tempMon-tempB*(tempB\tempMon);
% If linear model dimension is greater than 1,
% fit linear model of dimension-1 to the residuals.
% Take this as the higher order terms of the linear model.
%
% Also normalize each basis vector to max power of first
% component, and make sure convention is that this max
% is positive.
if (cal.nPrimaryBases > 1)
residB = FindLinMod(residMon,cal.nPrimaryBases-1);
for j = 1:cal.nPrimaryBases-1
residB(:,j) = maxPow*residB(:,j)/max(abs(residB(:,j)));
[nil,index] = max(abs(residB(:,j)));
if (residB(index,j) < 0)
residB(:,j) = -residB(:,j);
end
end
monB = [tempB residB];
else
monB = tempB;
end
% Zero means build one dimensional linear model just taking max measurement
% as the spectrum.
else
cal.nPrimaryBases = 1;
monB = tempMon(:,cal.describe.nMeas);
end
% Find weights with respect to adjusted linear model and
% store
monW = FindModelWeights(tempMon,monB);
for j = 1:cal.nPrimaryBases
mGammaRaw(:,i+(j-1)*cal.nDevices) = (monW(j,:))';
Pmon(:,i+(j-1)*cal.nDevices) = monB(:,j);
end
end
% Update calibration structure.
cal.S_device = cal.describe.S;
cal.P_device = Pmon;
cal.T_device = WlsToT(cal.describe.S);
cal.rawdata.rawGammaTable = mGammaRaw;
cal.rawdata.monSVs = monSVs;
|