Skip to content
This repository was archived by the owner on Jun 19, 2024. It is now read-only.

Commit

Permalink
first draft of spherical transformation algorithm
Browse files Browse the repository at this point in the history
  • Loading branch information
hbartle committed Dec 13, 2018
1 parent 4f6992a commit 42499ef
Show file tree
Hide file tree
Showing 7 changed files with 399 additions and 0 deletions.
18 changes: 18 additions & 0 deletions misc_functions/associatedLegendre.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,18 @@
function [P] = associatedLegendre(n,x)
%ASSOCIATEDLEGENDRE
% Modified associated legendre function to duplicate negative m values
%
% Input Arguments:
%
% n Degree
% x Argument
%
% Output Arguments
%
% P Polynomial Value

p = legendre(n,x);
P = [flipud(p(2:end,:));p]' ;

end

41 changes: 41 additions & 0 deletions misc_functions/associatedLegendreDerivative.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,41 @@
function [dP] = associatedLegendreDerivative(n,x)
%ASSOCIATEDLEGENDREDERIVATIVE
% First derivative of associated legendre function using recursion formula
% and boundary conditions
%

m = -n:n;
M = repmat(m,length(x),1);
X = repmat(x,1,length(m));
Pnm = associatedLegendre(n,x);

DP1 = Pnm;

% Handle m=0
i = find(m==0);
DP2 = [Pnm(:,i)./sqrt(1-x.^2) Pnm(:,i+1:end)];
DP2 = [fliplr(DP2(:,2:end)) DP2];

dP = 1./(1-X.^2) .* (n*x.*DP1 -(n+M).*(n-M+1).*sqrt(1-X.^2).*DP2);

% Handle boundaries at x = +/-1
i1 = find(abs(x+1)<0.0000001);
i2 = find(abs(x-1)<0.0000001);

k0 = find(abs(m)==0);
k1 = find(abs(m)==1);
k2 = find(abs(m)==2);
k3 = find(abs(m)>=3);

dP(i1,k0) = n*(n+1)/2;
dP(i1,k1) = Inf;
dP(i1,k2) = -(n-1)*n*(n+1)*(n+2)/4;
dP(i1,k3) = 0;

dP(i2,k0) = (-1)^(n+1)*n*(n+1)/2;
dP(i2,k1) = (-1)^n*Inf;
dP(i2,k2) = (-1)^n*(n-1)*n*(n+1)*(n+2)/4;
dP(i2,k3) = 0;

end

31 changes: 31 additions & 0 deletions misc_functions/rotateSphericalNFData.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
function [data_nf_rotated] = rotateSphericalNFData(data_nf)

% Create Results Table
p = length(data_nf.E);
data_nf_rotated = table(zeros(p,1),zeros(p,1),zeros(p,1),zeros(p,1),zeros(p,1),zeros(p,1),zeros(p,3),zeros(p,1));
data_nf_rotated.Properties.VariableNames = {'x','y','z','r','theta','phi','E','Eabs'};

data_nf_rotated.x = data_nf.x;
data_nf_rotated.y = data_nf.y;
data_nf_rotated.z = data_nf.z;

r = sqrt(data_nf.x.^2 + data_nf.y.^2 + data_nf.z.^2);
theta = acos(data_nf.z./r);
phi = atan2(data_nf.y,data_nf.x);
for i =1:length(data_nf.E)
R = [sin(theta(i))*cos(phi(i)) sin(theta(i))*sin(phi(i)) cos(theta(i));...
cos(theta(i))*cos(phi(i)) cos(theta(i))*sin(phi(i)) -sin(theta(i));...
-sin(phi(i)) cos(phi(i)) 0];

%E(1)-> Radial,E(2)->Theta,E(3)->Phi
data_nf_rotated.E(i,:) = (R*data_nf.E(i,:)')';
data_nf_rotated.r(i) = r(i);
data_nf_rotated.theta(i) = theta(i);
if phi(i) <0
data_nf_rotated.phi(i) = pi - phi(i);
else
data_nf_rotated.phi(i) = phi(i);
end
data_nf_rotated.Eabs(i) = data_nf.Eabs(i);
end

60 changes: 60 additions & 0 deletions misc_functions/sphericalVectorWaveFunction.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,60 @@
function [fr,ftheta,fphi, Y] = sphericalVectorWaveFunction(s,m,n,A,theta,phi,k)
%SPHERICALVECTORWAVEFUNCTION
%
% Input Arguments:
%
% s Parity
% m,n Iterators over
% A Radius
% theta Theta angle
% phi Phi angle
% k Wavenumber at given frequency
%
% Output Arguments:
%
% fr Radial Component of Spherical wave function
%
% ftheta Theta Component of Spherical wave function
%
% fphi Phi Component of Spherical wave function
%
% Y Factor that shows up when solving for the spherical wave
% expansion coefficients
%

% Helper Functions
z =@(n,r) sqrt(pi/(2*r)).* besselh(n+0.5,r);
dz = @(n,r) z(n-1,r) - (n+1)/r * z(n,r); % CHECK IF THIS IS CORRECT!!
delta = @(s,sigma) 0.5*(1+(-1)^(s+sigma));


%%%% Note %%%
% R,T and P provide values for each |m| to speed up the process
% This means the f output will be a matrix with m columns

% Radial Function
R = @(s,n,r) delta(s,1).*z(n,k*r) + delta(s,2).*1/(k*r).* dz(n,k*r);
% Theta Function
T = @(s,m,n,theta) delta(s,1)*1j*m./sin(theta) .* associatedLegendre(n,cos(theta))...
+ delta(s,2) .* associatedLegendreDerivative(n,cos(theta)); % Add derivative of Associated Legendre Function;

% Phi Function
P = @(m,phi) exp(1j*m.*repmat(phi,1,size(m,2)));



% Radial Component
fr = delta(s,2)*n*(n+1)/(k*A) *z(n,k*A) * associatedLegendre(n,cos(theta)).*P(m,theta);
% Theta Component
ftheta = repmat(R(s,n,A),1,size(m,2)).*T(s,m,n,theta).*P(m,phi);
% Phi Component
fphi = (-1)^s.*repmat(R(s,n,A),1,size(m,2)).*T(s+1,m,n,theta).*P(m,theta);



% Additional Factor that shows up when solving for the spherical wave
% expansion coefficients.
Y = 4*pi/(2*n + 1) *factorial(n+abs(m(1,:)))./factorial(n-abs(m(1,:))) *n*(n+1)*R(s,n,A).^2;

end

49 changes: 49 additions & 0 deletions misc_functions/sphericalVectorWaveFunctionFarField.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,49 @@
function [xtheta,xphi] = sphericalVectorWaveFunctionFarField(s,m,n,A,theta,phi,k)
%SPHERICALVECTORWAVEFUNCTION
% Far-Field approximation of the spherical vector wave function
%
% Input Arguments:
%
% s Parity
% m,n Iterators over
% A Radius
% theta Theta angle
% phi Phi angle
% k Wavenumber at given frequency
%
% Output Arguments:
%
% xtheta Theta Component of Spherical wave function
%
% xphi Phi Component of Spherical wave function
%


% Helper Functions
delta = @(s,sigma) 0.5*(1+(-1)^(s+sigma));


%%%% Note %%%
% T and P provide values for each |m| to speed up the process
% This means the f output will be a matrix with m columns

% Theta Function
T = @(s,m,n,theta) delta(s,1)*1j*m./sin(theta) .* associatedLegendre(n,cos(theta))...
+ delta(s,2).*associatedLegendreDerivative(n,cos(theta));
% Phi Function
P = @(m,phi) exp(1j*m.*repmat(phi,1,size(m,2)));




C = 1/(k*A) * (-1j)^(n+delta(s,1)) * exp(1j*k*A) *P(m,phi);
% Theta Component
xtheta = C .* T(s,m,n,theta);
% Phi Component
xphi = C .* (-1)^s .* T(s+1,m,n,theta);




end

125 changes: 125 additions & 0 deletions nf2ff_spherical.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,125 @@
%
% NF2FF Conversion using Spherical Scanner
%
clear
close all
clc

disp('************************************************')
disp(' Near-To-Far-Field Conversion')
disp(' Spherical Scanner')
disp('************************************************')

addpath('misc_functions')
addpath('plot_functions')
addpath('transformation_functions')
%% Load Data
disp('Load Data...')
f = 5e9;

setup = '../measurement/SphericalScan_HornAntenna/';


data_ff = readtable([setup 'Farfield/farfield (f=' num2str(f*1e-9) ') [1].txt']);

scans = dir(setup);
scan_names = {scans.name};
scan_names(ismember(scan_names,{'.','..','Farfield','Scan-NTheta-NPhi.txt'})) = [];

data_nf= cellfun(@(scan_names) readtable([setup,scan_names,'/NearFieldProbeResults' num2str(f*1e-9) 'GHz.txt']),scan_names,'UniformOutput',false);

% Rearrange farfield table
data_ff = table(data_ff.Var1*pi/180,data_ff.Var2*pi/180,data_ff.Var3,data_ff.Var4,data_ff.Var6);
data_ff.Properties.VariableNames = {'theta' 'phi' 'Eabs' 'Ethetaabs' 'Ephiabs'};

% Rearrange nearfield table
data_nf = cellfun(@rearrangeTables,data_nf,'UniformOutput',false);
data_nf = cellfun(@rotateSphericalNFData,data_nf,'UniformOutput',false);

% Select measurements to process
data_nf = data_nf([2]);
scan_names = scan_names([2]);

disp('Done!')
%% NF2FF transformation
disp('NF2FF Transformation...')
delta_theta=1;
theta_range = (0:delta_theta:180)*pi/180;
delta_phi = 1;
phi_range = (0:delta_phi:360)*pi/180;


data_nf2ff = cellfun(@(data_nf) nf2ff_spherical_manual(data_nf,f,theta_range,phi_range),data_nf,'Uniformoutput',false);

disp('Done!')
%% Plots
disp('Plotting...')
close all

normalized = true;
logarithmic = true;

% Phi=0 cut
figure('name','Far-Field Cuts,Phi=0°','numbertitle','off',...
'units','normalized','outerposition',[0 0 1 1]);
plotFFPhiCut(data_ff,0,normalized,logarithmic)
cellfun(@(data_nf2ff) plotNFPhiCutCylindrical(data_nf2ff,0,normalized,logarithmic),data_nf2ff)
grid on
xlabel('Theta [°]')
if logarithmic == true
ylim([-50 0])
ylabel('E-Field Pattern [dB]')
elseif normalized == true
ylim([0 1])
ylabel('E-Field Pattern [-]')
end
ylabel('E-Field Pattern [-]')
title('Far-Field Cut Phi=0°')
legend(['Far-Field',scan_names])


figure('name','Far-Field Error,Phi=0°','numbertitle','off',...
'units','normalized','outerposition',[0 0 1 1]);
cellfun(@(data_nf2ff) plotDiffPhiCutCylindrical(data_nf2ff,data_ff,0,theta_range),data_nf2ff)
grid on
xlabel('Theta [°]')
ylabel('Difference to Reference Far-Field [dB]')
title('Difference to Reference Far-Field, Phi=0°')
legend(scan_names)


% Theta=90 cut
figure('name','Far-Field Cuts,Theta=90°','numbertitle','off',...
'units','normalized','outerposition',[0 0 1 1]);
plotFFThetaCut(data_ff,pi/2,normalized,logarithmic)
cellfun(@(data_nf2ff) plotNFThetaCutCylindrical(data_nf2ff,pi/2,normalized,logarithmic),data_nf2ff)
grid on
xlabel('Phi [°]')
if logarithmic == true
ylim([-50 0])
ylabel('E-Field Pattern [dB]')
elseif normalized == true
ylim([0 1])
ylabel('E-Field Pattern [-]')
end
ylabel('E-Field Pattern [-]')
title('Far-Field Cut Theta=90°')
legend(['Far-Field',scan_names])


figure('name','Far-Field Error,Theta=90°','numbertitle','off',...
'units','normalized','outerposition',[0 0 1 1]);
cellfun(@(data_nf2ff) plotDiffThetaCutCylindrical(data_nf2ff,data_ff,pi/2),data_nf2ff)
grid on
xlabel('Phi [°]')
ylabel('Difference to Reference Far-Field [dB]')
title('Difference to Reference Far-Field, Theta=90°')
legend(scan_names)


disp('Done!')





Loading

0 comments on commit 42499ef

Please sign in to comment.