Modelling a Physical Channel in the Time Domain

This example uses the frequency domain lyceanem.models.time_domain.calculate_scattering() function to predict the time domain response for a given excitation signal and environment included in the model. This model allows for a very wide range of antennas and antenna arrays to be considered, but for simplicity only horn antennas will be included in this example. The simplest case would be a single source point and single receive point, rather than an aperture antenna such as a horn.

import numpy as np

Frequency and Mesh Resolution

sampling_freq = 23e9
model_time = 1e-7
num_samples = int(model_time * (sampling_freq))

# simulate receiver noise
bandwidth = 8e9
kb = 1.38065e-23
receiver_impedence = 50
thermal_noise_power = 4 * kb * 293.15 * receiver_impedence * bandwidth
noise_power = -80  # dbw
mean_noise = 0

model_freq = 10e9
wavelength = 3e8 / model_freq

Setup transmitters and receivers

import lyceanem.geometry.targets as TL
import lyceanem.geometry.geometryfunctions as GF

transmit_horn_structure, transmitting_antenna_surface_coords = TL.meshedHorn(
    58e-3, 58e-3, 128e-3, 2e-3, 0.21, wavelength * 0.5
)
receive_horn_structure, receiving_antenna_surface_coords = TL.meshedHorn(
    58e-3, 58e-3, 128e-3, 2e-3, 0.21, wavelength * 0.5
)

Position Transmitter

rotate the transmitting antenna to the desired orientation, and then translate to final position. lyceanem.geometry.geometryfunctions.translate_mesh(), lyceanem.geometry.geometryfunctions.mesh_rotate() and lyceanem.geometry.geometryfunctions.mesh_transform() are included, allowing translation, rotation, and transformation of the meshio objects as required.

rotation_vector1 = np.radians(np.asarray([90.0, 0.0, 0.0]))
rotation_vector2 = np.radians(np.asarray([0.0, 0.0, -90.0]))


transmit_horn_structure = GF.mesh_rotate(transmit_horn_structure, rotation_vector1)
transmit_horn_structure = GF.mesh_rotate(transmit_horn_structure, rotation_vector2)
transmit_horn_structure = GF.mesh_translate(
    transmit_horn_structure, np.asarray([2.695, 0, 0])
)
transmitting_antenna_surface_coords = GF.mesh_rotate(
    transmitting_antenna_surface_coords, rotation_vector1
)
transmitting_antenna_surface_coords = GF.mesh_rotate(
    transmitting_antenna_surface_coords, rotation_vector2
)
transmitting_antenna_surface_coords = GF.mesh_translate(
    transmitting_antenna_surface_coords, np.asarray([2.695, 0, 0])
)

Position Receiver

rotate the receiving horn to desired orientation and translate to final position.

receive_horn_structure = GF.mesh_rotate(receive_horn_structure, rotation_vector1)
receive_horn_structure = GF.mesh_translate(
    receive_horn_structure, np.asarray([0, 1.427, 0])
)
receiving_antenna_surface_coords = GF.mesh_rotate(
    receiving_antenna_surface_coords, rotation_vector1
)
receiving_antenna_surface_coords = GF.mesh_translate(
    receiving_antenna_surface_coords, np.asarray([0, 1.427, 0])
)

Create Scattering Plate

Create a Scattering plate a source of multipath reflections

reflectorplate, scatter_points = TL.meshedReflector(
    0.3, 0.3, 6e-3, wavelength * 0.5, sides="front"
)
position_vector = np.asarray([29e-3, 0.0, 0])
rotation_vector1 = np.radians(np.asarray([0.0, 90.0, 0.0]))
scatter_points = GF.mesh_rotate(scatter_points, rotation_vector1)
reflectorplate = GF.mesh_rotate(reflectorplate, rotation_vector1)
reflectorplate = GF.mesh_translate(reflectorplate, position_vector)
scatter_points = GF.mesh_translate(scatter_points, position_vector)

meshing reflector
args 0.3 0.3 0.006
majorsize 0.3
minorsize 0.3
thickness 0.006

Specify Reflection Angle

Rotate the scattering plate to the optimum angle for reflection from the transmitting to receiving horn

plate_orientation_angle = 45.0

rotation_vector = np.radians(np.asarray([0.0, 0.0, plate_orientation_angle]))
scatter_points = GF.mesh_rotate(scatter_points, rotation_vector)
reflectorplate = GF.mesh_rotate(reflectorplate, rotation_vector)
from lyceanem.base_classes import structures

blockers = structures([reflectorplate, receive_horn_structure, transmit_horn_structure])

Visualise the Scene Geometry

import pyvista as pv


## plot the mesh
plotter = pv.Plotter()
plotter.add_mesh(pv.from_meshio(reflectorplate), color="white", show_edges=True)
plotter.add_mesh(pv.from_meshio(receive_horn_structure), color="blue", show_edges=True)
plotter.add_mesh(pv.from_meshio(transmit_horn_structure), color="red", show_edges=True)
plotter.add_axes_at_origin()
plotter.show()

Static Scene

Interactive Scene

Specify desired Transmit Polarisation

The transmit polarisation has a significant effect on the channel characteristics. In this example the transmit horn will be vertically polarised, (e-vector aligned with the z direction)

desired_E_axis = np.zeros((1, 3), dtype=np.float32)
desired_E_axis[0, 1] = 1.0

Time Domain Scattering

import scipy.signal as sig
import lyceanem.models.time_domain as TD
from lyceanem.base_classes import structures


angle_values = np.linspace(0, 90, 91)
angle_increment = np.diff(angle_values)[0]
responsex = np.zeros((len(angle_values)), dtype="complex")
responsey = np.zeros((len(angle_values)), dtype="complex")
responsez = np.zeros((len(angle_values)), dtype="complex")

plate_orientation_angle = -45.0

rotation_vector = np.radians(
    np.asarray([0.0, 0.0, plate_orientation_angle + angle_increment])
)
scatter_points = GF.mesh_rotate(scatter_points, rotation_vector)
reflectorplate = GF.mesh_rotate(reflectorplate, rotation_vector)

wake_times = np.zeros((len(angle_values)))
Ex = np.zeros((len(angle_values), num_samples))
Ey = np.zeros((len(angle_values), num_samples))
Ez = np.zeros((len(angle_values), num_samples))

for angle_inc in range(len(angle_values)):
    rotation_vector = np.radians(np.asarray([0.0, 0.0, angle_increment]))
    scatter_points = GF.mesh_rotate(scatter_points, rotation_vector)
    reflectorplate = GF.mesh_rotate(reflectorplate, rotation_vector)
    blockers = structures(
        [reflectorplate, transmit_horn_structure, receive_horn_structure]
    )
    pulse_time = 5e-9
    output_power = 0.01  # dBwatts
    powerdbm = 10 * np.log10(output_power) + 30
    v_transmit = ((10 ** (powerdbm / 20)) * receiver_impedence) ** 0.5
    output_amplitude_rms = v_transmit / (1 / np.sqrt(2))
    output_amplitude_peak = v_transmit

    desired_E_axis = np.zeros((1, 3), dtype=np.float32)
    desired_E_axis[0, 1] = 1.0
    noise_volts_peak = (10 ** (noise_power / 10) * receiver_impedence) * 0.5

    excitation_signal = output_amplitude_rms * sig.chirp(
        np.linspace(0, pulse_time, int(pulse_time * sampling_freq)),
        model_freq - bandwidth,
        pulse_time,
        model_freq,
        method="linear",
        phi=0,
        vertex_zero=True,
    ) + np.random.normal(mean_noise, noise_volts_peak, int(pulse_time * sampling_freq))
    (
        Ex[angle_inc, :],
        Ey[angle_inc, :],
        Ez[angle_inc, :],
        wake_times[angle_inc],
    ) = TD.calculate_scattering(
        transmitting_antenna_surface_coords,
        receiving_antenna_surface_coords,
        excitation_signal,
        blockers,
        desired_E_axis,
        scatter_points=scatter_points,
        wavelength=wavelength,
        scattering=1,
        elements=False,
        sampling_freq=sampling_freq,
        num_samples=num_samples,
        beta=(2 * np.pi) / wavelength,
    )

    noise_volts = np.random.normal(mean_noise, noise_volts_peak, num_samples)
    Ex[angle_inc, :] = Ex[angle_inc, :] + noise_volts
    Ey[angle_inc, :] = Ey[angle_inc, :] + noise_volts
    Ez[angle_inc, :] = Ez[angle_inc, :] + noise_volts

C:\Users\lycea\miniconda3\envs\CudaDevelopment\Lib\site-packages\lyceanem\electromagnetics\empropagation.py:3719: ComplexWarning: Casting complex values to real discards the imaginary part
  uvn_axes[2, :] = point_vector
C:\Users\lycea\miniconda3\envs\CudaDevelopment\Lib\site-packages\lyceanem\electromagnetics\empropagation.py:3736: ComplexWarning: Casting complex values to real discards the imaginary part
  uvn_axes[0, :] = np.cross(local_axes[2, :], point_vector) / np.linalg.norm(
C:\Users\lycea\miniconda3\envs\CudaDevelopment\Lib\site-packages\lyceanem\electromagnetics\empropagation.py:3758: ComplexWarning: Casting complex values to real discards the imaginary part
  uvn_axes[1, :] = np.cross(point_vector, uvn_axes[0, :]) / np.linalg.norm(
C:\Users\lycea\miniconda3\envs\CudaDevelopment\Lib\site-packages\numba_cuda\numba\cuda\dispatcher.py:693: NumbaPerformanceWarning: Grid size 25 will likely result in GPU under-utilization due to low occupancy.
  warn(NumbaPerformanceWarning(msg))
C:\Users\lycea\miniconda3\envs\CudaDevelopment\Lib\site-packages\lyceanem\electromagnetics\empropagation.py:3719: ComplexWarning: Casting complex values to real discards the imaginary part
  uvn_axes[2, :] = point_vector
C:\Users\lycea\miniconda3\envs\CudaDevelopment\Lib\site-packages\lyceanem\electromagnetics\empropagation.py:3736: ComplexWarning: Casting complex values to real discards the imaginary part
  uvn_axes[0, :] = np.cross(local_axes[2, :], point_vector) / np.linalg.norm(
C:\Users\lycea\miniconda3\envs\CudaDevelopment\Lib\site-packages\lyceanem\electromagnetics\empropagation.py:3758: ComplexWarning: Casting complex values to real discards the imaginary part
  uvn_axes[1, :] = np.cross(point_vector, uvn_axes[0, :]) / np.linalg.norm(
C:\Users\lycea\miniconda3\envs\CudaDevelopment\Lib\site-packages\numba_cuda\numba\cuda\dispatcher.py:693: NumbaPerformanceWarning: Grid size 115 will likely result in GPU under-utilization due to low occupancy.
  warn(NumbaPerformanceWarning(msg))
C:\Users\lycea\miniconda3\envs\CudaDevelopment\Lib\site-packages\numba_cuda\numba\cuda\dispatcher.py:693: NumbaPerformanceWarning: Grid size 20 will likely result in GPU under-utilization due to low occupancy.
  warn(NumbaPerformanceWarning(msg))

Plot Normalised Response

Using matplotlib, plot the results

import matplotlib.pyplot as plt

time_index = np.linspace(0, model_time * 1e9, num_samples)
time, anglegrid = np.meshgrid(time_index[:1801], angle_values - 45)
norm_max = np.nanmax(
    np.array(
        [
            np.nanmax(10 * np.log10((Ex**2) / receiver_impedence)),
            np.nanmax(10 * np.log10((Ey**2) / receiver_impedence)),
            np.nanmax(10 * np.log10((Ez**2) / receiver_impedence)),
        ]
    )
)

fig2, ax2 = plt.subplots(constrained_layout=True)
origin = "lower"
# Now make a contour plot with the levels specified,
# and with the colormap generated automatically from a list
# of colors.

levels = np.linspace(-80, 0, 41)

CS = ax2.contourf(
    anglegrid,
    time,
    10 * np.log10((Ez[:, :1801] ** 2) / receiver_impedence) - norm_max,
    levels,
    origin=origin,
    extend="both",
)
cbar = fig2.colorbar(CS)
cbar.ax.set_ylabel("Received Power (dBm)")

ax2.set_ylim(0, 30)
ax2.set_xlim(-45, 45)

ax2.set_xticks(np.linspace(-45, 45, 7))
ax2.set_yticks(np.linspace(0, 30, 16))

ax2.set_xlabel("Rotation Angle (degrees)")
ax2.set_ylabel("Time of Flight (ns)")
ax2.set_title("Received Power vs Time for rotating Plate (24GHz)")

from scipy.fft import fft, fftfreq
import scipy

xf = fftfreq(len(time_index), 1 / sampling_freq)[: len(time_index) // 2]
frequency_index = np.where(xf == model_freq)
input_signal = excitation_signal * (output_amplitude_peak)
inputfft = fft(input_signal)
input_freq = fftfreq(120, 1 / sampling_freq)[:60]
freqfuncabs = scipy.interpolate.interp1d(input_freq, np.abs(inputfft[:60]))
freqfuncangle = scipy.interpolate.interp1d(input_freq, np.angle(inputfft[:60]))
newinput = freqfuncabs(xf[frequency_index]) * np.exp(freqfuncangle(xf[frequency_index]))
Exf = fft(Ex)
Eyf = fft(Ey)
Ezf = fft(Ez)

Frequency Specific Results

The time of flight plot is useful to displaying the output of the model, giving a understanding about what is physically happening in the channel, but to get an idea of the behaviour in the frequency domain we need to use a fourier transform to move from time and voltages to frequency.

s21x = 20 * np.log10(np.abs(Exf[:, frequency_index] / newinput)).ravel()
s21y = 20 * np.log10(np.abs(Eyf[:, frequency_index] / newinput)).ravel()
s21z = 20 * np.log10(np.abs(Ezf[:, frequency_index] / newinput)).ravel()
tdangles = np.linspace(-45, 45, 91)
fig, ax = plt.subplots()
ax.plot(tdangles, s21x - np.nanmax(s21z), label="Ex")
ax.plot(tdangles, s21y - np.nanmax(s21z), label="Ey")
ax.plot(tdangles, s21z - np.nanmax(s21z), label="Ez")
plt.xlabel("$\\theta_{N}$ (degrees)")
plt.ylabel("Normalised Level (dB)")
ax.set_ylim(-60.0, 0)
ax.set_xlim(np.min(angle_values) - 45, np.max(angle_values) - 45)
ax.set_xticks(np.linspace(np.min(angle_values) - 45, np.max(angle_values) - 45, 19))
ax.set_yticks(np.linspace(-60, 0.0, 21))
legend = ax.legend(loc="upper right", shadow=True)
plt.grid()
plt.title("$S_{21}$ at 16GHz")
plt.show()