Skip to content
Snippets Groups Projects
Commit 10b44351 authored by Eric Kooistra's avatar Eric Kooistra
Browse files

Minor updates.

parent 7ab94702
No related branches found
No related tags found
No related merge requests found
Pipeline #34868 passed
Stage: simulation
Stage: synthesis
Hide whitespace changes
Inline Side-by-side
applications/lofar2/model/lofar2_station_sdp_firmware_model.ipynb
+ 132
− 34
View file @ 10b44351
%% Cell type:markdown id:82c10597 tags:
# LOFAR2.0 Station SDP Firmware quantization model
Author: Eric Kooistra, 18 May 2022
Purpose: Model the expected signal levels in the SDP firmware and clarify calculations in [1].
References:
1. https://support.astron.nl/confluence/display/L2M/L4+SDPFW+Decision%3A+LOFAR2.0+SDP+Firmware+Quantization+Model
2. Understanding digital signal processing, R.G. Lyons
%% Cell type:code id:2b477516 tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
```
%% Cell type:markdown id:c2cc6c7a tags:
## 1 SDP Parameters
%% Cell type:code id:e1b6fa12 tags:
``` python
# SDP
N_complex = 2
N_fft = 1024
N_sub = N_fft / N_complex
f_adc = 200e6 # Hz
f_sub = f_adc / N_fft
T_int = 1 # s
N_int = f_adc * T_int
N_int_sub = f_sub * T_int
print(f"N_int = {N_int:.0f}")
print(f"N_int_sub = {N_int_sub:.1f}")
```
%% Output
N_int = 200000000
N_int_sub = 195312.5
%% Cell type:code id:eb325c9c tags:
``` python
# Signal data widths and full scale
W_adc = 14
W_subband = 18
W_crosslet = 16
W_beamlet_sum = 18
W_beamlet = 8
W_statistic = 64
FS = 2**(W_adc - 1) # full scale
sigma_fs_sine = FS / np.sqrt(2)
print("FS =", FS)
print(f"sigma_fs_sine = {sigma_fs_sine:.1f}")
```
%% Output
FS = 8192
sigma_fs_sine = 5792.6
%% Cell type:code id:3e71626f tags:
``` python
# Widths of coefficients and weights
W_fir_coef = 16
W_sub_weight = 16
W_sub_weight_fraction = 13
W_sub_weight_magnitude = W_sub_weight - W_sub_weight_fraction - 1
W_bf_weight = 16
W_bf_weight_fraction = 14
W_bf_weight_magnitude = W_bf_weight - W_bf_weight_fraction -1
W_beamlet_scale = 16
W_beamlet_scale_fraction = 15
W_beamlet_scale_magnitude = W_beamlet_scale - W_beamlet_scale_fraction
Unit_sub_weight = 2**W_sub_weight_fraction
Unit_bf_weight = 2**W_bf_weight_fraction
Unit_beamlet_scale = 2**W_beamlet_scale_fraction
print("Unit_sub_weight =", Unit_sub_weight)
print("Unit_bf_weight =", Unit_bf_weight)
print("Unit_beamlet_scale =", Unit_beamlet_scale)
```
%% Output
Unit_sub_weight = 8192
Unit_bf_weight = 16384
Unit_beamlet_scale = 32768
%% Cell type:code id:def6eba7 tags:
``` python
# DFT of real sine input --> show that:
G_fft_real_input_sine = 0.5
G_fft_real_input_dc = 1.0
# . DFT size
N_points = 1024
N_bins = N_points // 2 + 1 # positive frequency bins including DC and F_s/2
N_bins = N_points // 2 + 1 # positive frequency bins including DC and f_s/2
# . select a bin
i_bin = 200 # bin index in range(N_points // 2 )
i_bin = 200 # bin index in range(N_bins)
# . time and frequency axis
f_s = f_adc # sample frequency
f_s = 1 # normalized sample frequency
T_s = 1 / f_s # sample period
T_fft = N_points * T_s # DFT period
t_axis = np.linspace(0, T_fft, N_points, endpoint=False)
f_axis = np.linspace(0, f_s, N_points, endpoint=False)
f_axis_fft = f_axis - f_s/2 # fftshift axis
f_axis_rfft = f_axis[0:N_bins] # positive frequency bins
f_bin = i_bin / N_points * f_s # bin frequency
# . create sine at bin, use cos to see DC at i_bin = 0
x = np.cos(2 * np.pi * f_bin * t_axis)
# . DFT using complex input fft()
X_fft = np.fft.fftshift(np.fft.fft(x) / N_points)
# . DFT using real input rfft()
X_rfft = np.fft.rfft(x) / N_points
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.title('DFT of real sine using fft')
plt.plot(f_axis_fft, abs(X_fft))
plt.grid()
plt.subplot(1, 2, 2)
plt.title('DFT of real sine using rfft')
plt.plot(f_axis_rfft, abs(X_rfft))
plt.grid()
print("Conclusion: G_fft_real_input_sine =", G_fft_real_input_sine)
```
%% Output
Conclusion: G_fft_real_input_sine = 0.5
%% Cell type:code id:0ec00484 tags:
``` python
# Subband filterbank (F_sub)
# . FIR filter DC gain
G_fir_dc = 0.994817
# . Signal level bit growth to accomodate processing gain of FFT
W_sub_proc = np.log2(np.sqrt(N_sub))
W_sub_gain = 4
W_sub_gain = 4 # use W_sub_gain instead of W_sub_proc
# Subband equalizer (E_sub)
subband_weight_gain = 1.0
subband_weight_phase = 0
subband_weight_re = int(subband_weight_gain * Unit_sub_weight * np.cos(subband_weight_phase))
subband_weight_im = int(subband_weight_gain * Unit_sub_weight * np.sin(subband_weight_phase))
print("subband_weight_gain =", subband_weight_gain)
print("subband_weight_phase =", subband_weight_phase)
print(f"subband_weight_re = {subband_weight_re:d}")
print(f"subband_weight_im = {subband_weight_im:d}")
print()
# . Expected factor between subband amplitude and real signal input amplitude
G_subband = G_fir_dc * G_fft_real_input_sine * 2**W_sub_gain * subband_weight_gain
print(f"G_subband = {G_fir_dc} * {G_fft_real_input_sine} * 2**{W_sub_gain} * {subband_weight_gain} = {G_subband}")
print(" . G_fir_dc =", G_fir_dc)
print(" . G_fft_real_input_sine =", G_fft_real_input_sine)
print(" . W_sub_proc =", W_sub_proc)
print(" . W_sub_gain =", W_sub_gain)
print(" . subband_weight_gain =", subband_weight_gain)
```
%% Output
subband_weight_gain = 1.0
subband_weight_phase = 0
subband_weight_re = 8192
subband_weight_im = 0
G_subband = 0.994817 * 0.5 * 2**4 * 1.0 = 7.958536
. G_fir_dc = 0.994817
. G_fft_real_input_sine = 0.5
. W_sub_gain = 4
. subband_weight_gain = 1.0
. G_fir_dc = 0.994817
. G_fft_real_input_sine = 0.5
. W_sub_proc = 4.5
. W_sub_gain = 4
. subband_weight_gain = 1.0
%% Cell type:code id:4d197368 tags:
``` python
# Digital beamformer (BF)
N_ant = 10
# Assume all N_ant use same BF weight
beamlet_weight_gain = 1.0
beamlet_weight_phase = 0
beamlet_weight_re = int(beamlet_weight_gain * Unit_bf_weight * np.cos(beamlet_weight_phase))
beamlet_weight_im = int(beamlet_weight_gain * Unit_bf_weight * np.sin(beamlet_weight_phase))
print("beamlet_weight_gain =", beamlet_weight_gain)
print("beamlet_weight_phase =", beamlet_weight_phase)
print(f"beamlet_weight_re = {beamlet_weight_re:d}")
print(f"beamlet_weight_im = {beamlet_weight_im:d}")
print()
si_types = ["coherent", "incoherent"]
for si_type in si_types:
# . BF processing gain
if si_type == "coherent":
bf_proc = N_ant
else:
bf_proc = np.sqrt(N_ant)
# . Normalize BF weights to get BF DC gain is 1.0
beamlet_weight_gain = 1 / bf_proc
# . Expected factor between beamlet amplitude and real signal input amplitude
G_beamlet_sum = N_ant * beamlet_weight_gain * G_subband
print(f"BF for {si_type} input:")
print(f" . bf_proc = {bf_proc:.2f} for N_ant = {N_ant}")
print()
```
%% Output
beamlet_weight_gain = 1.0
beamlet_weight_phase = 0
beamlet_weight_re = 16384
beamlet_weight_im = 0
BF for coherent input:
. W_bf_proc = 10.00 for N_ant = 10
. bf_proc = 10.00 for N_ant = 10
BF for incoherent input:
. W_bf_proc = 3.16 for N_ant = 10
. bf_proc = 3.16 for N_ant = 10
%% Cell type:markdown id:d942fcc6 tags:
## 2 Quantization model
%% Cell type:code id:f66c5028 tags:
``` python
# Bit
# . Each bit yields a factor 2 in voltage, so a factor 2**2 = 4 in power
P_bit = 2**2
P_bit_dB = 10 * np.log10(P_bit)
print(f"P_bit_dB = {P_bit_dB:.2f} dB")
```
%% Output
P_bit_dB = 6.02 dB
%% Cell type:code id:a9fca052 tags:
``` python
# Quantization noise
# . The quantization noise power is q**2 * 1 / 12, so the standard deviation
# of the quantization noise is q * sqrt(1 / 12) < q = one LSbit
# . The quantization noise power is at a level of -10.79 dB or -1.8 bit.
# . The 0 dB power level corresponds to the power of one LSbit, so q**2
P_quant = 1 / 12 # for W >> 1 [2]
P_quant_dB = 10 * np.log10(P_quant)
sigma_quant = np.sqrt(P_quant)
print()
print(f"P_quant = {P_quant:.6f}")
print(f"P_quant_dB = {P_quant_dB:.2f} dB = {P_quant_dB / P_bit_dB:.1f} bit")
print(f"sigma_quant = {sigma_quant:.2f} q")
```
%% Output
P_quant = 0.083333
P_quant_dB = -10.79 dB = -1.8 bit
sigma_quant = 0.29 q
%% Cell type:code id:d9972b6b tags:
``` python
# System noise
n = np.arange(1,9)
sigma_sys = n
sigma_sys = n # = n * q, so n LSbits
P_sys = sigma_sys**2
P_tot = P_sys + P_quant
sigma_tot = np.sqrt(P_tot)
plt.figure()
plt.plot(n, (P_tot / P_sys - 1) * 100)
plt.title("Increase in total noise power due to sampling")
plt.xlabel("sigma_sys [q]")
plt.ylabel("[%]")
plt.grid()
```
%% Output
%% Cell type:code id:be2d952f tags:
``` python
# FS sine
P_fs_sine = FS**2 / 2
P_fs_sine_dB = 10 * np.log10(P_fs_sine)
print("W_adc = {W_adc} bits")
print("FS =", FS)
print(f"sigma_fs_sine = {sigma_fs_sine:.1f} q")
print(f"P_fs_sine_dB = {P_fs_sine_dB:.2f} dB = {P_fs_sine_dB / P_bit_dB:.1f} bit")
```
%% Output
W_adc = {W_adc} bits
FS = 8192
sigma_fs_sine = 5792.6 q
P_fs_sine_dB = 75.26 dB = 12.5 bit
%% Cell type:code id:a9e7fabc tags:
``` python
# SNR
SNR = P_fs_sine / P_quant
SNR_dB = 10 * np.log10(SNR)
print()
print(f"SNR_dB = P_fs_sine_dB - P_quant_dB = {P_fs_sine_dB:.2f} - {P_quant_dB:.2f} = {SNR_dB:.2f} dB")
```
%% Output
SNR_dB = P_fs_sine_dB - P_quant_dB = 75.26 - -10.79 = 86.05 dB
%% Cell type:code id:92852a53 tags:
``` python
# -50 dbFS level
Power_50dBFS = P_fs_sine_dB - 50
sigma_50dBFS = 10**(Power_50dBFS / 20)
ampl_50dBFS = sigma_50dBFS * np.sqrt(2)
print(f"Power at -50dBFS = {Power_50dBFS:.2f} dB, so sigma = {sigma_50dBFS:.1f} q, ampl = {ampl_50dBFS:.1f} q")
# Assume sigma = 16 q
sigma = 16
Power = sigma**2
Power_dB = 10 * np.log10(Power)
print()
print(f"sigma = {sigma:.0f} q corresponds to:")
print(f" . Power = {Power_dB:.2f} dB, so at {Power_dB - P_fs_sine_dB:.1f} dBFS")
print(f" . Range 3 sigma = +-{3 * sigma:.0f} q")
print(f" . Sine with amplitude A = sigma * sqrt(2) = {np.sqrt(2) * sigma:.1f} q")
```
%% Output
Power at -50dBFS = 25.26 dB, so sigma = 18.3 q, ampl = 25.9 q
sigma = 16 q corresponds to:
. Power = 24.08 dB, so at -51.2 dBFS
. Range 3 sigma = +-48 q
. Sine with amplitude A = sigma * sqrt(2) = 22.6 q
%% Cell type:markdown id:77bb14cc tags:
## 3 Expected signal levels in the SDP FW
%% Cell type:code id:a04af043 tags:
``` python
# ADC power statistic (AST)
sigma = sigma_fs_sine
sigma_bits = np.log2(sigma)
P_ast = (sigma)**2 * N_int
print(f"ADC sigma = {sigma:6.1f} q = {sigma_bits:4.1f} bits: P_ast = {P_ast:e}, uses {np.log2(P_ast):.1f} bits, is 0 dBFS = FS sine")
sigma = sigma_50dBFS
sigma_bits = np.log2(sigma)
P_ast = (sigma)**2 * N_int
print(f"ADC sigma = {sigma:6.1f} q = {sigma_bits:4.1f} bits: P_ast = {P_ast:e}, uses {np.log2(P_ast):.1f} bits, is -50dBFS")
sigma = 16
sigma_bits = np.log2(sigma)
P_ast = (sigma)**2 * N_int
print(f"ADC sigma = {sigma:6.1f} q = {sigma_bits:4.1f} bits: P_ast = {P_ast:e}, uses {np.log2(P_ast):.1f} bits")
```
%% Output
ADC sigma = 5792.6 q = 12.5 bits: P_ast = 6.710886e+15, uses 52.6 bits, is 0 dBFS = FS sine
ADC sigma = 18.3 q = 4.2 bits: P_ast = 6.710886e+10, uses 36.0 bits, is -50dBFS
ADC sigma = 16.0 q = 4.0 bits: P_ast = 5.120000e+10, uses 35.6 bits
%% Cell type:code id:0b2ac36c tags:
``` python
# Subband filterbank (F_sub)
ampl_sub_fs = FS * G_subband # subband amplitude for FS signal input sine
SST_fs = ampl_sub_fs**2 * N_int_sub
ampl_sub_50dBFS = ampl_50dBFS * G_subband # subband amplitude -50dBFS signal input sine
SST_50dBFS = ampl_sub_50dBFS**2 * N_int_sub
ampl_si_16q = 16 # [q], so 16 q is 4 bits amplitude
ampl_sub_16q = ampl_si_16q * G_subband # subband amplitude for signal input sine with A = 16 q
SST_ampl_16q = ampl_sub_16q**2 * N_int_sub
sigma_si_16q = 16 * np.sqrt(2) # [q], so A = 16 * sqrt(2) q for sigma = 16 q
sigma_sub_16q = sigma_si_16q * G_subband # subband sigma for arbitrary signal input with sigma = 16 q
SST_sigma_16q = sigma_sub_16q**2 * N_int_sub
print("Signal input level --> Expected subband level and SST level:")
print()
print(" ampl_si ampl_sub #bits SST #bits")
print(f"{FS:8.1f} {ampl_sub_fs:10.1f} {np.log2(ampl_sub_fs):8.4f} {SST_fs:e} {np.log2(SST_fs):.1f}")
print(f"{ampl_50dBFS:8.1f} {ampl_sub_50dBFS:10.1f} {np.log2(ampl_sub_50dBFS):8.4f} {SST_50dBFS:e} {np.log2(SST_50dBFS):.1f}")
print(f"{ampl_si_16q:8.1f} {ampl_sub_16q:10.1f} {np.log2(ampl_sub_16q):8.4f} {SST_ampl_16q:e} {np.log2(SST_ampl_16q):.1f}")
print()
print("sigma_si sigma_sub #bits SST #bits")
print(f"{sigma_si_16q:8.1f} {sigma_sub_16q:10.1f} {np.log2(sigma_sub_16q):8.4f} {SST_sigma_16q:e} {np.log2(SST_sigma_16q):.1f}")
```
%% Output
Signal input level --> Expected subband level and SST level:
ampl_si ampl_sub #bits SST #bits
8192.0 65196.3 15.9925 8.301877e+14 49.6
25.9 206.2 7.6877 8.301877e+09 33.0
16.0 127.3 6.9925 3.166915e+09 31.6
sigma_si sigma_sub #bits SST #bits
22.6 180.1 7.4925 6.333830e+09 32.6
%% Cell type:code id:f0b09a83 tags:
``` python
# Digital beamformer (BF)
# . is a coherent beamformer = voltage beamformer
# . uses BF weights to form beamlets from sum of weighted subbands
```
%% Cell type:code id:8713e865 tags:
``` python
```
......
applications/lofar2/model/signal_statistics.ipynb
+ 77
− 87
View file @ 10b44351
Source diff could not be displayed: it is too large. Options to address this: view the blob.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment