In [1]:
import re
import glob
Tomaž Šolc
Jožef Stefan Institute
Jamova 39, 1000 Ljubljana, Slovenia
tomaz.solc@ijs.si
In the following, each method is marked with an abbreviation as follows:
In [2]:
Pdmin = 0.9
In [3]:
Pfa = 0.1
Given $P_{fa}$ we can determine $\gamma_0$ from $\gamma$ complementary cumulative density function (CCDF) for the given detector with no input signal.
To measure the distribution of $\gamma$, the RF output of the generator was switched off and $N_p$ measurements of $\gamma$ taken for each of the chosen detectors.
In [4]:
Np = 1000
Given these measurements we approximate the CCDF and determine the value of $\gamma_0$ by interpolation.
In [5]:
def get_ccdf(x):
xs = array(x)
xs.sort()
N = float(len(xs))
P = arange(N)/N
return xs, P
def get_gamma0(gammaN):
gammaN, Pd = get_ccdf(gammaN)
gamma0 = interp(1. - Pfa, Pd, gammaN)
return gamma0
For example, energy detector using $f_s=1\mathrm{MHz}$ and $N_s=25000\mathrm{samples}$ yields the following value for the threshold:
In [6]:
gammaN = loadtxt("../measurements/pd/sim/sim_micsoft_gaussian_noise_m100dbm_fs1mhz_Ns25ks_ed_off.dat")
gamma0 = get_gamma0(gammaN)
gammaNs, Pd = get_ccdf(gammaN)
plot(gammaNs, 1.-Pd)
plot([gamma0, gamma0], [0, 1], 'b--')
a = axis()
plot([a[0], a[1]], [Pfa, Pfa], 'b--')
title("CCDF")
xlabel("$\gamma_0$")
ylabel("$P_{fa}$")
grid()
In [7]:
print "gamma0 =", gamma0, "@ Pfa =", Pfa
In [8]:
def iterate_campaign(path):
for fn in glob.glob(path):
g = re.search("_m([0-9_]+)dbm\.dat$", fn)
if g:
Pg = -float(g.group(1).replace('_', '.'))
gamma = loadtxt(fn)
yield Pg, gamma
print "Pg =", array(sorted(Pg for Pg, gamma
in iterate_campaign("../measurements/pd/sim/sim_micsoft_gaussian_noise_m100dbm_fs1mhz_Ns25ks_ed_*.dat")))
If $P_g$ is the output power, then
$P_{in}$ = $P_g$ - $A$
Probability of detection $P_d$ is calculated for each $P_{in}$ and each method by counting the number of $\gamma$ measurements that were above the calculated threshold $\gamma_0$.
In [9]:
def get_campaign_g(path, gamma0):
Pg = []
Pd = []
for Pg0, gamma in iterate_campaign(path):
Pg.append(Pg0)
Pd.append(mean(gamma > gamma0))
Pg = array(Pg)
Pd = array(Pd)
Pga = Pg.argsort()
Pd = Pd[Pga]
Pg = Pg[Pga]
return Pg, Pd
def get_campaign(path, gamma0):
Pg, Pd = get_campaign_g(path, gamma0)
return Pg, Pd
def get_Pinmin(path, Pdth=Pdmin):
gamma0 = None
for fn in glob.glob(path):
if fn.endswith("_off.dat"):
gammaN = loadtxt(fn)
gamma0 = get_gamma0(gammaN)
break
Pin, Pd = get_campaign(path, gamma0)
Pinmin = interp(Pdth, Pd, Pin)
return Pinmin
Given minimum $P_d$, $P_{dmin}$, the minimum signal power $P_{inmin}$ for the spectrum sensing method can be determined from measurements using interpolation.
For example, for energy detector:
In [10]:
Pinmin = get_Pinmin("../measurements/pd/sim/sim_micsoft_gaussian_noise_m100dbm_fs1mhz_Ns25ks_ed_*.dat")
print "Pinmin =", Pinmin, "@ Pdmin =", Pdmin
In [11]:
Pin, Pd = get_campaign("../measurements/pd/sim/sim_micsoft_gaussian_noise_m100dbm_fs1mhz_Ns25ks_ed_*.dat", gamma0)
plot(Pin, Pd, 'o-')
plot([Pinmin, Pinmin], [0, 1], 'b--')
x1, x2, y1, y2 = axis()
plot([x1, x2], [Pdmin, Pdmin], 'b--')
xlabel("$P_{in}$ [dBm]")
ylabel("$P_d$")
grid()
In [12]:
method_ps = set()
for path in glob.glob("../measurements/pd/sim/sim_micsoft*.dat"):
g = re.search("ks_(.*)_m", path)
if g:
method_ps.add(g.group(1))
method_ps = sorted(method_ps)
In [13]:
def get_batch_Pinmin_permethod(path, Pdth=Pdmin):
Pinmin_permethod = {}
for method_p in method_ps:
path2 = path.replace("*", "%s_*" % (method_p,))
Pinmin_permethod[method_p] = get_Pinmin(path2, Pdmin)
return Pinmin_permethod
In [14]:
from matplotlib import cm
def plot_comparisson(Pinmin_permethod):
cmap = cm.get_cmap('gist_rainbow')
left = []
height = []
labels = []
colors = []
methods = ['ed', 'scf', 'cav', 'cfn', 'mac', 'mme', 'eme', 'agm', 'met']
left1 = 0.
for n, method in enumerate(methods):
if method == 'ed':
method_ps = [ 'ismtv', 'ed']
elif method == 'scf':
method_ps = [ "%s_Np%d" % (method, Np) for Np in [64, 128] ]
else:
method_ps = [ "%s_l%d" % (method, l) for l in xrange(5, 25, 5) ]
for m, method_p in enumerate(method_ps):
try:
Pinmin = Pinmin_permethod[method_p]
except KeyError:
continue
left.append(left1)
height.append(Pinmin)
label = method_p.upper()
label = re.sub("_L(\d+)", r" $L=\1$", label)
label = re.sub("_NP(\d+)", r" $N'=\1$", label)
labels.append(label)
c = cmap(float(n)/(len(methods)-1) + m*.007)
# Color
colors.append(c)
# BW
#colors.append(.85)
left1 += 1.
left1 += .9
height = array(height)
left = array(left)
# print
#figure(figsize=(13, 5))
# screen
figure(figsize=(12, 5))
#o = 100
#bar(left, -height-o, color=colors, bottom=o)
bar(left, height, color=colors)
for x, y in zip(left, height):
text(x+.4, y-1., "%.1f" % y, horizontalalignment="center", rotation=70, backgroundcolor='w')
xticks(left+.4, labels, rotation=60, horizontalalignment='right')
ylabel("$P_{in}$ [dBm] @ $P_d = 0.9$")
#ylim(100, 118)
ylim(-130, -110)
xlim(-.5, left1-.5)
tight_layout()
grid()
Lower detectable power is better
In [15]:
Pinmin_permethod_25ms = get_batch_Pinmin_permethod("../measurements/pd/sim/sim_micsoft_gaussian_noise_m100dbm_fs1mhz_Ns25ks_*.dat")
plot_comparisson(Pinmin_permethod_25ms)
title("$t=25.0$ms (simulation $f_s=1$MHz, $N_s=25$ksamples)");
savefig("figures/pin_min_comparison_25ms_sim.png", dpi=300)
#savefig("pin_min_comparison_25ms.eps")
None
With a higher sampling rate, performance seems to decrease. Probably because less of the signal is covered with 25.000 samples.
In [16]:
Pinmin_permethod_13ms = get_batch_Pinmin_permethod("../measurements/pd/sim/sim_micsoft_gaussian_noise_m100dbm_fs2mhz_Ns25ks_*.dat")
plot_comparisson(Pinmin_permethod_13ms)
title("$t=12.5$ms (simulation $f_s=2$MHz, $N_s=25$ksamples)");
savefig("figures/pin_min_comparison_13ms_sim.png", dpi=300)
#savefig("pin_min_comparison_13ms.eps")
None
In [17]:
Pinmin_permethod_10ms = get_batch_Pinmin_permethod("../measurements/pd/sim/sim_micsoft_gaussian_noise_m100dbm_fs10mhz_Ns100ks_*.dat")
plot_comparisson(Pinmin_permethod_10ms)
title("$t=10.0$ms (simulation $f_s=10$MHz, $N_s=100$ksamples)");
savefig("figures/pin_min_comparison_10ms_sim.png", dpi=300)
#savefig("pin_min_comparison_10ms.eps")
None
We compare the probability of false alarm $P_{fa}$ for different levels of noise. First, the threshold $\gamma_0$ is calculated for the same noise level as when comparing $P_d$ (-100 dBm). Then, the noise level is increased and $P_{fa}$ plotted versus noise level.
Noise level is shown compared to the reference -100 dBm noise level.
In [18]:
def get_noise_threshold(path, Pfath = 2*Pfa):
gammaN = loadtxt(path.replace("*", "m100_0dbm"))
gamma0 = get_gamma0(gammaN)
Pin, Pd = get_campaign(path, gamma0)
ratio_max = interp(Pfath, Pd, Pin) + 100
return ratio_max
def plot_noise_threshold(path):
gammaN = loadtxt(path.replace("*", "m100_0dbm"))
gamma0 = get_gamma0(gammaN)
Pin, Pd = get_campaign(path, gamma0)
ratio_max = get_noise_threshold(path)
print "ratio_max =", ratio_max
plot(Pin+100, Pd, 'o-')
x1, x2, y1, y2 = axis()
plot([x1, x2], [2*Pfa, 2*Pfa], 'b--')
plot([ratio_max, ratio_max], [0, 1], 'b--')
xlabel("$P_{noise}'/P_{noise}$ [dB]")
ylabel("$P_{fa}'$")
grid()
In [19]:
plot_noise_threshold("../measurements/pd/sim/sim_noise_fs2mhz_Ns25ks_ed_*.dat")
title("ED (simulation, $f_s=2$MHz, $N_s=25$ksamples)");
In [20]:
plot_noise_threshold("../measurements/pd/sim/sim_noise_fs2mhz_Ns25ks_eme_l10_*.dat")
title("EME L=10 (simulation, $f_s=2$MHz, $N_s=25$ksamples)");
In [21]:
plot_noise_threshold("../measurements/pd/sim/sim_noise_fs2mhz_Ns25ks_cam_Np128_*.dat")
title("SCF N'=128 (simulation, $f_s=2$MHz, $N_s=25$ksamples)");
In [21]: