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

First draft.

parent 5e9d9fb0
No related branches found
No related tags found
No related merge requests found
Pipeline #30793 passed
Stage: simulation
Stage: synthesis
Hide whitespace changes
Inline Side-by-side
applications/lofar2/model/lofar2_station_sdp_firmware_model.ipynb 0 → 100644
+ 727
− 0
View file @ 901b19ba
%% 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}")
```
%% Output
N_int = 200000000
%% 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)
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:0ec00484 tags:
``` python
# Signal level bit growth to accomodate processing gain
W_sub_proc = math.log2(math.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}")
```
%% Output
W_sub_proc = 4.5
W_bf_proc = 3.29 for N_ant = 96
%% 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)
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)
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)
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)
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)
print(f"Power at -50dBFS = {Power_50dBFS:.2f} dB, so sigma = {sigma_50dBFS:.1f} q")
# Assume sigma = 16 q
sigma = 16
Power = sigma**2
Power_dB = 10 * math.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")
```
%% 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
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:
* 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))
```
%% Cell type:code id:74edfe32 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()
```
%% Output
SNR input = -33.979 dB
%% Cell type:code id:8713e865 tags:
``` python
```
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