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

Updates.

parent ce96e0f1
Branches
No related tags found
No related merge requests found
Pipeline #35725 passed
Stage: simulation
Stage: synthesis
Expand all Show whitespace changes
Inline Side-by-side
applications/lofar2/model/lofar2_station_sdp_firmware_model.ipynb
+ 80
− 40
View file @ 2aefd0ba
%% Cell type:markdown id:82c10597 tags:
# LOFAR2.0 Station SDP Firmware quantization model
Author: Eric Kooistra, Aug 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
# General
N_complex = 2
N_sidebands = 2
# SDP
N_fft = 1024 # number of time points, number of frequency bins
N_sub = N_fft / N_sidebands # number of subbands, DC and positive frequeny bins
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:0ec00484 tags:
``` python
# Gain factor between subband level and signal input level in the subband filterbank (F_sub)
# . for coherent WG sine based on amplitude
# . for white noise based on std
# . FIR filter DC gain
G_fir_dc = 0.994817 # actual DC gain of FIR filter in LOFAR
G_fir_dc = 1
# DFT gain for real input dc, sine and noise
G_fft_real_input_dc = 1 # coherent, one bin
G_fft_real_input_sine = 1 / N_sidebands # coherent, so proportional to 1/N
G_fft_real_input_noise = 1 / np.sqrt(N_fft) # incoherent, so proportional to 1/sqrt(N)
# . Signal level bit growth to accomodate processing gain of FFT
W_sub_proc = np.log2(np.sqrt(N_sub))
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(f"subband_weight_gain = {subband_weight_gain} (Unit_sub_weight = {Unit_sub_weight})")
print(f"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 from real signal input amplitude to subband amplitude
G_subband_ampl = G_fir_dc * G_fft_real_input_sine * 2**W_sub_gain * subband_weight_gain
print("Coherent WG sine input:")
print(f" G_subband_ampl = {G_fir_dc} * {G_fft_real_input_sine} * 2**{W_sub_gain} * {subband_weight_gain} \
= {G_subband_ampl} = {np.log2(G_subband_ampl):.2f} bits")
print(" . G_fir_dc =", G_fir_dc)
print(" . G_fft_real_input_sine =", G_fft_real_input_sine, "=", 1 / N_sidebands)
print(" . W_sub_gain =", W_sub_gain)
print(" . subband_weight_gain =", subband_weight_gain)
print()
# Expected factor from real signal input white noise sigma to subband amplitude
G_subband_sigma = G_fir_dc * G_fft_real_input_noise * 2**W_sub_gain * subband_weight_gain
print("Incoherent white noise input:")
print(f" G_subband_sigma = {G_fir_dc} * {G_fft_real_input_noise} * 2**{W_sub_gain} * {subband_weight_gain} \
= {G_subband_sigma} = {np.log2(G_subband_sigma):.2f} bits")
print(" . G_fir_dc =", G_fir_dc)
print(" . G_fft_real_input_noise =", G_fft_real_input_noise, "=", 1 / np.sqrt(N_fft))
print(" . W_sub_gain =", W_sub_gain)
print(" . subband_weight_gain =", subband_weight_gain)
print()
```
%% Output
subband_weight_gain = 1.0 (Unit_sub_weight = 8192)
subband_weight_phase = 0
subband_weight_re = 8192
subband_weight_im = 0
Coherent WG sine input:
G_subband_ampl = 1 * 0.5 * 2**4 * 1.0 = 8.0 = 3.00 bits
. G_fir_dc = 1
. G_fft_real_input_sine = 0.5 = 0.5
. W_sub_gain = 4
. subband_weight_gain = 1.0
Incoherent white noise input:
G_subband_sigma = 1 * 0.03125 * 2**4 * 1.0 = 0.5 = -1.00 bits
. G_fir_dc = 1
. G_fft_real_input_noise = 0.03125 = 0.03125
. W_sub_gain = 4
. subband_weight_gain = 1.0
%% Cell type:code id:ac73d7e3 tags:
``` python
# Gain factor G_beamlet_sum between beamlet and signal input in the digital beamformer (BF)
# . coherent input is same signal (e.g. sine or noise) on all inputs
# . incoherent input is different signal (noise) on all inputs
# . 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("Same BF weight for all inputs:")
print(f". beamlet_weight_gain = {beamlet_weight_gain} (Unit_bf_weight = {Unit_bf_weight})")
print(f". 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()
# Use range of N_ant for number of signal inputs in the BF
N_ant_arr = np.array([1, 12, 24, 48, 96])
print(f"N_ant_arr = {N_ant_arr}")
print()
# BF processing gain for N_ant coherent inputs and for N_ant incoherent inputs
bf_proc_coh = N_ant_arr
bf_proc_coh_bits = np.log2(bf_proc_coh)
bf_proc_incoh = np.sqrt(N_ant_arr)
bf_proc_incoh_bits = np.log2(bf_proc_incoh)
print("Gain from subband to beamlet:")
print()
for ni, na in enumerate(N_ant_arr):
print(f" N_ant = {na:2d} : bf_proc_coh = {bf_proc_coh[ni]:5.2f} = {np.log2(bf_proc_coh[ni]):.1f} bits")
print()
for ni, na in enumerate(N_ant_arr):
print(f" N_ant = {na:2d} : bf_proc_incoh = {bf_proc_incoh[ni]:5.2f} = {np.log2(bf_proc_incoh[ni]):.1f} bits")
print()
# Expected factor from real signal input amplitude to beamlet amplitude
G_beamlet_sum_ampl = G_subband_ampl * bf_proc_coh * beamlet_weight_gain
# Expected factor from real signal input sigma to beamlet sigma
G_beamlet_sum_sigma = G_subband_sigma * bf_proc_incoh * beamlet_weight_gain
print("Gain from signal input to beamlet:")
print()
for ni, na in enumerate(N_ant_arr):
print(f" N_ant = {na:2d} : G_beamlet_sum_ampl = {G_subband_ampl:.2f} * " \
f"{bf_proc_coh[ni]:5.2f} * {beamlet_weight_gain} " \
f"= {G_beamlet_sum_ampl[ni]:7.2f} = {np.log2(G_beamlet_sum_ampl[ni]):.1f} bits")
print()
for ni, na in enumerate(N_ant_arr):
print(f" N_ant = {na:2d} : G_beamlet_sum_sigma = {G_subband_sigma:.2f} * " \
f"{bf_proc_incoh[ni]:5.2f} * {beamlet_weight_gain} " \
f"= {G_beamlet_sum_sigma[ni]:6.2f} = {np.log2(G_beamlet_sum_sigma[ni]):.1f} bits")
print()
```
%% Output
Same BF weight for all inputs:
. beamlet_weight_gain = 1.0 (Unit_bf_weight = 16384)
. beamlet_weight_phase = 0
. beamlet_weight_re = 16384
. beamlet_weight_im = 0
N_ant_arr = [ 1 12 24 48 96]
Gain from subband to beamlet:
N_ant = 1 : bf_proc_coh = 1.00 = 0.0 bits
N_ant = 12 : bf_proc_coh = 12.00 = 3.6 bits
N_ant = 24 : bf_proc_coh = 24.00 = 4.6 bits
N_ant = 48 : bf_proc_coh = 48.00 = 5.6 bits
N_ant = 96 : bf_proc_coh = 96.00 = 6.6 bits
N_ant = 1 : bf_proc_incoh = 1.00 = 0.0 bits
N_ant = 12 : bf_proc_incoh = 3.46 = 1.8 bits
N_ant = 24 : bf_proc_incoh = 4.90 = 2.3 bits
N_ant = 48 : bf_proc_incoh = 6.93 = 2.8 bits
N_ant = 96 : bf_proc_incoh = 9.80 = 3.3 bits
Gain from signal input to beamlet:
N_ant = 1 : G_beamlet_sum_ampl = 8.00 * 1.00 * 1.0 = 8.00 = 3.0 bits
N_ant = 12 : G_beamlet_sum_ampl = 8.00 * 12.00 * 1.0 = 96.00 = 6.6 bits
N_ant = 24 : G_beamlet_sum_ampl = 8.00 * 24.00 * 1.0 = 192.00 = 7.6 bits
N_ant = 48 : G_beamlet_sum_ampl = 8.00 * 48.00 * 1.0 = 384.00 = 8.6 bits
N_ant = 96 : G_beamlet_sum_ampl = 8.00 * 96.00 * 1.0 = 768.00 = 9.6 bits
N_ant = 1 : G_beamlet_sum_sigma = 0.50 * 1.00 * 1.0 = 0.50 = -1.0 bits
N_ant = 12 : G_beamlet_sum_sigma = 0.50 * 3.46 * 1.0 = 1.73 = 0.8 bits
N_ant = 24 : G_beamlet_sum_sigma = 0.50 * 4.90 * 1.0 = 2.45 = 1.3 bits
N_ant = 48 : G_beamlet_sum_sigma = 0.50 * 6.93 * 1.0 = 3.46 = 1.8 bits
N_ant = 96 : G_beamlet_sum_sigma = 0.50 * 9.80 * 1.0 = 4.90 = 2.3 bits
%% Cell type:code id:98f1917e tags:
``` python
# Maximum coherent signal input amplitude
si_ampl_max = 2**(W_beamlet_sum - 1) / G_beamlet_sum_ampl
for ni, na in enumerate(N_ant_arr):
print(f"N_ant = {na:2d} : si_ampl_max = {si_ampl_max[ni] / FS:f} " \
f"= {si_ampl_max[ni]:6.0f} = {np.log2(si_ampl_max[ni]):5.1f} bits")
print()
```
%% Output
N_ant = 1 : si_ampl_max = 2.000000 = 16384 = 14.0 bits
N_ant = 12 : si_ampl_max = 0.166667 = 1365 = 10.4 bits
N_ant = 24 : si_ampl_max = 0.083333 = 683 = 9.4 bits
N_ant = 48 : si_ampl_max = 0.041667 = 341 = 8.4 bits
N_ant = 96 : si_ampl_max = 0.020833 = 171 = 7.4 bits
%% 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 or 0 bit 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
# Impact of quantization on the system noise
# The quantization noise has sigma_quant = 0.29 q, this increases the system noise.
# The system noise has sigma_sys = n * q. For n = 2 the quantization increases the
# total power by 2% (so sigma_sys increase by sqrt(2 %) is about 1 %).
n = np.arange(1,9)
sigma_sys = n # = n * q, so sigma of 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 quantization")
plt.xlabel("sigma_sys [q]")
plt.ylabel("[%]")
plt.grid()
```
%% Output
%% Cell type:code id:be2d952f tags:
``` python
# Full scale (FS) sine
P_fs_sine = FS**2 / 2
P_fs_sine_dB = 10 * np.log10(P_fs_sine)
print(f"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 = 14 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
# Signal to noise ratio (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
# Signal level relative to FS sine
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 corresponds to:")
print(f" . sigma = {sigma_50dBFS:.1f} q")
print(f" . Noise range 3 sigma = +-{3 * sigma_50dBFS:.0f} q")
print(f" . Sine with amplitude A = = sigma * sqrt(2) = {ampl_50dBFS:.1f} q")
# Assume signal with sigma = 16 q
sigma_16q = 16
power_16q = sigma_16q**2
power_16q_dB = 10 * np.log10(power_16q)
dBFS_16q = power_16q_dB - P_fs_sine_dB
print()
print(f"sigma = {sigma_16q:.0f} q corresponds to:")
print(f" . Power = {power_16q_dB:.2f} dB, so at {dBFS_16q:.1f} dBFS")
print(f" . Noise range 3 sigma = +-{3 * sigma_16q:.0f} q")
print(f" . Sine with amplitude A = sigma * sqrt(2) = {np.sqrt(2) * sigma_16q:.1f} q")
```
%% Output
Power at -50dBFS = 25.26 dB corresponds to:
. sigma = 18.3 q
. Noise range 3 sigma = +-55 q
. Sine with amplitude A = = sigma * sqrt(2) = 25.9 q
sigma = 16 q corresponds to:
. Power = 24.08 dB, so at -51.2 dBFS
. Noise 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:markdown id:f7fff7a0 tags:
## 3.1 Signal input power and DC level
%% Cell type:code id:a04af043 tags:
``` python
# Signal input power statistic for ADC / WG (AST)
si_sigma = sigma_fs_sine
si_sigma_bits = np.log2(si_sigma)
P_ast = (si_sigma)**2 * N_int
print(f"SI sigma = {si_sigma:6.1f} q = {si_sigma_bits:4.1f} bits: P_ast = {P_ast:e},\
uses {np.log2(P_ast):.1f} bits, is 0 dBFS = FS sine")
si_sigma = sigma_50dBFS
si_sigma_bits = np.log2(si_sigma)
P_ast = (si_sigma)**2 * N_int
print(f"SI sigma = {si_sigma:6.1f} q = {si_sigma_bits:4.1f} bits: P_ast = {P_ast:e},\
uses {np.log2(P_ast):.1f} bits, is -50 dBFS")
si_sigma = sigma_16q
si_sigma_bits = np.log2(si_sigma)
P_ast = (si_sigma)**2 * N_int
print(f"SI sigma = {si_sigma:6.1f} q = {si_sigma_bits:4.1f} bits: P_ast = {P_ast:e},\
uses {np.log2(P_ast):.1f} bits, is {dBFS_16q:.1f} dBFS")
```
%% Output
SI sigma = 5792.6 q = 12.5 bits: P_ast = 6.710886e+15,uses 52.6 bits, is 0 dBFS = FS sine
SI sigma = 18.3 q = 4.2 bits: P_ast = 6.710886e+10,uses 36.0 bits, is -50 dBFS
SI sigma = 16.0 q = 4.0 bits: P_ast = 5.120000e+10,uses 35.6 bits, is -51.2 dBFS
%% Cell type:markdown id:7ce94d23 tags:
From measured P_ast and DC_ast to signal input sigma in q units:
* si_rms = sqrt(P_ast / N_int)
* si_mean = DC_ast / N_int
* si_sigma = sqrt(si_rms^2 - si_mean^2)
%% Cell type:markdown id:5cb555b6 tags:
## 3.2 Subband statistics (SST)
%% Cell type:markdown id:f842d856 tags:
For a complex signal (like subbands and beamlets):
* power complex = power real + power imag = (std real)^2 + (std imag)^2
* power real = power imag = power complex / 2
* std real = std imag = std complex / sqrt(2)
* std complex = sqrt(power complex)
* ampl real = ampl imag = std complex = std real * sqrt(2) = std imag * sqrt(2)
%% Cell type:code id:5ba30659 tags:
``` python
# Subband level and SST level for coherent (WG sine) input
# . use ampl real = ampl imag = std complex = sqrt(power complex) to calculate SST for sine input
si_ampl_fs = FS
si_ampl_fs_bits = np.log2(si_ampl_fs)
si_sigma_fs = si_ampl_fs / np.sqrt(2)
si_sigma_fs_bits = np.log2(si_sigma)
sub_ampl_fs = si_ampl_fs * G_subband_ampl # subband amplitude for FS signal input sine
sub_ampl_fs_bits = np.log2(sub_ampl_fs)
SST_fs = sub_ampl_fs**2 * N_int_sub
SST_fs_dB = 10 * np.log10(SST_fs)
SST_fs_bits = np.log2(SST_fs)
si_ampl_fs4 = FS / 4
si_ampl_fs4_bits = np.log2(si_ampl_fs4)
si_sigma_fs4 = si_ampl_fs4 / np.sqrt(2)
si_sigma_fs4_bits = np.log2(si_sigma)
sub_ampl_fs4 = si_ampl_fs4 * G_subband_ampl # subband amplitude for FS signal input sine
sub_ampl_fs4_bits = np.log2(sub_ampl_fs4)
SST_fs4 = sub_ampl_fs4**2 * N_int_sub
SST_fs4_dB = 10 * np.log10(SST_fs4)
SST_fs4_bits = np.log2(SST_fs4)
si_ampl_50dBFS = ampl_50dBFS
si_ampl_50dBFS_bits = np.log2(si_ampl_50dBFS)
si_sigma_50dBFS = sigma_50dBFS
si_sigma_50dBFS_bits = np.log2(si_sigma_50dBFS)
sub_ampl_50dBFS = si_ampl_50dBFS * G_subband_ampl # subband amplitude -50dBFS signal input sine
sub_ampl_50dBFS_bits = np.log2(sub_ampl_50dBFS)
SST_50dBFS = sub_ampl_50dBFS**2 * N_int_sub
SST_50dBFS_dB = 10 * np.log10(SST_50dBFS)
SST_50dBFS_bits = np.log2(SST_50dBFS)
si_ampl_s16q = sigma_16q * np.sqrt(2)
si_ampl_s16q_bits = np.log2(si_ampl_s16q)
si_sigma_s16q = sigma_16q
si_sigma_s16q_bits = np.log2(sigma_16q) # = 16
sub_ampl_s16q = si_ampl_s16q * G_subband_ampl # subband amplitude for signal input sine with sigma = 16 q
sub_ampl_s16q_bits = np.log2(sub_ampl_s16q)
SST_s16q = sub_ampl_s16q**2 * N_int_sub
SST_s16q_dB = 10 * np.log10(SST_s16q)
SST_s16q_bits = np.log2(SST_s16q)
print("Coherent (WG sine) signal input level --> Expected subband level and SST level:")
print()
print(" si_ampl si_sigma sub_sigma = SST")
print(" sub_ampl")
print(" value #bits value #bits value #bits value dB #bits")
print(f"{si_ampl_fs:9.1f} {si_ampl_fs_bits:6.1f} " \
f"{si_sigma_fs:9.1f} {si_sigma_fs_bits:6.1f} " \
f"{sub_ampl_fs:12.1f} {sub_ampl_fs_bits:6.1f} " \
f"{SST_fs:15e} {SST_fs_dB:6.1f} {SST_fs_bits:6.1f}, "
f"at 0 dBFS (= FS sine)")
print(f"{si_ampl_fs4:9.1f} {si_ampl_fs4_bits:6.1f} " \
f"{si_sigma_fs4:9.1f} {si_sigma_fs4_bits:6.1f} " \
f"{sub_ampl_fs4:12.1f} {sub_ampl_fs4_bits:6.1f} " \
f"{SST_fs4:15e} {SST_fs4_dB:6.1f} {SST_fs4_bits:6.1f}, "
f"at {20*np.log10(1 / 4**2):.1f} dBFS (= FS / 4)")
print(f"{si_ampl_50dBFS:9.1f} {si_ampl_50dBFS_bits:6.1f} " \
f"{si_sigma_50dBFS:9.1f} {si_sigma_50dBFS_bits:6.1f} " \
f"{sub_ampl_50dBFS:12.1f} {sub_ampl_50dBFS_bits:6.1f} " \
f"{SST_50dBFS:15e} {SST_50dBFS_dB:6.1f} {SST_50dBFS_bits:6.1f}, "
f"at -50 dBFS (= FS / {10**(50/20):.0f})")
print(f"{si_ampl_s16q:9.1f} {si_ampl_s16q_bits:6.1f} " \
f"{si_sigma_s16q:9.1f} {si_sigma_s16q_bits:6.1f} " \
f"{sub_ampl_s16q:12.1f} {sub_ampl_s16q_bits:6.1f} "\
f"{SST_s16q:15e} {SST_s16q_dB:6.1f} {SST_s16q_bits:6.1f}, "\
f"at {dBFS_16q:.1f} dBFS (= FS / {10**(-dBFS_16q/20):.0f})")
```
%% Output
Coherent (WG sine) signal input level --> Expected subband level and SST level:
si_ampl si_sigma sub_sigma = SST
sub_ampl
value #bits value #bits value #bits value dB #bits
8192.0 13.0 5792.6 4.0 65536.0 16.0 8.388608e+14 149.2 49.6, at 0 dBFS (= FS sine)
2048.0 11.0 1448.2 4.0 16384.0 14.0 5.242880e+13 137.2 45.6, at -24.1 dBFS (= FS / 4)
25.9 4.7 18.3 4.2 207.2 7.7 8.388608e+09 99.2 33.0, at -50 dBFS (= FS / 316)
22.6 4.5 16.0 4.0 181.0 7.5 6.400000e+09 98.1 32.6, at -51.2 dBFS (= FS / 362)
%% Cell type:code id:5ec1330a tags:
``` python
# Subband level and SST level for incoherent white noise input
# . the signal input power is equally distributed over all N_sub subbands
# . use std complex to calculate SST for noise input
# . use std real = std imag = std complex / sqrt(2) = sqrt(power complex / 2) to calculate subband level
# si_std = FS / 4
ni_sigma_fs4 = FS / 4
ni_sigma_fs4_bits = np.log2(ni_sigma_fs4)
nsub_sigma_fs4 = ni_sigma_fs4 * G_subband_sigma
nsub_sigma_fs4_bits = np.log2(nsub_sigma_fs4)
nsub_sigma_re_fs4 = nsub_sigma_fs4 / np.sqrt(N_complex)
nsub_sigma_re_fs4_bits = np.log2(nsub_sigma_re_fs4)
nSST_fs4 = nsub_sigma_fs4**2 * N_int_sub
nSST_fs4_dB = 10 * np.log10(nSST_fs4)
nSST_fs4_bits = np.log2(nSST_fs4)
# si_std = 16
ni_sigma_s16q = sigma_16q
ni_sigma_s16q_bits = np.log2(ni_sigma_s16q) # = 16
nsub_sigma_s16q = ni_sigma_s16q * G_subband_sigma
nsub_sigma_s16q_bits = np.log2(nsub_sigma_s16q)
nsub_sigma_re_s16q = nsub_sigma_s16q / np.sqrt(N_complex)
nsub_sigma_re_s16q_bits = np.log2(nsub_sigma_re_s16q)
nSST_s16q = nsub_sigma_s16q**2 * N_int_sub
nSST_s16q_dB = 10 * np.log10(nSST_s16q)
nSST_s16q_bits = np.log2(nSST_s16q)
print("Incoherent white noise input level --> Expected subband level and SST level:")
print()
print(" si_sigma sub_sigma sub_sigma_re = SST")
print(" sub_sigma_im")
print(" value #bits value #bits value #bits value dB #bits")
print(f"{ni_sigma_fs4:9.1f} {ni_sigma_fs4_bits:6.1f} " \
f"{nsub_sigma_fs4:10.1f} {nsub_sigma_fs4_bits:6.1f} " \
f"{nsub_sigma_re_fs4:13.1f} {nsub_sigma_re_fs4_bits:6.1f} " \
f"{nSST_fs4:15e} {nSST_fs4_dB:6.1f} {nSST_fs4_bits:6.1f}")
print(f"{ni_sigma_s16q:9.1f} {ni_sigma_s16q_bits:6.1f} " \
f"{nsub_sigma_s16q:10.1f} {nsub_sigma_s16q_bits:6.1f} " \
f"{nsub_sigma_re_s16q:13.1f} {nsub_sigma_re_s16q_bits:6.1f} " \
f"{nSST_s16q:15e} {nSST_s16q_dB:6.1f} {nSST_s16q_bits:6.1f}, at {dBFS_16q:.1f} dBFS")
```
%% Output
Incoherent white noise input level --> Expected subband level and SST level:
si_sigma sub_sigma sub_sigma_re = SST
sub_sigma_im
value #bits value #bits value #bits value dB #bits
2048.0 11.0 1024.0 10.0 724.1 9.5 2.048000e+11 113.1 37.6
16.0 4.0 8.0 3.0 5.7 2.5 1.250000e+07 71.0 23.6, at -51.2 dBFS
%% Cell type:markdown id:d6c867ae tags:
Conclusion:
* For FS sine input the subband amplitude is 16 bits, so including the sign bit this fits in W_subband = 18b. The one spare bit is to fit special test signals (e.g. first harmonic of FS square wave input)
* For XST the W_crosslet = 16b subband samples can only fit sine signal input <= 0.5 FS
* For sigma = FS / 4 white noise input the subband sigma uses 10 bits, so 9.5 bits for the subband real and imaginary parts. The 4 sigma just fits in FS and corresponds to 2 bits, so including the sign bit the 4 sigma range of the subband real and imag fits in 1 + 9.5 + 2 = 12.5 bits.
%% Cell type:code id:33f37393 tags:
``` python
# From measured SST to SSTq in q^2 units and subband and signal input
# . depends on whether the input was a coherent sine like signal or a white noise like signal
# Fsub gain implementation factors
print(f"G_subband_ampl = {G_subband_ampl} = {np.log2(G_subband_ampl)} bits")
print(f"G_subband_sigma = {G_subband_sigma} = {np.log2(G_subband_sigma)} bits")
print()
# Coherent (WG sine) signal: from SST to subband amplitude and signal input amplitude in q units
sub_SST = SST_fs4 # SST in WG sine frequency subband for si_ampl = FS / 4 = 2048
si_ampl_exp = si_ampl_fs4
sub_power = sub_SST / N_int_sub
sub_ampl = np.sqrt(sub_power) # ampl real = ampl imag = std complex = sqrt(power complex)
si_ampl = sub_ampl / G_subband_ampl
print(f"sub_SST = {sub_SST:e}")
print(f". sub_power = {sub_power}")
print(f". sub_ampl = {sub_ampl}")
print(f". si_ampl = {si_ampl} (si_ampl_exp = {si_ampl_exp})")
print()
# Incoherent white noise signal: from SST to subband sigma and signal input sigma in q units:
sub_SST = nSST_fs4 # SST in any subband for si_sigma = FS / 4 = 2048
si_sigma_exp = ni_sigma_fs4
sub_power = sub_SST / N_int_sub
sub_sigma = np.sqrt(sub_power) # std complex = sqrt(power)
sub_sigma_re = sub_sigma / np.sqrt(N_complex)
sub_sigma_im = sub_sigma / np.sqrt(N_complex)
si_sigma = sub_sigma / G_subband_sigma
print(f"sub_SST = {sub_SST:e}")
print(f". sub_power = {sub_power}")
print(f". sub_sigma = {sub_sigma}")
print(f". sub_sigma_re = {sub_sigma_re}")
print(f". sub_sigma_im = {sub_sigma_im}")
print(f". si_sigma = {si_sigma} (si_sigma_exp = {si_sigma_exp})")
```
%% Output
G_subband_ampl = 8.0 = 3.0 bits
G_subband_sigma = 0.5 = -1.0 bits
sub_SST = 5.242880e+13
. sub_power = 268435456.0
. sub_ampl = 16384.0
. si_ampl = 2048.0 (si_ampl_exp = 2048.0)
sub_SST = 2.048000e+11
. sub_power = 1048576.0
. sub_sigma = 1024.0
. sub_sigma_re = 724.0773439350246
. sub_sigma_im = 724.0773439350246
. si_sigma = 2048.0 (si_sigma_exp = 2048.0)
%% Cell type:markdown id:66d49365 tags:
## 3.4 Crosslet statistics (XST)
The crosslet statistics have W_crosslet = 16b, but use the same LSbit level as the subbands. The subbands have W_subband = 18b and the maximum subband sine amplitude is 16b. Therefore the maximum sine input for no XST overflow is A = 0.5. If subband_weight = 1.0 then the auto correlations of the XST are eqaul to the SST.
%% Cell type:markdown id:ba543d00 tags:
## 3.5 Beamlet statistics (BST)
%% 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
# Beamlet_sum level and BST level for coherent (WG sine) input (similar as for SST)
beamlet_ampl_fs = si_ampl_fs * G_beamlet_sum_ampl # beamlet amplitude for FS signal input sine
beamlet_ampl_fs_bits = np.log2(beamlet_ampl_fs)
BST_fs = beamlet_ampl_fs**2 * N_int_sub
BST_fs_dB = 10 * np.log10(BST_fs)
BST_fs_bits = np.log2(BST_fs)
beamlet_ampl_fs4 = si_ampl_fs4 * G_beamlet_sum_ampl # beamlet amplitude for FS signal input sine
beamlet_ampl_fs4_bits = np.log2(beamlet_ampl_fs4)
BST_fs4 = beamlet_ampl_fs4**2 * N_int_sub
BST_fs4_dB = 10 * np.log10(BST_fs4)
BST_fs4_bits = np.log2(BST_fs4)
beamlet_ampl_50dBFS = si_ampl_50dBFS * G_beamlet_sum_ampl # beamlet amplitude -50dBFS signal input sine
beamlet_ampl_50dBFS_bits = np.log2(beamlet_ampl_50dBFS)
BST_50dBFS = beamlet_ampl_50dBFS**2 * N_int_sub
BST_50dBFS_dB = 10 * np.log10(BST_50dBFS)
BST_50dBFS_bits = np.log2(BST_50dBFS)
beamlet_ampl_s16q = si_ampl_s16q * G_beamlet_sum_ampl # beamlet amplitude for signal input sine with sigma = 16 q
beamlet_ampl_s16q_bits = np.log2(beamlet_ampl_s16q)
BST_s16q = beamlet_ampl_s16q**2 * N_int_sub
BST_s16q_dB = 10 * np.log10(BST_s16q)
BST_s16q_bits = np.log2(BST_s16q)
print("Coherent (WG sine) signal input level --> Expected beamlet level and BST level:")
print()
print("N_ant si_ampl si_sigma beamlet_sigma = BST")
print(" beamlet_ampl")
print(" value #bits value #bits value #bits value dB #bits")
for ni, na in enumerate(N_ant_arr):
print(f"{na:5d} " \
f"{si_ampl_fs:8.1f} {si_ampl_fs_bits:6.1f} " \
f"{si_sigma_fs:9.1f} {si_sigma_fs_bits:6.1f} " \
f"{beamlet_ampl_fs[ni]:12.1f} {beamlet_ampl_fs_bits[ni]:6.1f} " \
f"{BST_fs[ni]:15e} {BST_fs_dB[ni]:6.1f} {BST_fs_bits[ni]:6.1f}, " \
f"at 0 dBFS (= FS sine)")
print()
for ni, na in enumerate(N_ant_arr):
print(f"{na:5d} " \
f"{si_ampl_fs4:8.1f} {si_ampl_fs4_bits:6.1f} " \
f"{si_sigma_fs4:9.1f} {si_sigma_fs4_bits:6.1f} " \
f"{beamlet_ampl_fs4[ni]:12.1f} {beamlet_ampl_fs4_bits[ni]:6.1f} " \
f"{BST_fs4[ni]:15e} {BST_fs4_dB[ni]:6.1f} {BST_fs4_bits[ni]:6.1f}, " \
f"at {20*np.log10(1 / 4**2):.1f} dBFS (= FS / 4)")
print()
for ni, na in enumerate(N_ant_arr):
print(f"{na:5d} " \
f"{si_ampl_50dBFS:8.1f} {si_ampl_50dBFS_bits:6.1f} " \
f"{si_sigma_50dBFS:9.1f} {si_sigma_50dBFS_bits:6.1f} " \
f"{beamlet_ampl_50dBFS[ni]:12.1f} {beamlet_ampl_50dBFS_bits[ni]:6.1f} " \
f"{BST_50dBFS[ni]:15e} {BST_50dBFS_dB[ni]:6.1f} {BST_50dBFS_bits[ni]:6.1f}, "
f"at -50 dBFS (= FS / {10**(50/20):.0f})")
print()
for ni, na in enumerate(N_ant_arr):
print(f"{na:5d} " \
f"{si_ampl_s16q:8.1f} {si_ampl_s16q_bits:6.1f} " \
f"{si_sigma_s16q:9.1f} {si_sigma_s16q_bits:6.1f} " \
f"{beamlet_ampl_s16q[ni]:12.1f} {beamlet_ampl_s16q_bits[ni]:6.1f} " \
f"{BST_s16q[ni]:15e} {BST_s16q_dB[ni]:6.1f} {BST_s16q_bits[ni]:6.1f}, " \
f"at {dBFS_16q:.1f} dBFS (= FS / {10**(-dBFS_16q/20):.0f})")
```
%% Output
Coherent (WG sine) signal input level --> Expected beamlet level and BST level:
N_ant si_ampl si_sigma beamlet_sigma = BST
beamlet_ampl
value #bits value #bits value #bits value dB #bits
1 8192.0 13.0 5792.6 4.0 65536.0 16.0 8.388608e+14 149.2 49.6, at 0 dBFS (= FS sine)
12 8192.0 13.0 5792.6 4.0 786432.0 19.6 1.207960e+17 170.8 56.7, at 0 dBFS (= FS sine)
24 8192.0 13.0 5792.6 4.0 1572864.0 20.6 4.831838e+17 176.8 58.7, at 0 dBFS (= FS sine)
48 8192.0 13.0 5792.6 4.0 3145728.0 21.6 1.932735e+18 182.9 60.7, at 0 dBFS (= FS sine)
96 8192.0 13.0 5792.6 4.0 6291456.0 22.6 7.730941e+18 188.9 62.7, at 0 dBFS (= FS sine)
1 2048.0 11.0 1448.2 4.0 16384.0 14.0 5.242880e+13 137.2 45.6, at -24.1 dBFS (= FS / 4)
12 2048.0 11.0 1448.2 4.0 196608.0 17.6 7.549747e+15 158.8 52.7, at -24.1 dBFS (= FS / 4)
24 2048.0 11.0 1448.2 4.0 393216.0 18.6 3.019899e+16 164.8 54.7, at -24.1 dBFS (= FS / 4)
48 2048.0 11.0 1448.2 4.0 786432.0 19.6 1.207960e+17 170.8 56.7, at -24.1 dBFS (= FS / 4)
96 2048.0 11.0 1448.2 4.0 1572864.0 20.6 4.831838e+17 176.8 58.7, at -24.1 dBFS (= FS / 4)
1 25.9 4.7 18.3 4.2 207.2 7.7 8.388608e+09 99.2 33.0, at -50 dBFS (= FS / 316)
12 25.9 4.7 18.3 4.2 2486.9 11.3 1.207960e+12 120.8 40.1, at -50 dBFS (= FS / 316)
24 25.9 4.7 18.3 4.2 4973.8 12.3 4.831838e+12 126.8 42.1, at -50 dBFS (= FS / 316)
48 25.9 4.7 18.3 4.2 9947.7 13.3 1.932735e+13 132.9 44.1, at -50 dBFS (= FS / 316)
96 25.9 4.7 18.3 4.2 19895.3 14.3 7.730941e+13 138.9 46.1, at -50 dBFS (= FS / 316)
1 22.6 4.5 16.0 4.0 181.0 7.5 6.400000e+09 98.1 32.6, at -51.2 dBFS (= FS / 362)
12 22.6 4.5 16.0 4.0 2172.2 11.1 9.216000e+11 119.6 39.7, at -51.2 dBFS (= FS / 362)
24 22.6 4.5 16.0 4.0 4344.5 12.1 3.686400e+12 125.7 41.7, at -51.2 dBFS (= FS / 362)
48 22.6 4.5 16.0 4.0 8688.9 13.1 1.474560e+13 131.7 43.7, at -51.2 dBFS (= FS / 362)
96 22.6 4.5 16.0 4.0 17377.9 14.1 5.898240e+13 137.7 45.7, at -51.2 dBFS (= FS / 362)
%% Cell type:code id:06c7b393 tags:
``` python
# Beamlet level and BST level for incoherent white noise input (similar as for SST)
# si_std = FS / 4
nbeamlet_sigma_fs4 = ni_sigma_fs4 * G_beamlet_sum_sigma
nbeamlet_sigma_fs4_bits = np.log2(nbeamlet_sigma_fs4)
nbeamlet_sigma_re_fs4 = nbeamlet_sigma_fs4 / np.sqrt(N_complex)
nbeamlet_sigma_re_fs4_bits = np.log2(nbeamlet_sigma_re_fs4)
nBST_fs4 = nbeamlet_sigma_fs4**2 * N_int_sub
nBST_fs4_dB = 10 * np.log10(nBST_fs4)
nBST_fs4_bits = np.log2(nBST_fs4)
# si_std = 16
nbeamlet_sigma_s16q = ni_sigma_s16q * G_beamlet_sum_sigma
nbeamlet_sigma_s16q_bits = np.log2(nbeamlet_sigma_s16q)
nbeamlet_sigma_re_s16q = nbeamlet_sigma_s16q / np.sqrt(N_complex)
nbeamlet_sigma_re_s16q_bits = np.log2(nbeamlet_sigma_re_s16q)
nBST_s16q = nbeamlet_sigma_s16q**2 * N_int_sub
nBST_s16q_dB = 10 * np.log10(nBST_s16q)
nBST_s16q_bits = np.log2(nBST_s16q)
print("Incoherent white noise input level --> Expected beamlet level and BST level:")
print()
print("N_ant si_sigma beamlet_sigma beamlet_sigma_re = BST")
print(" beamlet_sigma_im")
print(" value #bits value #bits value #bits value dB #bits")
for ni, na in enumerate(N_ant_arr):
print(f"{na:5d} " \
f"{ni_sigma_fs4:9.1f} {ni_sigma_fs4_bits:6.1f} " \
f"{nbeamlet_sigma_fs4[ni]:10.1f} {nbeamlet_sigma_fs4_bits[ni]:6.1f} " \
f"{nbeamlet_sigma_re_fs4[ni]:16.1f} {nbeamlet_sigma_re_fs4_bits[ni]:6.1f} " \
f"{nBST_fs4[ni]:15e} {nBST_fs4_dB[ni]:6.1f} {nBST_fs4_bits[ni]:6.1f}")
print()
for ni, na in enumerate(N_ant_arr):
print(f"{na:5d} " \
f"{ni_sigma_s16q:9.1f} {ni_sigma_s16q_bits:6.1f} " \
f"{nbeamlet_sigma_s16q[ni]:10.1f} {nbeamlet_sigma_s16q_bits[ni]:6.1f} " \
f"{nbeamlet_sigma_re_s16q[ni]:16.1f} {nbeamlet_sigma_re_s16q_bits[ni]:6.1f} " \
f"{nBST_s16q[ni]:15e} {nBST_s16q_dB[ni]:6.1f} {nBST_s16q_bits[ni]:6.1f}, at {dBFS_16q:.1f} dBFS")
```
%% Output
Incoherent white noise input level --> Expected beamlet level and BST level:
N_ant si_sigma beamlet_sigma beamlet_sigma_re = BST
beamlet_sigma_im
value #bits value #bits value #bits value dB #bits
1 2048.0 11.0 1024.0 10.0 724.1 9.5 2.048000e+11 113.1 37.6
12 2048.0 11.0 3547.2 11.8 2508.3 11.3 2.457600e+12 123.9 41.2
24 2048.0 11.0 5016.6 12.3 3547.2 11.8 4.915200e+12 126.9 42.2
48 2048.0 11.0 7094.5 12.8 5016.6 12.3 9.830400e+12 129.9 43.2
96 2048.0 11.0 10033.1 13.3 7094.5 12.8 1.966080e+13 132.9 44.2
1 16.0 4.0 8.0 3.0 5.7 2.5 1.250000e+07 71.0 23.6, at -51.2 dBFS
12 16.0 4.0 27.7 4.8 19.6 4.3 1.500000e+08 81.8 27.2, at -51.2 dBFS
24 16.0 4.0 39.2 5.3 27.7 4.8 3.000000e+08 84.8 28.2, at -51.2 dBFS
48 16.0 4.0 55.4 5.8 39.2 5.3 6.000000e+08 87.8 29.2, at -51.2 dBFS
96 16.0 4.0 78.4 6.3 55.4 5.8 1.200000e+09 90.8 30.2, at -51.2 dBFS
%% Cell type:markdown id:8829fd41 tags:
If subband_weight = 1.0 and beamlet_weight = 1.0 and N_ant = 1 then the BST are eqaul to the SST
%% Cell type:markdown id:cdb20624 tags:
## 3.6 Beamlet output
The beamlet output is W_beamlet = 8 bit. The beamlet has a sign bit about 1 bit for the sigma and about 2 bits to fit a range of 4 sigma, so about 4 bits can carry the noise signal. The extra 4 bits are for some RFI and to fit differences in subband noise level due to the antenna and RCU2 band filter shape. The subband noise level can be equalized using the subband weights, to have more dynamic range for RFI the beamlet.
Choosing FPGA_beamlet_output_scale_RW = 1 sqrt(N_ant) makes that the beamlet output level for noise input is equal to that of N_ant = 1 for all N_ant. The BST can be used to check whether the beamlet output will fit.
%% Cell type:markdown id:d2086ec5 tags:
# Appendix 1: DFT of real input DC, sine and noise
# Appendix 1: DFT of real input DC, sine, block and noise
%% Cell type:code id:def6eba7 tags:
``` python
# Show DFT bin levels for real input
# . Sine with A = 1, so power A^2 / 2, yields phasor in both side bands, each with amplitude
# A / N_sidebands = 0.5, because both with equal power, such that P_sine = A^2/2. The factor
# is 1 / N and not 1 / sqrt(N)is because the input sine is coherent.
# . For DC at bin frequency 0 the real input gain is 1.0, because DC has only one phasor
# component of frequency 0.
# . For white noise the power in each bin is a factor N_fft less, so the std() is
# sqrt(N_fft) less, because the white noise is incoherent.
#
# . Input level
ampl = 1.0
sigma = ampl / np.sqrt(2)
sigma = 4
# . DFT size
N_bins = N_fft // 2 + 1 # positive frequency bins including DC and f_s/2
# . select a bin
i_bin = 200 # bin index in range(N_bins)
dc_bin = 0 # DC
# . 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_fft * T_s # DFT period
t_axis = np.linspace(0, T_fft, N_fft, endpoint=False)
f_axis = np.linspace(0, f_s, N_fft, endpoint=False)
f_axis_fft = f_axis - f_s/2 # fftshift axis
f_axis_rfft = f_axis[0:N_bins] # positive frequency bins
bw_bin = f_s / N_fft # bin band width
f_bin = i_bin * bw_bin # bin frequency
# . create sine at bin + DC, use cos to see DC at i_bin = 0
s = ampl * np.cos(2 * np.pi * f_bin * t_axis)
dc = ampl * np.cos(2 * np.pi * dc_bin * t_axis) # equivalent to dc = ampl
noise = np.random.randn(N_fft)
noise *= sigma / np.std(noise) # apply requested sigma
b = ampl * np.sign(s) # block wave, sign: -1 if x < 0, 0 if x==0, 1 if x > 0
x = s + dc
y = noise
# . DFT using complex input fft()
S_fft = np.fft.fftshift(np.fft.fft(s) / N_fft)
B_fft = np.fft.fftshift(np.fft.fft(b) / N_fft)
X_fft = np.fft.fftshift(np.fft.fft(x) / N_fft)
Y_fft = np.fft.fftshift(np.fft.fft(y) / N_fft)
# . DFT using real input rfft()
S_rfft = np.fft.rfft(s) / N_fft
B_rfft = np.fft.rfft(b) / N_fft
X_rfft = np.fft.rfft(x) / N_fft
Y_rfft = np.fft.rfft(y) / N_fft
# Plot sine spectrum
# . DSB = double sideband
# . SSB = single sideband (= DC + positive frequencies)
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.title('DFT of real input sine using fft (DSB)')
plt.plot(f_axis_fft, abs(X_fft))
plt.grid()
plt.subplot(1, 2, 2)
plt.title('DFT of real input sine using rfft (SSB)')
plt.plot(f_axis_rfft, abs(X_rfft))
plt.grid()
# Plot block spectrum
# . DSB = double sideband
# . SSB = single sideband (= DC + positive frequencies)
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.title('DFT of real input block using fft (DSB)')
plt.plot(f_axis_fft, abs(B_fft))
plt.grid()
plt.subplot(1, 2, 2)
plt.title('DFT of real input block using rfft (SSB)')
plt.plot(f_axis_rfft, abs(B_rfft))
plt.grid()
# Plot noise spectrum
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.title('DFT of real input noise using fft (DSB)')
plt.plot(f_axis_fft, abs(Y_fft))
plt.grid()
plt.subplot(1, 2, 2)
plt.title('DFT of real input noise using rfft (SSB)')
plt.plot(f_axis_rfft, abs(Y_rfft))
plt.grid()
print("The DFT of the sine plot shows:")
print(f". G_fft_real_input_dc = {G_fft_real_input_dc}")
print(f". G_fft_real_input_sine = {G_fft_real_input_sine}")
print()
S_max = max(abs(S_rfft))
B_max = max(abs(B_rfft))
print(f"The DFT of the block plot shows that the first harmonic has an amplitude of 4/pi * A/2 = {B_max}, \
which is larger than A / 2 = {S_max} for sine input. Hence the bin samples need 1 bit more than for a full \
scale sine, because to also fit e.g. this harmonic of a block wave.")
```
%% Output
The DFT of the sine plot shows:
. G_fft_real_input_dc = 1
. G_fft_real_input_sine = 0.5
The DFT of the block plot shows that the first harmonic has an amplitude of 4/pi * A/2 = 0.6364949321522198, which is larger than A / 2 = 0.5000000000000007 for sine input. Hence the bin samples need 1 bit more than for a full scale sine, because to also fit e.g. this harmonic of a block wave.
%% Cell type:code id:2e386180 tags:
``` python
# Log amplitude, sigma and power of the sine bins
# . For sine input the bin is a constant phasor when the sine frequency corresponds to
# the center of the bin and a rotating phasor when the sine frequency is off center,
# the phasor then rotates with the delta frequency = bin frequency - sine frequency.
# . The input sine has ampl = A and power A**2 / 2.
# . The real input results in N_sidebands = 2 bins.
# . The amplitude of the phasor is equal to the std() of a rotating phasor.
# . The amplitude of the bin phasor is A/2, so the bin phasor A/2 exp(jwt) has power
# (A/2)**2
# . The total power in the N_sidebands = 2 bins is eual to the input power as expected
# . input sine
sin_std = np.std(s)
sine_power = np.sum(s**2) / N_fft
print(f"sine ampl = {ampl:.4f}")
print(f"sine sigma = {sin_std:.4f} (= {sigma:.4f})")
print(f"sine power = {sin_std**2:.4f} (= {sine_power:.4f})")
print()
# . fft bin
bin_ampl = np.max(np.abs(S_rfft))
bin_re = np.max(S_rfft.real)
bin_im = np.max(S_rfft.imag)
bin_power = bin_ampl**2
bins_power = bin_power * N_sidebands
# . Model bin_std using a rotating phasor with bw_bin frequency to have one complete
# period in t_axis
bin_phasor = bin_ampl * np.exp(2 * np.pi * 1j * bw_bin * t_axis)
bin_std = np.std(bin_phasor)
plt.figure(figsize=(10, 4))
plt.title('Bin phasor to model bin_std')
plt.plot(t_axis, bin_phasor.real, 'r', t_axis, bin_phasor.imag, 'b')
plt.grid()
print(f"sine bin ampl = {bin_ampl:.4f}")
print(f"sine bin re = {bin_re:.4f}")
print(f"sine bin im = {bin_im:.4f}")
print(f"sine bin std = {bin_std:.4f}")
print(f"sine bin power = {bin_power:.4f}")
print(f"sine all bins power = {bins_power:.4f}")
```
%% Output
sine ampl = 1.0000
sine sigma = 0.7071 (= 0.7071)
sine sigma = 0.7071 (= 4.0000)
sine power = 0.5000 (= 0.5000)
sine bin ampl = 0.5000
sine bin re = 0.5000
sine bin im = 0.0000
sine bin std = 0.5000
sine bin power = 0.2500
sine all bins power = 0.5000
%% Cell type:code id:97e9a32d tags:
``` python
# Log sigma and power of the noise bins (incoherent signal input)
# . input noise
noise_std = np.std(y)
noise_power = np.sum(y**2) / N_fft
print(f"noise sigma = {noise_std:.4f} (= {sigma:.4f})")
print(f"noise power = {noise_std**2:.4f} (= {noise_power:.4f})")
print()
# . fft bin
# . The white noise will appear equally in all bins, therefore the bin_std can
# be modelled by averaging over all bins. This however does cause small
# differences in fft input and output power, due to the DC bin (?)
bin_std = np.std(Y_rfft)
bin_power = bin_std**2
bin_re_std = np.std(Y_rfft.real)
bin_re_power = bin_re_std**2
bin_im_std = np.std(Y_rfft.imag)
bin_im_power = bin_im_std**2
bins_power = bin_power * N_fft
print(f"N_fft = {N_fft}")
print(f"sqrt(N_fft) = {np.sqrt(N_fft)}")
print(f"sigma / std(Y_fft) = {sigma / np.std(Y_fft):f}")
print(f"sigma / std(Y_rfft) = {sigma / bin_std:f}")
print()
print(f"noise bin std (fft) = {np.std(Y_fft):f}")
print(f"noise bin std (rfft) = {bin_std:f}")
print(f"noise bin.re std = {bin_re_std:f}")
print(f"noise bin.im std = {bin_im_std:f}")
print(f"noise bin power = {bin_power:f}")
print(f"noise bin.re power + bin.im power = {bin_re_power + bin_im_power:f}")
print(f"noise bins power = {bins_power:f} (= {noise_power:f})")
print()
print("The ratio of real input noise std and DFT bin noise std shows:")
print(f". G_fft_real_input_noise = {G_fft_real_input_noise} = (1 / sqrt({N_fft}))")
```
%% Output
noise sigma = 0.7071 (= 0.7071)
noise power = 0.5000 (= 0.5001)
noise sigma = 4.0000 (= 4.0000)
noise power = 16.0000 (= 16.0014)
N_fft = 1024
sqrt(N_fft) = 32.0
sigma / std(Y_fft) = 31.996921
sigma / std(Y_rfft) = 32.018608
sigma / std(Y_fft) = 32.047266
sigma / std(Y_rfft) = 32.082825
noise bin std (fft) = 0.022099
noise bin std (rfft) = 0.022084
noise bin.re std = 0.015340
noise bin.im std = 0.015887
noise bin power = 0.000488
noise bin.re power + bin.im power = 0.000488
noise bins power = 0.499419 (= 0.500120)
noise bin std (fft) = 0.124816
noise bin std (rfft) = 0.124677
noise bin.re std = 0.090000
noise bin.im std = 0.086281
noise bin power = 0.015544
noise bin.re power + bin.im power = 0.015544
noise bins power = 15.917496 (= 16.001391)
The ratio of real input noise std and DFT bin noise std shows:
. G_fft_real_input_noise = 0.03125 = (1 / sqrt(1024))
%% Cell type:markdown id:6d9f356a tags:
Conclusion:
* For coherent sine input is easiest to calculate power via the amplitude of the single bin phasor
* For incoherent white noise input is easiest to calculate power via the std of all bins
%% Cell type:code id:17aee663 tags:
``` python
# DFT of square wave
```
......
doc/erko_howto_tools.txt
+ 7
− 2
View file @ 2aefd0ba
......@@ -417,7 +417,9 @@ then:
> ping 10.87.0.186 (= dop386)
> ssh -X dop386.astron.nl
> ssh -X kooistra@10.87.0.186 # dop386 from ASTROM
> ssh -X kooistra@10.87.0.186 # dop386 from ASTRON
> sshdop386 # dop386 terminal, alias in .bashrc
> dop386t # dop386 terminal, alias in .bashrc
> ssh -J bastion.astron.nl kooistra@10.87.0.186 # dop386 from home
1) gitlab
......@@ -481,12 +483,15 @@ then:
64 sdp_rw.py --host 10.99.0.250 --port 4842 -r xst_integration_interval
71 sdp_rw.py --host 10.99.0.250 --port 4842 -r signal_input_bsn
4) tcpdump
> ifconfig # to find ethernet port
> sudo tcpdump -vvXXSnelfi enp5s0 port 5001 # for UDP only
> sudo tcpdump -vvXXSnelfi enp5s0 port 5001 > tcpdump.txt (> is new file, >> is append)
5) opcua-client & # to verify sdp_rw.py
- opc.tcp://10.99.0.250:4842 RPi4 connect
- write via RW point attribute window
- right mouse subscribe/unsubsribe to data change
*******************************************************************************
* Polarion:
......
doc/sdp_useful_commands_erko.txt 0 → 100644
+ 503
− 0
View file @ 2aefd0ba
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment