fresnel.m

near field localization

function fresnel(ks, doa)
%
% by Joseph Handfield
% This function calculates near field source fields and inputs the field
% into the Capon function. The DOA, range and intensity of each sensed
% source is returned to this function for MSE analysis with the simulated
% input.
%______________________________________________________________________
dataType = 1; % 1 = Actual Correlation, 2 = Theoretical Correlation (not used)
msg2 = strcat(['Your sensor source D must be <= ' num2str(2*(ks/2)^2) 'wl.']);
disp(msg2);
% The following line can locate the point sources of the test signal along
% an arc of fixed range equal to sfnd(2). Otherwise the sources are located
% along a parallel source plane/line at a distance away determined by sfnd(2).
arc = 0; % input('Fixed range mode for test signal? 1=Yes:');
% sfnd will be a 1 X 2 array of system parameters
sfnd = input('Enter source field width and distance from sensor to source line.');
% loc will be n wide depending on number of sources
loc = input('Enter source locations in array form.');
% The next step will translate the test source locations if a certain DOA is desired
loc = loc + sfnd(2)*tan(doa)
if isempty(loc), error('No location data'); end
% if any(abs(loc) > sfnd(1)/2), error('Need a wider source field'); end
% If simulating non unity intensities for the sources, uncomment next
% line and remove []
val = [1/2 1 1 1/2]; % input('Enter the respective source intensities.');
sources = length(loc); % How many source points are being simulated
if isempty(val) | (length(val)~= length(loc))
val = ones(1,sources);
end
index = -ks:ks;
snsrs = length(index);
phases2 = [];
sensArray = [];
senstotal = 0;
covSum = zeros(snsrs);
rate = zeros(snsrs,sources);
sens = zeros(snsrs,sources);
angleA = atan(loc./sfnd(2))
% The following lines are true for angles that pertain to source line spacing
if arc == 1
% Lock in broadside distance from sensor to source line sfnd(2) = range
rate = pi./(16*sfnd(2))*(index').^2*(cos(angleA)).^2 ;
else
% Normal mode : cos(theta)/sfnd(2) = stretching radius along source line
rate = pi./(16*sfnd(2))*(index').^2*(cos(angleA)).^3 ;
end
noisy = input('Enter SNR for additive noise, or enter _ for no noise:');
records = 140; %input('Number of records to average?');
% Now to generate samples every 1/10 of a period
timeSamp = exp(i*2*pi/10*(1:10));
if dataType == 1
for p= 1:records
phases = [];
sens = 0;
for u = 1:sources
% Separate random phases will be used to 'unlock' simulated sources
rphase = (rand(1)*2 - 1) * pi;
part = exp(i*(-pi/2*sin(angleA(u))*index' + rate(:,u)));
% In the next step, a random phase and time samples are
% generated for the source signal making 'sens' an array
% that is # of sensors -by- 10 time samples
npart = part*timeSamp*exp(i*rphase)*val(u);
if isempty(noisy) ~= 1
npart = addNoise(npart,[snsrs 10],noisy);
end
sens = sens + npart;
% Extra sources are added on with different phases
phases = [phases rphase];
end
corrAvg = 0;
% Now to time average the correlation of each time sample:
for ts = 1:10
corrAvg = corrAvg + sens(:,ts)*ctranspose(sens(:,ts));
end
% Now, to time average the sum of the correlations:
covSum = covSum + corrAvg/10;
phases2 = [phases2; phases];
end
% Now, to ensemble average the sum of the time correlations at random phases:
covFinal = covSum/records;
end
% Now to calculate the Capon weights and frequency matrix
srcSep = input('What is the distance between positions on the source line?');
% Capon2 takes the correlation matrix and the last sensor output over time
fresp = capon2(index,sens,sfnd,covFinal,srcSep);
if fresp ~= 0
% These lines produce a crude source line reconstruction, so
% errors can be computed.
spos = sfnd(1)/srcSep;
pvals = length(loc);
srcLine = zeros(spos,1);
spacing = (0:spos)*srcSep - sfnd(1)/2;
for yval = 1:pvals
p = find(spacing == loc(yval));
srcLine(p) = val(yval);
p=0;
end
avals = size(fresp,1);
reconst = zeros(spos,1);
for ansval = 1:avals
resp = find(spacing == sfnd(2)*tan(fresp(ansval,1)));
reconst(resp) = fresp(ansval,3);
end
errs = reconst - srcLine;
disp('The mean squared error is:');
meanSqEr = mse(errs)
disp('end'); disp(' ');
end
prompt1 = input('Another try?');
if prompt1 == 1
fresnel(ks, doa);
end
function xout = addNoise(ex,xn,snr)
%
% by Joseph Handfield
% ex = values at xv source positions
% xn = number of reconstruction points & number of time samples
% snr = signal to noise ratio
if nargin < 3, error('Not enough input arguments.'); end
% add some random noise according to the snr term
noiseVar = max(max(abs(ex))) / 10^(snr/10); % From Reisenfeld & Aboutanios
% The below takes into consideration that the real and imaginary noise proportion is
random
% randPct = rand(1);
noise = sqrt(noiseVar)*randn(xn(1)*2, xn(2));
realnoise = noise(1:xn(1),:);
imagnoise = noise(xn(1)+1:xn(1)*2,:)*i;
xout = ex + realnoise + imagnoise;
function sArray = capon2(index,sens,sfnd,cov,srcSep)
%
% Written 5/9/07 by: Joseph Handfield
% Objective of this function is to execte the Capon Spectral Analysis
% method on input field data to locate nearfield sources by their
% incident angle from the central sensor and their radial distance
% from that sensor. This is done using a manifold array of specially-built
% frequency filters that combines with the input field signal at various
% values of theta and r. The ideal value(s) of xc are approximately one.
%_______________________________________
% 'guess' will tell the algorithm to search the gradient of the spectra for
% a number of peaks
guess = 3; %input('How many sources are you assuming?');
if isempty(guess) | guess == 0, error('No Sources'); end
fov= sfnd(1); distL = sfnd(2);
spos = fov/srcSep;
offs = 0; %input('offset to source line spacing?');
if isempty(offs), offs = 0; end
plotmode = 1; % Enter 1 for mesh/contour plot, 2 for normal 2D plot of a fixed range arc
% Corresponding radii in wavelengths, change the below fraction to modify sampling
distance
if plotmode == 1
rad = ((0:50)/1000 + 0.98)*sqrt(distL^2 + offs^2);
else
rad = distL;
end
spacing = (0:spos)*srcSep - fov/2 + offs;
% thetas is set to angles that pertain to the spacing of the source line
thetas = atan(spacing/distL);
invcov = inv(cov);
counter = 0;
%This is here in case diagonal loading is needed (singular correlation matrix)
if any(invcov) == inf
while estimate < 0.5e-16
counter = counter + 1;
cov = cov + eye(length(cov))*0.2; %Amount added can be adjusted to converge
quicker
estimate = rcond(cov)
end
invcov = pinv(cov);
end
% The next few lines generate the MVDR (Capon?) spectrum
xc = zeros(length(thetas),length(rad));
minvar = zeros(length(thetas),length(rad));
for a = 1:length(thetas)
for b = 1:length(rad)
% These calcultions assume that the sensor spacing is lambda/4
rate = pi/(16*rad(b))*(cos(thetas(a))*index').^2;
man = exp(i*(-pi/2*sin(thetas(a))*index' + rate));
minvar(a,b) = 1/ (ctranspose(man)*invcov*man);
weights = ctranspose(sens)*(invcov*man*minvar(a,b));
xc(a,b) = sum(abs(weights));
end
end
% Sobel operators, to create a gradient plot of the spectra
if plotmode == 1
B1 = [-1 -2 -1; 0 0 0; 1 2 1];
B2 = [-1 0 1; -2 0 2; -1 0 1];
u1 = filter2(B1,abs(xc),'same');
u2 = filter2(B2,abs(xc),'same');
gradPlot = abs(u1) + abs(u2);
disp('The Capon Matrix has solution(s) of:');
% Acccording to value in guess, this picks the maximum valued Capon peaks
% this was changed to locate all peaks when some may be many orders of magnitude higher
temp = log(abs(xc));
argBkgd = median(median(temp));
sxc = size(xc);
aset = [];
p=0; niter = 0;
manA = [];
while p ~= guess
niter = niter + 1;
soln = max(max(temp));
[a,b] = find(temp == soln);
if a ~= sxc(1)
rp = gradPlot(a,b) < 0.5*gradPlot(a+1,b);
else rp = 0; end
if a > 1
lp = gradPlot(a,b) < 0.5*gradPlot(a-1,b);
else lp = 0; end
if b ~= sxc(2)
up = gradPlot(a,b) < 0.5*gradPlot(a,b+1);
else up = 0; end
if b > 1
dp = gradPlot(a,b) < 0.5*gradPlot(a,b-1);
else dp = 0; end
if (rp & lp) & (up & dp)
soln
disp(['at theta = ' num2str(thetas(a)) ...
', and at a radius = ' num2str(rad(b))]);
rateA = pi/(16*rad(b))*(cos(thetas(a))*index').^2;
manA = [manA exp(i*(-pi/2*sin(thetas(a))*index' + rate))];
p = p+1;
% If the solution positions are to be recorded
aset = [aset; thetas(a) rad(b) 0];
end
if niter > 3*guess
disp('Valid Capon solutions are found');
p=guess;
end
temp(a,b) = argBkgd; % Reduce the solution peak so others can be located
end
pic = input('Would you like to create a mesh plot? Enter title.');
if isempty(pic) == 1
% do nothing
else
% figure;
% logMVDR = log(abs(xc));
% mesh(rad,thetas,logMVDR);
% axis tight;
% ylabel('DOA - radians');
% xlabel('range (wavelengths)');
% zlabel('Log Spectral Intensity');
% graph = strcat(['MVDR Filter Response ' pic]);
% title(graph);
% For Power Spectral Estimator (denominator of Capon solution)
figure;
mesh(rad,thetas,abs(minvar));
ylabel('DOA - radians');
xlabel('range (wavelengths)');
zlabel('Spectral Intensity');
axis tight;
graph2 = strcat(['Power Spectral Density of ' pic]);
title(graph2);
% figure;
% maxLog = fix(max(max(logMVDR)));
% levels = (-10:0) + maxLog;
% [cs, h]= contourf(rad,thetas,logMVDR,levels); clabel(cs,h);
% in2 = input('To save jpeg, type 1');
% if in2 == 1
% print(gcf, '-djpeg', pic);
% clear temp;
% save(pic);
% end
end
% Gaussian elimination division to find intensities of field points
if length(manA) ~= 0
intens = max(abs(sens'))/manA'
aset(:,3) = intens';
end
sArray = aset;
else
figure;
plot(thetas,log(abs(xc)));
xlabel('DOA - radians');
ylabel('Log Spectral Intensity');
sArray = 0;
end