APPENDIX E
MATLAB/OCTAVE .M-FILE SOURCE CODE
::::::::::::::
linbrd.m
::::::::::::::
function [spec] = linbrd(d,lb)
%Linbrd
% Accepts Tecmag spectra raw data and outputs processed spectra.
% the tecmag vector is in d, lb Hz line broadening is applied.
%
% RJM 8 1 97
sr=d(:,1); % Real
si=d(:,2); % Imag
t=.001*d(:,3); % tecmag time is in msec
fid= sr + j*si;
fidb=fid-mean(fid);
fidlb=fidb.*exp(-lb*pi*t);
spec=fftshift(fft(fidlb));
%That's all folks.
endfunction
::::::::::::::
autophas.m
::::::::::::::
function [y0,y1,th0,th1] = autophas(x);
% AUTOPHAS
% [y0,y1,th0,th1] = autophas(x);
% Autophase performs automatic zero and first order phase correction
% on the complex spectrum in x. The criterion is to maximize the sum of
% the real components of x. y0 is the spectrum with zero order correction
% applied, y1 is the spectrum with first order correction. th0 and th1 are the
% corresponding phase angle and phase coefficient.
%
% Roger J. McNichols
% Texas A&M University
% 12-19-97
n = length(x);
x = x(:);
% zero order correction.
a0 = sum((real(x)));
b0 = sum(imag(x));
th0=atan(-b0/a0);
y0 = -x*exp(j*th0);
% First order correction...
w = [linspace(-pi,pi,n)]';
a1 = sum(real(y0).*cos(w));
b1 = sum(real(y0).*sin(w));
c1 = sum(imag(y0).*sin(w));
d1 = sum(imag(y0).*cos(w));
th1 = -atan((a1-d1)/(b1+c1));
y1 = y0.*exp(j*w*th1);
y0 = y0;
y1 = y1;
% That's all folks!
endfunction
::::::::::::::
linreg.m
::::::::::::::
function [a,b,R,yhat,resid] = linreg(x,y)
% LINREG
% [a,b,R,yhat,resid] = linreg(x,y) performs a linear regression to fit
% y = a*x + b in the least squares sense. Observations
% should appear in y, the independent variable in x. Returned value
% R is the correlation coefficient (not r^2).
% If x and y are both NxM matrices, then M linear regressions are
% performed fitting columns of y to corresponding columns of x.
% a,b, and R are then rowvectors of length M.
%
% The optional return arguments yhat and resid are the fits (predicted
% y values) and the residuals (y-yhat). Both yhat and resid are NxM
% matrices.
%
%
% Roger McNichols
% Texas A&M University
% 11/10/94
% Check to make sure matrices are of same size. If not an error will occur.
z = x+y;
% calculate regression coefficients
[n,m] = size(x);
xbar = ones(n,1)*mean(x);
ybar = ones(n,1)*mean(y);
Sxx = sum(x.^2) - 1/n*(sum(x)).^2;
Syy = sum(y.^2) - 1/n*(sum(y)).^2;
Sxy = sum(x.*y) - 1/n*(sum(x).*sum(y));
a = Sxy./Sxx;
b = mean(y) - a.*mean(x);
R = Sxy./((Sxx.^.5).*(Syy.^.5));
% Calculate residuals, fits, etc.
yhat = ones(n,1)*b + (ones(n,1)*a).*x;
resid = y-yhat;
::::::::::::::
smooth.m
::::::::::::::
function [y] = smooth(x,n)
% SMOOTH
% y=smooth(x,n) performs n-point smoothing on spectra contained in the
% vector x. If x is a matrix, then smoothing is performed column-wise.
% Before smoothing, x is padded with its edge values on both ends so that the
% output spectra y are of the same length as x. This should be taken
% into consideration when using information near the ends of the spectra.
%
% RJM 5/8/94
% Marcel's old smoothing routine is now called msmooth.
p=floor(n/2);
q=n-p;
[i,j]=size(x);
xt = [ ones(p,1)*x(1,:);x;ones(q,1)*x(i,:)];
sux=sum(xt(1:n,:));
y=sux/n;
for k = 2:i,
sux=sux-xt(k-1,:)+xt(k+p,:);
y=[y; sux/n];
end
::::::::::::::
loada.m
::::::::::::::
function x = loada(filename,c)
% LOADA
% var = loada('filename',r) performs an ascii text file load in a
% way similar to MATLAB'S 'load -a' command. The OCTAVE ascii format
% is not required. Rather ascii data are read directly as rows and
% columns of a matrix having c columns.
%
% RJM 2/17/98
fp=fopen(filename,'r');
if (fp == -1)
printf("Invalid file.\n");
return
endif
x = fscanf(fp,'%f',[c Inf]);
x = x';
fclose(fp);
endfunction
% That's all folks!