1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
|
function [delay, fidelity, nf2ff_out] = DelayFidelity(nf2ff, port, path, weight_theta, weight_phi, theta, phi, f_0, f_c, varargin)
% [delay, fidelity] = DelayFidelity(nf2ff, port, path, weight_theta, weight_phi, theta, phi, f_lo, f_hi, varargin)
%
%
% This function calculates the time delay from the source port to the phase center of the antenna and the fidelity.
% The fidelity is the similarity between the excitation pulse and the radiated pulse (normalized scalar product).
% The resolution of the delay will be equal to or better than ((f_0 + f_c)*Oversampling)^-1 when using Gaussian excitation.
% Oversampling is an input parameter to InitFDTD. The rows of delay and fidelity correspond to theta and the columns to phi.
%
% input:
% nf2ff: return value of CreateNF2FFBox.
% port: return value of AddLumpedPort
% path: path of the simulation results.
% weight_theta: weight of the E_theta component
% weight_phi: weight of the E_phi component
% -> with both (possibly complex) parameters any polarization can be examined
% theta: theta values to be simulated
% phi: phi values to be simulated
% f_0: center frequency of SetGaussExcite
% f_c: cutoff frequency of SetGaussExcite
%
% variable input:
% 'Center': phase center of the antenna for CalcNF2FF
% 'Radius': radius for CalcNF2FF
% 'Mode': mode CalcNF2FF
%
% example:
% theta = [-180:10:180] * pi / 180;
% phi = [0, 90] * pi / 180;
% % use circular right handed polarization
% [delay, fidelity] = DelayFidelity2(nf2ff, port, Sim_Path, -1i, 1, theta, phi, f_0, f_c, 'Mode', 1);
% figure
% polar(theta.', delay(:,1) * 3e11); % delay in mm
% figure
% polar(theta', (fidelity(:,1)-0.95)/0.05); % last 5 percent of fidelity
%
% Author: Georg Michel
C0 = 299792458;
center = [0, 0, 0];
radius = 1;
nf2ff_mode = 0;
for n=1:2:numel(varargin)
if (strcmp(varargin{n},'Center')==1);
center = varargin{n+1};
elseif (strcmp(varargin{n},'Radius')==1);
radius = varargin{n+1};
elseif (strcmp(varargin{n},'Mode')==1);
nf2ff_mode = varargin{n+1};
end
end
port_ut = load(fullfile(path, port.U_filename));
port_it = load(fullfile(path, port.I_filename));
dt = port_ut(2,1) - port_ut(1,1);
fftsize = 2^(nextpow2(size(port_ut)(1)) + 1);
df = 1 / (dt * fftsize);
uport = fft(port_ut(:, 2), fftsize)(1:fftsize/2+1);
iport = fft(port_it(:, 2), fftsize)(1:fftsize/2+1);
fport = df * (0:fftsize/2);
f_ind = find(fport > (f_0 - f_c ) & fport < (f_0 + f_c));
disp(["frequencies: ", num2str(numel(f_ind))]);
exc_f = uport.' + iport.' * port.Feed_R; %excitation in freq domain
exc_f(!f_ind) = 0;
exc_f /= sqrt(exc_f * exc_f'); % normalization (transposing also conjugates)
nf2ff = CalcNF2FF(nf2ff, path, fport(f_ind), theta, phi, ...
'Center', center, 'Radius', radius, 'Mode', nf2ff_mode);
radfield = weight_theta * cell2mat(nf2ff.E_theta) + weight_phi * cell2mat(nf2ff.E_phi); % rows: theta(f1), columns: phi(f1), phi(f2), ...phi(fn)
radfield = reshape(radfield, [length(nf2ff.theta), length(nf2ff.phi), length(nf2ff.freq)]);
correction = reshape(exp(-2i*pi*nf2ff.r/C0*nf2ff.freq), 1,1,numel(nf2ff.freq)); %dimensions: theta, phi, frequencies
radfield = radfield./correction; % correct for radius delay
% normalize radfield
radnorm = sqrt(dot(radfield, radfield, 3));
radfield ./= radnorm;
%initialize radiated field in fully populated frequency domain
rad_f = zeros([numel(nf2ff.theta), numel(nf2ff.phi), numel(fport)]);
rad_f(:, :, f_ind) = radfield; % assign selected frequencies
exc_f = reshape(exc_f, [1,1,numel(exc_f)]); %make exc_f confomant with rad_f
cr_f = rad_f .* conj(exc_f); % calculate cross correlation
% calculate the cross correlation in time domain (analytic signal)
cr = ifft(cr_f(:, :, 1:end-1), [], 3) * (numel(fport) -1); % twice the FFT normalization (sqrt^2) because product of two normalized functions
%search for the maxiumum of the envelope
[fidelity, delay_ind] = max(abs(cr), [], 3);
delay = (delay_ind - 1) * dt * 2; % double time step because of single-sided FFT
nf2ff_out = nf2ff; %possibly needed for plotting the far field and other things
disp(["DelayFidelity: delay resolution = ", num2str(dt*2e9), "ns"]);
return;
|
