iq imbalance compensation

bug_fixes_integration_tx
David Michaeli 2021-12-06 18:26:47 +02:00
rodzic d465c5b425
commit 120ea6b79b
3 zmienionych plików z 241 dodań i 72 usunięć

Wyświetl plik

@ -1,12 +1,15 @@
from PySimpleGUI.PySimpleGUI import Canvas, Column
import time
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg, NavigationToolbar2Tk
import SoapySDR
from SoapySDR import SOAPY_SDR_RX, SOAPY_SDR_TX, SOAPY_SDR_CS16
from PySimpleGUI import Window, WIN_CLOSED, Slider, Button, theme, Text, Radio, Image
from PySimpleGUI import Window, WIN_CLOSED, Slider, Button, theme, Text, Radio, Image, InputText, Canvas
from numpy.lib.arraypad import pad
##
## WINDOW FUNCTIONS
def draw_figure_with_toolbar(canvas, fig, canvas_toolbar):
if canvas.children:
for child in canvas.winfo_children():
@ -33,124 +36,246 @@ def create_window():
layout = [
[Column(
layout=[
[Text('Receiver LO')],
[Text('Transmitter Fc')]
[Text('Receiver LO [Hz]:'), InputText('868000000', key='RxFreq'), Text('RX Bandwidth: 2.5 MHz')],
[Text('Transmitter Fc [Hz]:'), Text(key='TxCWFreq')],
]
)]
)],
[Canvas(key='controls_cv')],
[Column(
layout=[
[Canvas(key='fig_cv', size=(400*2, 400))]
[Canvas(key='fig_cv', size=(300*2, 300))]
],
background_color='#DAE0E6',
pad=(0,0)
)],
[Button("Correct", size=(10,1)), Button("Exit", size=(10,1))],
[Text('RSSI:', key='rssi'), Text('Alpha:', key='alpha'), Text('Cos(Phi):', key='cos_phi')],
[Canvas(key='controls_params_cv')],
[Column(
layout=[
[Canvas(key='params_cv', size=(200*3, 200))]
],
background_color='#DAE0E6',
pad=(0,0)
)],
[Button("Correct", size=(10,1)), Button("Exit", size=(10,1)), Button("Run", size=(10,1))],
]
window = Window("CaribouLite PlayGround", layout, location=(800,400))
return window
def update_iq_graphs(window, I,Q,I_new, Q_new):
plt.figure(1).clf()
fig = plt.gcf()
DPI = fig.get_dpi()
fig.set_size_inches(304*2/float(DPI), 304/float(DPI))
plt.subplot(121)
plt.scatter(I, Q)
plt.title('I/Q original')
plt.subplot(122)
plt.scatter(I_new, Q_new)
plt.title('I/Q corrected')
plt.grid()
draw_figure_with_toolbar(window['fig_cv'].TKCanvas, fig, window['controls_cv'].TKCanvas)
def estimate_correct(rx_buff):
s_real = rx_buff[::2].astype(np.float32)
s_imag = rx_buff[1::2].astype(np.float32)
def update_est_graphs(window, freq_diff, alpha_vec, cos_phi_vec, rssi_vec):
plt.figure(2).clf()
fig = plt.gcf()
DPI = fig.get_dpi()
fig.set_size_inches(204*3/float(DPI), 204/float(DPI))
plt.subplot(131)
plt.plot(freq_diff, alpha_vec)
plt.title('Alpha')
plt.subplot(132)
plt.plot(freq_diff, cos_phi_vec)
plt.title('cos(phi)')
plt.subplot(133)
plt.plot(freq_diff, rssi_vec)
plt.title('rssi')
plt.grid()
draw_figure_with_toolbar(window['params_cv'].TKCanvas, fig, window['controls_params_cv'].TKCanvas)
c = s_real + 1j * s_imag
c_no_bias = c - c.mean()
##
## I/Q Correction Functions
def fix_iq_blind(x):
z = x - np.mean(x)
p_in = np.var(z)
# scale Q to have unit amplitude
Q_amp = np.sqrt(2*np.mean(s_imag**2))
c /= Q_amp
I, Q = c.real, c.imag
theta1 = (-1) * np.mean(np.sign(z.real) * z.imag)
theta2 = np.mean(np.abs(z.real))
theta3 = np.mean(np.abs(z.imag))
c1 = theta1/theta2
c2 = np.sqrt( ( (theta3**2 - theta1**2) / theta2**2 ) )
g = theta3/theta2
phi = np.arcsin(theta1/theta3)
wI = c2 * z.real
wQ = c1 * z.real + z.imag
z_out = wI + 1j * wQ
return (z_out, c1, c2, p_in)
def fix_iq_imbalance(x):
# remove DC and save input power
z = x - np.mean(x)
p_in = np.var(z)
# scale Q to have unit amplitude (remember we're assuming a single input tone)
Q_amp = np.sqrt(2*np.mean(z.imag**2)) # <== should be z.imag**2????
z /= Q_amp
I, Q = z.real, z.imag
alpha_est = np.sqrt(2*np.mean(I**2))
sin_phi_est = (2/alpha_est)*np.mean(I*Q)
cos_phi_est = np.sqrt(1-sin_phi_est**2)
cos_phi_est = np.sqrt(1 - sin_phi_est**2)
print('phase error: ', np.arccos(cos_phi_est)*360/3.1415, ' degrees')
print('amplitude error: ', 20*np.log10(alpha_est), ' dB' )
I_new_p = (1/alpha_est)*I
Q_new_p = (-sin_phi_est/alpha_est)*I + Q
I_new_p = (1/alpha_est) * I
Q_new_p = (-sin_phi_est/alpha_est) * I + Q
y = (I_new_p + 1j*Q_new_p) / cos_phi_est
y = (I_new_p + 1j*Q_new_p)/cos_phi_est
return (s_real, s_imag, alpha_est, sin_phi_est, cos_phi_est, I_new_p, Q_new_p)
print ('phase error:', np.arccos(cos_phi_est)*360/2/np.pi, 'degrees')
print ('amplitude error:', 20*np.log10(alpha_est), 'dB')
z_out = y*np.sqrt(p_in/np.var(y))
return (z_out, alpha_est, cos_phi_est, p_in)
##
## Soapy Control functions
def setup_receiver(sdr, channel, freq_hz):
use_agc = False # Use or don't use the AGC
# The wide channel parameters
sdr.setGainMode(SOAPY_SDR_RX, channel, use_agc) # Set the gain mode
sdr.setGain(SOAPY_SDR_RX, channel, 70) # Set the gain
sdr.setFrequency(SOAPY_SDR_RX, channel, freq_hz) # Tune the LO
sdr.setBandwidth(SOAPY_SDR_RX, channel, 1.6e6)
rx_stream = sdr.setupStream(SOAPY_SDR_RX, SOAPY_SDR_CS16, [channel]) # Setup data stream
return rx_stream
def update_receiver_freq(sdr, stream, channel, freq_hz):
sdr.setFrequency(SOAPY_SDR_RX, channel, freq_hz)
def setup_transmitter(sdr, channel, freq_hz):
sdr.setFrequency(SOAPY_SDR_TX, channel, 867.5e6)
sdr.setGain(SOAPY_SDR_TX, channel, 60)
tx_stream = sdr.setupStream(SOAPY_SDR_TX, SOAPY_SDR_CS16, [channel], dict(CW="1"))
return tx_stream
def update_transmitter_freq(sdr, stream, channel, freq_hz):
sdr.deactivateStream(stream)
sdr.setFrequency(SOAPY_SDR_TX, channel, freq_hz)
sdr.activateStream(stream)
##
## MAIN
def main():
# Data and Source Configuration
tx_chan = 0 # S1G = 0, 6G = 1
rx_chan = 1 # S1G = 0, 6G = 1
N = 16384 # Nu5mber of complex samples per transfer
fs = 4e6 # Radio sample Rate
freq = 868e6 # LO tuning frequency in Hz
use_agc = True # Use or don't use the AGC
gain = 60
n_reads = 0
tx_chan = 0 # S1G = 0
rx_chan = 1 # 6G = 1
N = 16384 # Number of complex samples per transfer
rx_buff = np.empty(2 * N, np.int16) # Create memory buffer for data stream
timeout_us = int(5e6)
# Initialize CaribouLite
# Initialize CaribouLite Soapy
sdr = SoapySDR.Device(dict(driver="Cariboulite")) # Create Cariboulite
# The wide channel parameters
sdr.setSampleRate(SOAPY_SDR_RX, rx_chan, fs) # Set sample rate
sdr.setGainMode(SOAPY_SDR_RX, rx_chan, use_agc) # Set the gain mode
sdr.setGain(SOAPY_SDR_RX, rx_chan, gain) # Set the gain
sdr.setFrequency(SOAPY_SDR_RX, rx_chan, freq) # Tune the LO
sdr.setBandwidth(SOAPY_SDR_RX, rx_chan, 2.5e6)
# The S1G channel parameter
sdr.setFrequency(SOAPY_SDR_TX, tx_chan, 868.5e6)
sdr.setGain(SOAPY_SDR_TX, tx_chan, 50)
cw_tx_stream = sdr.setupStream(SOAPY_SDR_TX, SOAPY_SDR_CS16, [tx_chan], dict(CW="1"))
sdr.activateStream(cw_tx_stream)
# Create data buffer and start streaming samples to it
rx_stream = sdr.setupStream(SOAPY_SDR_RX, SOAPY_SDR_CS16,
[rx_chan]) # Setup data stream
sdr.activateStream(rx_stream) # this turns the radio on
rx_stream = setup_receiver(sdr, rx_chan, 868e6)
tx_stream = setup_transmitter(sdr, tx_chan, 868.5e6)
sdr.activateStream(tx_stream)
sdr.activateStream(rx_stream)
# create the window
window = create_window()
while True:
event, values = window.read(timeout = 20)
if (event == "Exit" or event == WIN_CLOSED):
break;
elif event == "Correct":
print("clicked correct")
if (event == "Correct"):
n_reads += 1
sr = sdr.readStream(rx_stream, [rx_buff], N, timeoutUs=timeout_us)
freq = float(values['RxFreq'])
tx_freq = freq + 0.5e6
print("setting frequency ", freq)
window['TxCWFreq'].update(tx_freq)
update_transmitter_freq(sdr, tx_stream, tx_chan, tx_freq)
update_receiver_freq(sdr, rx_stream, rx_chan, freq)
time.sleep(0.5)
sr = sdr.readStream(rx_stream, [rx_buff], N, timeoutUs=int(5e6))
# Make sure that the proper number of samples was read
rc = sr.ret
if (rc != N):
print("Error Reading Samples from Device (error code = %d)!" % rc)
break;
I,Q, alpha_est, sin_phi_est, cos_phi_est, I_new, Q_new = estimate_correct(rx_buff)
s_real = rx_buff[::2].astype(np.float32)
s_imag = rx_buff[1::2].astype(np.float32)
z = s_real + 1j*s_imag
#(z_out, alpha_est, cos_phi_est, p_in) = fix_iq_imbalance(z)
(z_out, alpha_est, cos_phi_est, p_in) = fix_iq_blind(z)
rssi = 20*np.log10(p_in)
plt.figure(1).clf()
fig = plt.gcf()
DPI = fig.get_dpi()
fig.set_size_inches(404*2/float(DPI), 404/float(DPI))
plt.subplot(121)
plt.scatter(I, Q)
plt.title('I/Q original')
plt.subplot(122)
plt.scatter(I_new, Q_new)
plt.title('I/Q corrected')
plt.grid()
draw_figure_with_toolbar(window['fig_cv'].TKCanvas, fig, window['controls_cv'].TKCanvas)
window['rssi'].update('RSSI: %f dBm' % rssi)
window['alpha'].update('Alpha: %f' % alpha_est)
window['cos_phi'].update('Cos(Phi): %f' % cos_phi_est)
update_iq_graphs(window, s_real, s_imag, z_out.real, z_out.imag)
elif event == "Run":
print("clicked run")
bw = 2.5e6
rx_freq = float(values['RxFreq'])
update_receiver_freq(sdr, rx_stream, rx_chan, rx_freq)
freq_num = 60
tx_freq_vec = rx_freq + np.linspace(-bw/2, bw/2, freq_num)
alpha_est_vec = np.empty(tx_freq_vec.size, np.float32)
cos_phi_vec = np.empty(tx_freq_vec.size, np.float32)
rssi_vec = np.empty(tx_freq_vec.size, np.float32)
index = 0
for tx_freq in tx_freq_vec:
window['TxCWFreq'].update(tx_freq)
update_transmitter_freq(sdr, tx_stream, tx_chan, tx_freq)
time.sleep(0.1)
sr = sdr.readStream(rx_stream, [rx_buff], N, timeoutUs=int(5e6))
# Make sure that the proper number of samples was read
rc = sr.ret
if (rc != N):
print("Error Reading Samples from Device (error code = %d)!" % rc)
break;
s_real = rx_buff[::2].astype(np.float32)
s_imag = rx_buff[1::2].astype(np.float32)
z = s_real + 1j*s_imag
#(z_out, alpha_est, cos_phi_est, p_in) = fix_iq_imbalance(z)
(z_out, alpha_est, cos_phi_est, p_in) = fix_iq_blind(z)
alpha_est_vec[index] = 20*np.log10(alpha_est)
cos_phi_vec[index] = np.arccos(cos_phi_est)*360/2/np.pi
rssi_vec[index] = 20*np.log10(p_in)
window['rssi'].update('RSSI: %f dBm' % rssi_vec[index])
window['alpha'].update('Alpha: %f' % alpha_est_vec[index])
window['cos_phi'].update('Cos(Phi): %f' % cos_phi_vec[index])
update_iq_graphs(window, s_real, s_imag, z_out.real, z_out.imag)
index += 1
time.sleep(0.1)
update_est_graphs(window, tx_freq_vec-rx_freq, alpha_est_vec, cos_phi_vec, rssi_vec)
# Stop streaming and close the connection to the radio
window.close()
sdr.deactivateStream(rx_stream)
sdr.deactivateStream(cw_tx_stream)
sdr.deactivateStream(tx_stream)
sdr.closeStream(rx_stream)
sdr.closeStream(tx_stream)
@ -158,7 +283,9 @@ def main():
## Global Area
# Enumerate the Soapy devices
#res = SoapySDR.Device.enumerate()
#for result in res:
# print(result)
res = SoapySDR.Device.enumerate()
for result in res:
print(result)
# run the program
main()

Wyświetl plik

@ -240,6 +240,48 @@ int SampleQueue::ReadSamples(sample_complex_float* buffer, size_t num_elements,
buffer[i].q = (float)interm_native_buffer[i].q / max_val;
}
double sumI = 0.0;
double sumQ = 0.0;
for (int i = 0; i < res; i++)
{
sumI += buffer[i].i;
sumQ += buffer[i].q;
}
sumI /= (double)res;
sumQ /= (double)res;
for (int i = 0; i < res; i++)
{
buffer[i].i -= sumI;
buffer[i].q -= sumQ;
}
double theta1 = 0.0;
double theta2 = 0.0;
double theta3 = 0.0;
for (int i = 0; i < res; i++)
{
int sign_I = (buffer[i].i > 0 ? 1 : -1);
int sign_Q = (buffer[i].q > 0 ? 1 : -1);
theta1 += sign_I * buffer[i].q;
theta2 += sign_I * buffer[i].i;
theta3 += sign_Q * buffer[i].q;
}
theta1 = - theta1 / (double)res;
theta2 /= (double)res;
theta3 /= (double)res;
double c1 = theta1 / theta2;
double c2 = sqrt( (theta3*theta3 - theta1*theta1) / (theta2*theta2) );
for (int i = 0; i < res; i++)
{
float ii = buffer[i].i;
float qq = buffer[i].q;
buffer[i].i = c2 * ii;
buffer[i].q = c1 * ii + qq;
}
return res;
}