% Demonstration of the limits of Classical Least Squares (multilinear % regression) of two very closely-spaced "noiseless" overlapping Gaussian % peaks caused by the numerical precision of the computer software. This % demonstrates three different mathematical formulation of the % least-squares calculation that give different results when the numerical % precision limits of the computer are reached. But practically, the % difference between these methods is unlikey to be seen; even the tiniest % bit of added random noise (line 15) or signal instability produces a % far greater error. format compact format long g clear clf ADCbits=12; noise=0.0; % Also try small anounts of noise like .00000001; t=[1:10:1000]; amp=[1 .9]; % Amplitudes of the peaks pos=[500 511]; % Positions of the peaks wid=[300 300]; % Widths of the peaks % A = matrix containing one unit-amplidude peak in each of its rows A = zeros(length(pos),length(t)); for k=1:length(pos) A(k,:)=gaussian(t,pos(k),wid(k)); % You can use any peak function here end spectrum=amp*A; % Multiplies each row by the corresponding amplitude and adds them up clf % Simulation ObservedSpectrum = spectrum + noise.*randn(size(t)); % adds noise hold on % Regress ObservedSpectrum onto matrix of reference spectra. Can be done % bt several different mathematical methods that give the same results. % Using straight multilinear least-squares regression formulation tic MeasuredAmpInverse = ObservedSpectrum*A'*inv(A*A'); InverseTime=toc % Using pseudoinverse tic MeasuredAmpPseudoinverse = ObservedSpectrum*pinv(A); PseudoinverseTime=toc % Using matrix division tic MeasuredAmpMatrixDivision = ObservedSpectrum/A; MatrixDivisionTime=toc figure(1);plot(t,ObservedSpectrum,'m.') plot(t,amp(1).*A(1,:),t,amp(2).*A(2,:)) plot(t,MeasuredAmpInverse(1).*A(1,:),t,MeasuredAmpInverse(2).*A(2,:)) plot(t,MeasuredAmpInverse(1).*A(1,:)+MeasuredAmpInverse(2).*A(2,:),'k') title('Analysis of mixture of peaks by multiple regression') xlabel('Dots: Measured data; Colored lines are true and measured component spectra (blue and green)') hold off disp(' ') SeparationToWidthRatio=(pos(2)-pos(1))./mean(wid) disp(['Method 1, the inverse method: MeasuredAmp=ObservedSpectrum*A''*inv(A*A'')'] ) disp(MeasuredAmpInverse) ErrorInverse = MeasuredAmpInverse-amp; RelativePercentErrorInverse = 100.*ErrorInverse./amp disp('--------------------------------------------------------------------') disp('Method 2, the pseudo-inverse method: MeasuredAmp=ObservedSpectrum*pinv(A)' ) disp(MeasuredAmpPseudoinverse) ErrorPseudoinverse = MeasuredAmpPseudoinverse-amp; RelativePercentErrorPseudoinverse = 100.*ErrorPseudoinverse./amp disp('--------------------------------------------------------------------') disp('Method 2, the matrix division method: MeasuredAmp=ObservedSpectrum/A' ) disp(MeasuredAmpMatrixDivision) ErrorMatrixDivision = MeasuredAmpMatrixDivision-amp; RelativePercentErrorMatrixDivision = 100.*ErrorMatrixDivision./amp