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!