``````

In [1]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal
import scipy.stats
import scipy.special

``````
``````

In [2]:

x = np.fromfile('ed4yae_10GHz.c64', dtype = 'complex64')

``````
``````

In [3]:

avg = 100
x_pwr_avg = np.sum(np.abs(x[:x.size//100*100].reshape((-1,100)))**2, axis = 1)
plt.plot(10*np.log10(x_pwr_avg), '.')

``````
``````

Out[3]:

[<matplotlib.lines.Line2D at 0x7f60543de198>]

``````
``````

In [4]:

tone = 10*np.log10(x_pwr_avg) > -25
windows = list()
in_tone = False
for j, t in enumerate(tone):
if not in_tone and t:
left = j
in_tone = True
elif in_tone and not t:
windows.append((left, j))
in_tone = False
tone_windows = [w for w in windows if w[1] - w[0] > 50]

``````
``````

In [5]:

def estimate_rice(x):
r = np.average(x) / np.std(x)
theta = 1
for iteration in range(100):
I0 = scipy.special.iv(0, theta**2/4)
I1 = scipy.special.iv(1, theta**2/4)
xi = 2 + theta**2 - np.pi/8 * np.exp(-theta**2/2) * ((2 + theta**2) * I0 + theta**2 * I1)**2
old_theta = theta
theta = np.sqrt(xi * (1 + r**2) - 2)
if np.abs(theta - old_theta) < 1e-14:
break
sigma = np.std(x) / np.sqrt(xi)
nu = np.sqrt(np.average(x)**2 + (xi - 2) * sigma**2)
return nu, sigma

``````
``````

In [6]:

tone_parameters = np.array([estimate_rice(np.abs(x[(w[0] + 1) * avg : (w[1] - 1) * avg])) for w in tone_windows])
K = 0.5*tone_parameters[:,0]**2/tone_parameters[:,1]**2
Omega = tone_parameters[:,0]**2 + 2*tone_parameters[:,1]**2

``````
``````

In [7]:

fig, ax1 = plt.subplots(figsize = (12, 6))
ax1.plot(10*np.log10(K), 'o-', label = 'K')
ax2 = ax1.twinx()
ax2.plot(10*np.log10(Omega), 'x-', color = 'C1', label = '\$\Omega\$')
fig.legend()
plt.title('Rician model parameters')
ax1.set_ylabel('K (dB)', color = 'C0')
ax2.set_ylabel('\$\Omega\$ (dB)', color = 'C1')
ax1.set_xlabel('Measurement')
ax1.set_xticks(np.arange(K.size));

``````
``````

``````
``````

In [8]:

def plot_correlation(x, noncoherent = False, span = 1000):
if noncoherent:
x = np.abs(x)
x = x - np.average(x)
corr = scipy.signal.correlate(x, x, mode = 'full')
corr /= scipy.signal.correlate(np.ones(x.size), np.ones(x.size), mode = 'full')
corr /= np.max(np.abs(corr[corr.size//2-span:corr.size//2+span]))
t = np.arange(-span, span)*1e-3
plt.plot(t, np.abs(corr[corr.size//2-span:corr.size//2+span]))

``````
``````

In [9]:

plt.figure(figsize = (10,6), facecolor = 'w')
for w in tone_windows:
plot_correlation(x[(w[0] + 1) * avg : (w[1] - 1) * avg])
plt.title('Signal self-correlation')
plt.xlabel('Time (s)')
plt.ylabel('Correlation');

``````
``````

``````
``````

In [10]:

plt.figure(figsize = (10,6), facecolor = 'w')
for w in tone_windows:
plot_correlation(x[(w[0] + 1) * avg : (w[1] - 1) * avg], noncoherent = True, span = 1000)
plt.title('Amplitude noise self-correlation')
plt.xlabel('Time (s)')
plt.ylabel('Correlation');

``````
``````

``````
``````

In [11]:

for sig in [x[(w[0] + 1) * avg : (w[1] - 1) * avg] for w in tone_windows]:
plt.figure()
bins = plt.hist(np.abs(sig), bins = 100)[1]
tone_parameters = estimate_rice(np.abs(sig))
K = 0.5*tone_parameters[0]**2/tone_parameters[1]**2
Omega = tone_parameters[0]**2 + 2*tone_parameters[1]**2
pdf = scipy.stats.rice.pdf(bins, tone_parameters[0]/tone_parameters[1], scale = tone_parameters[1])
pdf *= sig.size / np.sum(pdf)
plt.plot(bins, pdf)
plt.title('Amplitude distribution')
plt.legend(['Rice distribution fit', 'Measured'])
plt.xlabel('Amplitude')
plt.ylabel('Number of samples')

``````
``````

``````
``````

In [12]:

def plot_spectrum(x, N = 1024):
w = np.hanning(N)
fspan = 50*N//1024
f_even = np.abs(np.fft.fftshift(np.fft.fft(x[:x.size//N*N].reshape((-1,N)) * w), axes = 1)[:,N//2-fspan:N//2+fspan])**2
f_odd = np.abs(np.fft.fftshift(np.fft.fft(x[N//2:x.size//N*N-N//2].reshape((-1,N)) * w), axes = 1)[:,N//2-fspan:N//2+fspan])**2
f = np.zeros((f_even.shape[0] + f_odd.shape[0], f_even.shape[1]))
f[::2,:] = f_even
f[1::2,:] = f_odd
plt.figure(figsize = (12,8), facecolor = 'w')
plt.imshow(np.log10(f).T, extent = (0, f.shape[0]*N//2*1e-3, -fspan*1e3/N, fspan*1e3/N), aspect = 'auto', origin = 'bottom', vmin = -3, vmax = 2.5, interpolation = 'lanczos')
plt.ylabel('Frequency (Hz)')
plt.xlabel('Time (s)')
plt.title('Spectrum of beacon tones')
return f

``````
``````

In [13]:

all_tones = np.concatenate([np.concatenate((x[(w[0] + 1) * avg : (w[1] - 1) * avg], np.zeros(1000))) for w in tone_windows])
plot_spectrum(all_tones);

``````
``````

``````