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

Separated statistics modelling from sdp firmware model.

parent 8e949cbe
No related branches found
No related tags found
No related merge requests found
Expand all Show whitespace changes
Inline Side-by-side
applications/lofar2/model/lofar2_station_sdp_firmware_model.ipynb
+ 225
− 407
View file @ 3d47bbf8
%% 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 math
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_ant = 96
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 / math.sqrt(2)
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
# . select a bin
i_bin = 200 # bin index in range(N_points // 2 )
# . 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
# Signal level bit growth to accomodate processing gain
W_sub_proc = math.log2(math.sqrt(N_sub))
# 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_bf_proc = math.log2(math.sqrt(N_ant))
print("W_sub_proc =", W_sub_proc)
print(f"W_bf_proc = {W_bf_proc:.2f} for N_ant = {N_ant}")
# 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_gain =", W_sub_gain)
print(" . subband_weight_gain =", subband_weight_gain)
```
%% Output
W_sub_proc = 4.5
W_bf_proc = 3.29 for N_ant = 96
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
%% 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 for incoherent input:
. W_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
P_bit = 2**2
P_bit_dB = 10 * math.log10(P_bit)
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
P_quant = 1 / 12 # for W >> 1 [2]
P_quant_dB = 10 * math.log10(P_quant)
sigma_quant = math.sqrt(P_quant)
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
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 * math.log10(P_fs_sine)
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 * math.log10(SNR)
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")
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 * math.log10(Power)
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) = {math.sqrt(2) * sigma:.1f} 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
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
# 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
```
%% Cell type:code id:a656367f tags:
``` python
```
%% Cell type:markdown id:27d0fe5a tags:
## 4 Signal statistics
%% Cell type:markdown id:4ddef2d8 tags:
### 4.1 Statistics basics:
* dc = mean # direct current
* sigma = std # standard deviation
* var = std**2 # variance
* mean power = var + mean**2
* rms = sqrt(mean power) = sqrt(var + mean**2)
Coherent and incoherent signals. With S signals, the std of their sum:
# 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
* increases by S for coherent signals
* increases by sqrt(S) for incoherent signals
%% Cell type:code id:9c55fb7b tags:
``` python
def rms(arr):
"""Root mean square of values in arr
rms = sqrt(mean powers) = sqrt(std**2 + mean**2)
The rms() also works for complex input thanks to using np.abs().
"""
return np.sqrt(np.mean(np.abs(arr)**2.0))
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}")
```
%% Cell type:code id:74edfe32 tags:
%% Cell type:code id:f0b09a83 tags:
``` python
N_samples = 10000
sigma = 0.5
var = sigma**2
dc = -0.2
# Signal input voltages
si = np.random.randn(N_samples)
si *= sigma / np.std(si) # apply requested sigma
si += dc # add offset
print(f"mean(si) = {np.mean(si):.6f}, expected {dc}")
print(f"std(si) = {np.std(si):.6f}, expected {sigma}")
print(f"rms(si) = {rms(si):.6f}, expected {np.sqrt(var + dc**2):.6f}")
```
%% Output
mean(si) = -0.204032, expected -0.2
std(si) = 0.500000, expected 0.5
rms(si) = 0.540027, expected 0.538516
%% Cell type:markdown id:17d333f1 tags:
### 4.2 Beamforming
In the beamformer the weak signal in the beamlet adds coherently and the sky
signals from other directions add incoherently. Hence the SNR of the weak
signal in the BF output improves by factor S/sqrt(S) = sqrt(S).
%% Cell type:code id:89845ec3 tags:
``` python
# Signal inputs
# . Warning: do use not too high values for N_samples and S_max to avoid too long compute time.
N_samples = 1000
S_max = 200
S_arr = np.arange(1, S_max + 1)
sigma_weak = 0.05
sigma_other = 0.5
si_weak = np.random.randn(N_samples)
si_weak *= sigma_weak / np.std(si_weak)
# Beamformer sum(S)
bf_weak_std_arr = []
bf_other_std_arr = []
bf_sys_std_arr = []
bf_SNR_dB_arr = []
for S in S_arr:
# The weak signal in the beamlet adds coherently for all signal inputs
bf_weak = S * si_weak
bf_weak_std = np.std(bf_weak)
bf_weak_std_arr.append(bf_weak_std)
# The other signals from other directions add incoherently
bf_other = np.zeros(N_samples)
for si in range(1, S + 1):
si_other = np.random.randn(N_samples)
si_other *= sigma_other / np.std(si_other)
bf_other += si_other
bf_other_std = np.std(bf_other)
bf_other_std_arr.append(bf_other_std)
# Total BF output
bf_sys_std = np.std(bf_weak + bf_other)
bf_sys_std_arr.append(bf_sys_std)
SNR_dB = 20 * np.log10(bf_weak_std / bf_other_std)
bf_SNR_dB_arr.append(SNR_dB)
plt.figure(1)
plt.plot(S_arr, bf_weak_std_arr, 'g', S_arr, bf_other_std_arr, 'b', S_arr, bf_sys_std_arr, 'r')
plt.title("Beamformer")
plt.xlabel("Number of signal inputs")
plt.ylabel("BF std")
plt.legend(['bf_weak', 'bf_other', 'bf_sys'])
plt.grid()
plt.figure(2)
plt.plot(S_arr, bf_SNR_dB_arr, 'r')
plt.title("Beamformer")
plt.xlabel("Number of signal inputs")
plt.ylabel("SNR [dB]")
plt.grid()
```
%% Output
%% Cell type:markdown id:71aa6647 tags:
**Conclusion:**
The beamformer improves the SNR of the weak signal by factor sqrt(S). For most very weak astronimical signals this is not enough to make them appear above the system noise, so then additional beamforming is needed or integration in time using a correlator.
%% Cell type:markdown id:84b8930c tags:
### 4.3 Correlation
%% Cell type:code id:470fd269 tags:
``` python
# Signal inputs
N_steps = 1000
N_min = 100
N_max = 100000
n_incr = (N_max / N_min)**(1 / N_steps)
N_arr = []
for s in range(N_steps + 1):
n = int(N_min * n_incr**s)
N_arr.append(n)
sigma_weak = 0.01
sigma_other = 0.5
SNR_dB = 20 * np.log10(sigma_weak / sigma_other)
print(f"SNR input = {SNR_dB:.3f} dB")
# Correlator mean(A * B)
cor_weak_arr = []
cor_other_arr = []
cor_sys_arr = []
cor_SNR_dB_arr = []
for N in N_arr:
si_weak = np.random.randn(N)
si_weak *= sigma_weak / np.std(si_weak)
# Signal input A
sA_other = np.random.randn(N)
sA_other *= sigma_other / np.std(sA_other)
sA_sys = sA_other + si_weak
# Signal input B
sB_other = np.random.randn(N)
sB_other *= sigma_other / np.std(sB_other)
sB_sys = sB_other + si_weak
# Correlate A and B
cor_weak = np.mean(si_weak * si_weak)
cor_weak_arr.append(cor_weak)
cor_other = np.mean(sA_other * sB_other)
cor_other_arr.append(cor_other)
cor_sys = np.mean(sA_sys * sB_sys)
cor_sys_arr.append(cor_sys)
#print(f"{N}, {cor_weak:9.6f}, {cor_other:9.6f}, {cor_sys:9.6f}")
SNR_dB = 10 * np.log10(np.abs(cor_weak / cor_other))
cor_SNR_dB_arr.append(SNR_dB)
#print(f"{N}, SNR output = {SNR_dB:.0f} dB")
plt.figure(1)
plt.plot(N_arr, cor_weak_arr, 'g', N_arr, cor_other_arr, 'b', N_arr, cor_sys_arr, 'r')
plt.title("Correlator")
plt.xlabel("Number of samples")
plt.ylabel("Cross power")
plt.legend(['cor_weak', 'cor_other', 'cor_sys'])
plt.grid()
plt.figure(2)
plt.plot(N_arr, cor_SNR_dB_arr, 'r')
plt.title("Correlator")
plt.xlabel("Number of samples")
plt.ylabel("SNR [dB]")
plt.grid()
# Digital beamformer (BF)
# . is a coherent beamformer = voltage beamformer
# . uses BF weights to form beamlets from sum of weighted subbands
```
%% Output
SNR input = -33.979 dB
%% Cell type:code id:8713e865 tags:
``` python
```
......
applications/lofar2/model/signal_statistics.ipynb 0 → 100644
+ 424
− 0
View file @ 3d47bbf8
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment