Monarchs and Matched Filtering¶
V. Hunter Adams¶
import numpy
import matplotlib.pyplot as plt
from scipy.interpolate import InterpolatedUnivariateSpline
from scipy.integrate import odeint
from IPython.display import Latex
from IPython.display import Image
import binascii
from numpy.linalg import pinv
import random
import scipy
%matplotlib inline
A Monarch in ISM Band¶
The Monarch transmits at 915 MHz, at a wavelength of approximately 32.8 cm. Let us consider the case of two antennas (transmit and receive antennas) in free space separated by a distance $R$.
Derivation of Friis Equation¶
Let us assume that $P_T$ Watts of total power is delivered to the transmit antenna, which (for the moment) is assumed to be omnidirectional and lossless. Furthermore, we'll assume that the receive antenna is in the far field of the the transmit antenna (a safe assumption for transmissions from orbit). As the signal propagates spherically out from the transmit antenna, the power density (watts per square meter) of the plane wave decreases. By the time the signal reaches the receive antenna at a distance $R$ away, the power density is given by:
Any losses and directionality of the transmit antenna can be absorbed by a gain $G_T$. A transmit gain greater than one for a lossless antenna means that it is transmitting in a preferrred direction, and that direction is towards the receive antenna. Augmenting the above equation:
Now consider the receive antenna. The aperature (i.e. effective area or receiving cross section) of an antenna is a measure of how effective an antenna is at receiving the power of radio waves. It is the area, oriented perpendicular to the direction of an incoming radio wave, that would intercept the same amount of power from that wave as is produced by the antenna receiving it. We can therefore augment the equation again to get received power:
The aperature for any antenna can also be expressed as:
(derivation based on thought experiment involving radiating resistive load in thermal equilibrium with an antenna, haven't grok'd this yet)
Rewriting again:
This is the Friis Transmission Formula. In the above equation, power is measured in linear units (watts). If, however, we convert to logarithmic units (decibels), we get the following:
def todBw(mW):
return 10*numpy.log10(mW) - 30
def todBm(mW):
return todBw(mW) + 30
def toMilliWatts(dBm):
return 10.**((dBm)/10.)
todBw(25)
Or, in decibels:
Reflecting on the Friis Equation¶
From the above equation, it can be seen that path loss increases as wavelength decreases (or as frequency increases). This is why cell phones generally operate at less than 2 GHz. The received power depends on the power of the transmitter, the gain of the transmitter, the gain of the receiver, the wavelength being used, and the distance between antennae. Let us now substitute in some realistic values for the Sprite.
Received Power from Sprite¶
The Monarch has a $25\text{mW} = -16 dBW$, isotropic ($G_T=0$) transmitter that communicates on the 32.8 cm ($\lambda=0.328$) band from orbit ($R=500000m$). On the receiving end, we'll assume a handheld yagi with a gain of approximately $G_R=7 dB$. Substituting and solving, we find:
Pt = 25 #mW
Gt = 0. #dBm (isotropic)
lam = 0.328 #meters
R = 500000. #meters
Gr = 7 #dBm
Coding up the Friis Equation:
def friisEquation(Pt=Pt, Gt=Gt, Gc=0, Gr=Gr, lam=lam, R=R):
return todBm(Pt) + Gt + Gc + Gr + todBm((lam/(4*numpy.pi*R))**2.)
friisEquation()
The noise power is given by the sum of the noise temperature and feed losses. The black body thermal noise is radiated from the environment, and the feed losses come from within the system.
$K_B$ is the Botzmann constant, $T=150K$ is the background temperature, and $B=64 kHz$ is the receiver bandwidth. $F_N$ is the receiver noise figure, assumed to be $9dBW$ (coming largely from the low noise amplifier). Solving:
Fn = 39. #dBm
T = 150. #Kelvin
B = 50000. #Bandwidth (Hz)
Kb = 1.38064852e-23 #Boltzmann (SI Units)
Code up the noise equation:
def noiseEquation(Fn=Fn, T=T, B=B, Kb=Kb):
return Fn + todBm(Kb*T*B)
noiseEquation()
This gives a signal to noise ratio (SNR) of approximately:
def getSNRdBm(Gc=0):
return friisEquation(Gc=Gc) - noiseEquation()
getSNRdBm()
If we were to represent received power and noise power in watts and find their ratio, the above result means that the ratio would be:
def getSNRRatio(Gc=0):
return toMilliWatts(getSNRdBm(Gc=Gc))
getSNRRatio()
Just how bad is this?¶
Consider the following signal, which varies in amplitude between -1 and 1:
Create a signal of the designated length:
def createSignal(totransmit):
totransmit=str(totransmit)
test= [ord(c) for c in totransmit]
test1 = []
for i in test:
test1.extend(['{0:08b}'.format(i)])
test2 = ''
for i in test1:
test2 += i
data = ' '
for i in test2:
data+=i+' '
signal= numpy.fromstring(data, sep=' ')
return (2*signal-1)
signal = createSignal('Greetings Earthlings')
Based on the variance of that signal and the signal to noise ratio, find the standard deviation of the additive Gaussian noise:
def getNoiseSTD(signal, Gc=0):
return numpy.sqrt(1./getSNRRatio(Gc=Gc))
Generate a noise vector by drawing from a Gaussian with the appropriate standard deviation:
def makeNoise(noiseSTD, signal):
return noiseSTD*numpy.random.randn(len(signal))
Create a function that compares the transmitted to the received signal:
def runMonteCarlo(runs=100, ylim=4.5, Gc=0):
noiseSTD = getNoiseSTD(signal, Gc=Gc)
for i in range(runs):
noise = makeNoise(noiseSTD, signal)
plt.plot(signal+noise, alpha=0.1)
plt.plot(signal, 'r',lw=2, label='Signal', alpha=1)
plt.title('Received Signal')
plt.legend(loc='best')
plt.ylim([-ylim,ylim])
plt.show()
With all noise in view:
runMonteCarlo()
Within a tighter range:
runMonteCarlo(ylim=1.5)
What's the Bit Error Rate?¶
From the above SNR, we can solve for the variance of the additive Gaussian noise:
thresh=0.5
def drawHistos(draws=1000, Gc=0):
plt.hist(getNoiseSTD(signal, Gc=Gc)*numpy.random.randn(draws)-1, bins=50,
alpha=0.5, label='p(-1)')
plt.hist(getNoiseSTD(signal, Gc=Gc)*numpy.random.randn(draws)+1, bins=50,
alpha=0.5, label='p( 1)')
plt.plot(numpy.ones(70)-thresh, numpy.arange(0,70,1), 'b--', label='threshold', lw=2)
plt.plot(numpy.ones(70)-(1+thresh), numpy.arange(0,70,1), 'g--', lw=2)
plt.title('Additive Gaussian Noise Distribution')
plt.legend(loc='best')
plt.show()
drawHistos(Gc=0)
Let us assume a thresholding of 0.5. That is to say, a signal above 0.5 is considered a 1, and a signal below -0.5 is considered a -1. These thresholds are plotted above. Let us furthermore assume that there is an equal probability of a 1 or a -1 being transmitted. The probability of a bit error is then given by:
PATH = "./"
Image(filename = PATH + "tree.png", width=600, height=600)
The probability of transmitting a 0, the probability of transmitting a 1, the mean value for 0 transmissions, and the mean value for 1 transmissions:
p0 = 0.5
p1 = 0.5
mu0 = -1.
mu1 = 1.
For a particular threshold, find the z-values:
We can find the z-values for each conditional probability:
def getZvals(sigma):
return (thresh-mu0)/sigma, (-thresh-mu1)/sigma
We can plot the z-values on the normalized Gaussian distribution to see the sections that we'll be integrating:
def showZVals(mu=0, Gc=0):
sigma = getNoiseSTD(signal, Gc=Gc)
z1, z2 = getZvals(sigma)
x = numpy.linspace(-5,5,1000)
curve = (1./(sigma*numpy.sqrt(2*numpy.pi)))*numpy.exp((-(x-mu)**2.)/(2*sigma**2.))
spl = InterpolatedUnivariateSpline(x, curve)
y1 = numpy.linspace(0, spl(z1), 1000)
y2 = numpy.linspace(0, spl(z2), 1000)
plt.plot(x, curve)
plt.plot(numpy.ones(len(y1))*z1, y1, 'r--', label='z-vals')
plt.plot(numpy.ones(len(y1))*z2, y2, 'r--')
plt.legend(loc='best')
plt.title('Z-Distribuion')
plt.show()
print('Z1: '+str(z1))
print('Z2: '+str(z2))
showZVals(Gc=0)
Write a function to integrate for the conditional probabilities:
def derivs(y, t, sigma, mu):
dydt = [(1./(sigma*numpy.sqrt(2*numpy.pi)))*numpy.exp((-(t-mu)**2.)/(2*sigma**2.))]
return dydt
def getIntegral(zval, Gc=0):
t = numpy.linspace(-9*getNoiseSTD(signal, Gc=Gc), zval, 5000)
y0 = [0.]
sigma = getNoiseSTD(signal, Gc=Gc)
sol = odeint(derivs, y0, t, args=(sigma, 0.))
return sol[:,0]
Make sure this is correct by verifying that the curve is normalized:
def testIntegral(Gc=0):
test=getIntegral(10)
sigma = getNoiseSTD(signal, Gc=Gc)
plt.plot(numpy.linspace(-9*sigma, 9*sigma, 5000), test)
plt.plot(numpy.linspace(-9*sigma, 9*sigma, 5000), numpy.ones(len(test)), 'r--')
plt.xlim([-10, 10])
plt.show()
testIntegral()
Use the above to write a function for the conditional probability:
def conditionalProb(sigma, Gc=0):
z = min(getZvals(sigma))
prob = getIntegral(z, Gc=Gc)[-1]
return prob
We can now find the total probability of a bit error:
def bitErrorRate(Gc, signal=signal):
sigma = getNoiseSTD(signal, Gc=Gc)
perr = conditionalProb(sigma, Gc=Gc)
return p0*perr + p1*perr
bitErrorRate(Gc=0)
This gives a total probability of bit error of:
This indicates that 26.7 percent of transmitted bits will be corrupted. This is an intolerable amount of bit error, which can be improved by adding coding gain to the transmission.
How much Coding Gain is Required?¶
Let's plot the bit error rate as a function of coding gain:
def BERvsGc():
Gcs = numpy.linspace(0, 16, 100)
BERs = numpy.zeros(len(Gcs))
for i in range(len(Gcs)):
BERs[i] = bitErrorRate(Gc=Gcs[i])
plt.plot(Gcs, BERs*100, label='bit error rate')
plt.xlabel('Gc (dBm)')
plt.ylabel('Bit Error Rate (percent)')
plt.title('Bit Error Rate vs. Coding Gain')
plt.yscale('log')
plt.legend(loc='best')
plt.show()
BERvsGc()
The Shannon-Hartley Theorem (Quick Sidestep)¶
Introducing the Theorem¶
The noisy channel coding theorem states that for any given degree of noise contamination of a communication channel, discrete data can be communicated error-free up to a computable maximum rate through the channel. If we apply this to the case of continuous-time analog communications in the presence of Gaussian noise, we arrive at the Shannon-Hartley Theorem.
This theorem establishes the channel capacity for the communication link, which is a bound on the maximum amount of error-free information that can be transmitted per unit time with a specified bandwidth. It is assumed that the signal power is bounded, and that the Gaussian noise is characterized by a known power or power spectral density.
Applied to the Sprite¶
If we substitute the Sprite's bandwidth and SNR into the above equation, we find:
def getShannon():
return B*numpy.log2(1+getSNRRatio())
getShannon()
This means that the maximum amount of error-free information that can be communicated via this noisy channel is 25 kilobits per second. If a line rate of less than 25 kbps is used, there exists a coding technique (involving error correction) which allows the probability of error at the receiver to be made arbitrarily small. The theorem tells us nothing about what this coding technique is, however.
Improving the Signal to Noise Ratio¶
A matched filter is obtained by correlating a known signal (or template) with an unknown signal to detect the presence of the template in the unknown signal. A matched filter is the optimal linear filter for maximizing the signal to noise ratio in the presence of additive stochastic noise. A proof of this is given below (from Wikipedia):
Consider the observed signal $x$, which is composed of the desirable signal $s$ and additive noise $\nu$:
We seek a filter $h$ of the form shown below that will maximize the signal to noise ratio.
Let us define the covariance matrix for the noise (which is Hermitian):
Consider the output, $y$, which is the inner product of our filter and the received signal:
To be clear, we're considering some linear filter $h$. The output of our filter is the inner product of $h$ and the received signal, which is composed of the desired signal $s$ and the noise $\nu$. In the above equation, we consider the filter output from the signal, $y_s$, and the filter output from the noise, $y_{\nu}$. To maximize the signal to noise ratio, we want to maximize the ratio below:
Expanding the denominator:
Rewriting the SNR:
We now pull a trick, rewriting the numerator to lead ourselves to the Cauchy-Schwartz inequality:
We want to find an upper limit on the above expression. First, we recognize a form of the Cauchy-Schwartz Inequality:
Applying this to the above expression:
Simplifying the term on the right:
This limit is achieved when:
Yielding the optimal filter:
which makes sense intuitively. The linear filter is essentially a sliding dot product, and the above expression states that the SNR is best when the filter (the signal you're looking for) is parallel to the desired signal. One can choose $\alpha$ in accordance to the desired power of the filter output in response to noise (often unity).
Solving for $\alpha$:
which yields the normalized filter:
SNR and Coding Gain¶
Let's consider this equation from above one more time:
The term on the right places an upper bound on the signal to noise ratio in terms of the noise covariance (something we have no control over) and the desired signal (something we do have control over). We are sending digital signals composed of 1's and 0's (converted to 1's and -1's for convenience). Let's consider how the bound on SNR changes with signal length:
Signal Length 1 (in $P_R=1$ Units)¶
Signal Length 2¶
$\vdots$¶
Signal Length N¶
What this means¶
For a signal of length 1, the SNR is given by:
For a signal of length $N$:
In dBm:
With matched filtering, we have effectively increased the gain of the transmitted signal. The gain is related to the length of the filter as shown above. To be clear:
where $N$ is the length of the code. By adding 10 dBm of coding gain, we can see that the conditional probabilities that lead to bit error become far smaller, and the z-values increase in magnitude:
drawHistos(Gc=10*numpy.log10(16))
showZVals(Gc=10*numpy.log10(16))
An Example¶
We want to transmit the message 'Greetings Earthlings.' Begin by converting this string to ASCII, and then to binary:
signal = createSignal('Greetings Earthlings')
print((signal+1)/2)
plt.plot(signal)
plt.ylim([-1.5, 1.5])
plt.title('Message in Binary')
plt.show()
Let's consider a matched filter of length 16, which will add $10\cdot \text{log}_{10}(16) \approx 12$dB of gain. We'll use the codes shown below (justified in a later section):
code0 = numpy.array([numpy.ones(32)])
code1 = -1*numpy.array([numpy.ones(32)])
Write a function that creates a message, encodes it, and noisifies it:
def encode(signal, code0, code1):
encoded_message = numpy.array([[]])
for i in signal:
if i==1:
encoded_message = numpy.hstack((encoded_message, code1))
else:
encoded_message = numpy.hstack((encoded_message, code0))
encoded_message = encoded_message[0]
return encoded_message
Write a function to observe the encoded signal:
def showEncoded():
encoded_message = encode(signal, code0, code1)
sigma = getNoiseSTD(encoded_message, Gc=0)
noise = makeNoise(sigma, encoded_message)
plt.plot(encoded_message+noise, label='encoded signal+noise');
plt.plot(encoded_message, 'r', lw=2, label='encoded signal')
plt.title('Encoded and Received Signals')
plt.xlabel('Sample Number')
plt.ylabel('Received Power')
plt.legend(loc='best')
plt.show()
showEncoded()
Write a function to decode the noisy message:
def decode(encoded_noisy_message, code0, code1):
zero_dot_products = []
one_dot_products = []
chiplength = len(code0[0])
messagelength = len(encoded_noisy_message)
sigma = getNoiseSTD(signal)
scale = 1./(sigma**2. * numpy.sqrt((1./sigma**2.)*chiplength))
for i in range(0, messagelength - chiplength):
zdot = numpy.dot(numpy.array([encoded_noisy_message[i:i+chiplength]]), scale*code0.T)
odot = numpy.dot(numpy.array([encoded_noisy_message[i:i+chiplength]]), scale*code1.T)
zero_dot_products.extend([zdot[0][0]])
one_dot_products.extend([odot[0][0]])
return numpy.array(one_dot_products), numpy.array(zero_dot_products)
Write a function that takes a decoded message and finds the difference in the dot products:
def testDecode():
encoded_message = encode(signal, code0, code1)
messagelen = len(encoded_message)
chiplen = len(code0[0])
sigma = getNoiseSTD(encoded_message, Gc=0)
noise = makeNoise(sigma, encoded_message)
encoded_noisy_message = encoded_message+noise
decoded1, decoded0 = decode(encoded_noisy_message, code0, code1)
difference = decoded1 - decoded0
plt.plot(encoded_noisy_message, label='Received Message', alpha=0.5)
plt.plot(decoded1 - decoded0, label='Decoded Message')
plt.grid('on')
plt.legend(loc='best')
plt.title('Raw and Filtered Signal')
plt.show()
return difference
decoded_message = testDecode();
Write a soft decoder to recover the transmitted message:
def softDecode(filtered_message, code0, code1):
message = []
messagestr = ''
messagelen = len(filtered_message)
chiplen = len(code0[0])
for i in numpy.arange(0, messagelen-1, chiplen):
message.extend([int(filtered_message[i]/abs(filtered_message[i]))])
message.extend([int(filtered_message[-1]/abs(filtered_message[-1]))])
for i in message:
messagestr += str(int((i+1)/2))
return numpy.array(message), messagestr
Write a function for converting back to ASCII --> characters.
def decodeCharacters(data):
message = ''
counter = 0
while counter < (len(data) - len(data)%8):
message+=chr(int(data[counter:counter+8],2))
counter+=8
return message
Decode, and show results:
soft1, soft2 = softDecode(decoded_message, code0, code1)
print('Received Message:')
print(soft2)
print('\n')
print('Difference between Sent and Received Messages:')
diff = soft1 - signal
print(diff)
print('\n')
print('Number of Errors:')
print(str(numpy.count_nonzero(diff))+'/'+str(len(signal))+' bits')
print('\n')
print('Received Message (converted back to characters)')
print(decodeCharacters(soft2))
Bayesian Analysis of PRN Codes¶
Suppose our PRN codes are of length N. As the raw binary stream comes in, we can dot each N bits in the stream with each of the PRN codes. Because all of the codes are orthogonal, a perfect transmission of a particular PRN code will yield a dot product of N with the correct PRN, and a dot product of zero with all other PRN's. As shown in the previous sections, however, the bit error rate for individual bits is quite high for the case of the Sprite. This means that a transmission will have some degree of corruption, and will likely have a dot product $<N$ with the correct PRN code and a dot product $>0$ with some of the incorrect PRN codes. Furthermore, it's not always the case that the Sprite is transmitting information. Random noise will have some power in some of the PRN codes (and, on rare occasion, it may have a great deal of power along a particular PRN code).
All of the PRN codes are orthogonal. If they are of length N, then we can construct N orthogonal codes using Hadamard recursion. If we normalize all of these PRN codes, we can think of them as the basis vectors for a frame that is rotated relative to the frame of the incoming data. By normalizing all of the PRN codes, stacking them in columns, and then combining those columns into a matrix, we form a rotation matrix between the data frame and the frame for which the PRN codes are simply the columns of an identity matrix. This is useful because it makes it very easy to see how much power a particular string of bits has in each PRN code.
A string of bits may have power in a particular PRN code because:
- The Sprite transitted that PRN code
- The Sprite transmitted a different PRN code, but a bit error put some power into that particular PRN code
- The Sprite did not transmit anything, and random noise put some power into that particular PRN code
Via a nonlinear transformation, we can transform each of these vectors (now represented in the PRN frame) such that their magnitude in each direction (along each PRN code) is precisely the probability that a Sprite transmitted that PRN code.
This is useful because it means the output of the filter (which is now nonlinear and time-varying) is a meaningful quantity. Rather than looking for spikes associated with dot products, we can say for every string of $N$ bits precisely how confident we are that that string of bits came from a particular Sprite vs. any other Sprite vs. noise.
def hadamardRecursion(n):
if n == 1:
return numpy.array([[1.]])
else:
top = numpy.hstack((hadamardRecursion(n-1), hadamardRecursion(n-1)))
bottom = numpy.hstack((hadamardRecursion(n-1), -1.*hadamardRecursion(n-1)))
return numpy.vstack((top, bottom))
$p(j|PRNx)$¶
Consider the case where PRNx, of length $N$, is being transmitted. Suppose that the value of the dot product associated with PRNx is $j$. Find the conditional probability $p(PRNx|j)$.
In order to get the value $j$, $n$ bit flips are required, where $n$ is:
In order to get $j$ for the dot product, any $n$ of the $N$ bits must flip. It doesn't matter which $n$. The number of ways that one can choose $n$ bits from the $N$ bits is given by:
All of these potential solutions have the same probability:
Where $B$ is the bit error rate. Since any of these possibilities will do, to get the probability of any of them happening we sum the probability of each of them occuring:
def pJgivenPRNx(N, j, B):
first_term = numpy.math.factorial(N)/(numpy.math.factorial(int((N-j)/2.))*numpy.math.factorial(int(N - ((N-j)/2.))))
second_term = B**((N-j)/2.)
third_term = (1-B)**(N - ((N-j)/2.))
return first_term*second_term*third_term
Write a function to visualize this conditional probability for all possible dot products $j$ as the bit error rate is varies between 0 and 1:
def testpJgivenPRNx(N):
berr = numpy.arange(0., .5025, 0.025)
for j in berr:
test = []
for i in numpy.arange(-N, N+1):
test.extend([pJgivenPRNx(N, i, j)])
plt.plot(numpy.arange(-N, N+1), test, '-', label='Bit Error: '+str(j))
plt.xlim([-N,N])
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0., ncol=3)
plt.ylabel('$p(j|PRNx)$')
plt.xlabel('j')
plt.title('$p(j|PRNx)$')
plt.show()
testpJgivenPRNx(128)
And the converse:
def pNotJgivenPRNx(N, j, B):
return 1. - pJgivenPRNx(N, j, B)
$p(j|\overline{PRNx})$¶
The other possibility is that no PRN code was being transmitted and random noise ended up with some power along PRNx, or that another PRN code was being transmitted. These cases should be considered slightly differently, since we no longer need to take into account bit error rate. Each bit has an equal probability of being a 1 or a -1. Because all of the PRN codes are composed of 1's and -1's, this means that each random bit has an equal probability of increasing the dot product with a particular PRN code or of decreasing the dot product with a particular PRN code. Rather than think about this case in terms of 1's and -1's, it makes sense to consider it in terms of constructive and destructive bits.
The total number of random arrangements of -1's and 1's in a bitstring of length $N$ is given by:
Any one particular arrangement has an equal probability of any other of being randomly drawn. For a particular arrangement $a_i$, the probability is:
To get dot product $j$ with PRNx, we require that the differnce between constructive bits (bits that increment the dot product by one) and destructive bits (bits that decrement by one) to be $j$. Every bit will either be constructive or destructive and has equal probability of being either. Denoting all constructive bits as c and all destructive bits and d, we have the following system of equations:
Solving for $c$:
So, given dot product $j$, we require $c$ constructive bits. How many ways are there to choose $\frac{j+N}{2}$ constructive bits from the $N$ total bits? This is an n choose k problem:
Each of these $k$ arrangements has equal probability of occuring, so the probability of any of them occuring is the sum of all of their individual probabilities:
def pJgivenNotPRN(N, j):
return numpy.math.factorial(N)/((2.**N) * numpy.math.factorial(int((j+N)/2.)) * numpy.math.factorial(int(N - ((j+N)/2.))))
Write a function to visualize the conditional probability as the bit error is varied:
def testpJgivenNotPRN(N):
test = []
for i in numpy.arange(-N, N+.1, .1):
test.extend([pJgivenNotPRN(N, i)])
plt.plot(numpy.arange(-N, N+.1, .1), test, label='$p(j|Not PRN)$')
plt.xlim([-N,N])
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.ylabel('$p(j|Not PRN)$')
plt.xlabel('j')
plt.title('$p(j|Not PRN)$')
plt.show()
testpJgivenNotPRN(16)
(Note, independent of bit error rate). And the converse:
def pNotJgivenNotPRN(N, j):
return 1. - pJgivenNotPRN(N, j)
Bayes' Theorem¶
At this point, we've found each of the conditional probabilites.
The quantity that we're interested in is the conditional probability of PRNx having been sent given the received signal $j$. Bayes' Theorem provides that quantity:
All of these are quantities for which we've already solved. Coding it up, and assuming equal probability of a PRN code being sent or not sent $(p(PRNx) = p(\overline{PRNx}$):
def pPRNxgivenJ(N, j, B, NPRN):
numerator = pJgivenPRNx(N, j, B)
denominator = pJgivenPRNx(N, j, B) + pJgivenNotPRN(N, j)
return numerator/denominator
Write a function to test this out for varying dot products and bit errors:
def testpPRNxgivenJ(N):
berr = numpy.arange(0.025, 1., 0.025)
for j in berr:
test = []
for i in numpy.arange(-N, N, .1):
test.extend([pPRNxgivenJ(N, i, j, N)])
if j<=0.5:
plt.plot(numpy.arange(-N, N, .1), test, label='Bit Error: '+str(j))
else:
plt.plot(numpy.arange(-N, N, .1), test, label='Bit Error: '+str(j), alpha=0.5)
plt.xlim([-N,N])
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0., ncol=3)
plt.ylabel('$p(PRNx|j)$')
plt.xlabel('j')
plt.title('$p(PRNx|j)$')
plt.ylim([-0.1, 1.1])
plt.show()
testpPRNxgivenJ(150)
Building A Filter Around This¶
The classic matched filter is a linear time invariant filter that outputs magnitudes of dot products. The filter that I want will instead output the conditional probability of a particular sequence of bits being from any particular PRN code. This means that the filter will time varying and nonlinear. Essentially, we're pumping the output of the linear filter through a nonlinear function.
def getPRNMatrix(N, NPRN):
mat = numpy.zeros((NPRN, N))
for i in range(NPRN):
mat[i,:] = numpy.random.randint(0, 2, N)*2 - 1
return mat
We'll have a number of Sprites all talking at random, potentially overlapping intervals. It'll be up to the filter to let us know which Sprites are talking, and how confident it is that what it's receiving is indeed a signal. Begin by writing a function that will return a time-domain signal with randomly dispersed signals and silences and the proper amount of noise.
def makeRandomSignal(NPRN, numsigs, overlays):
integers = numpy.zeros((2*numsigs, overlays))
for i in range(len(integers[:,0])):
if i%2==0:
integers[i,:] = random.sample(range(NPRN+1), overlays)
return integers.astype(int)
Write a function that takes this signal and represents it with PRN codes:
def signalToPRN(N, NPRN, integer_signal):
prn_matrix = getPRNMatrix(N, NPRN)
rows = len(integer_signal[:,0])
cols= len(integer_signal[0,:])
signal = numpy.zeros((N*rows, cols*2 + 4))
for i in range(rows):
signal[i*N,cols:-4] = integer_signal[i,:]
for j in range(cols):
if integer_signal[i,j] == 0:
signal[i*N:(i+1)*N,j] = numpy.zeros(N)
else:
signal[i*N:(i+1)*N,j] = prn_matrix[integer_signal[i,j]-1,:]
for i in range(1,cols):
signal[i:,i] = signal[0:-i,i]
signal[0:i,i] = numpy.zeros(i)
signal[i:,i+cols] = signal[0:-i,i+cols]
signal[0:i,i+cols] = numpy.zeros(i)
noisestd = getNoiseSTD(1, Gc=0)
signal[:,-4] = numpy.sum(signal[:,cols:-2], axis=1)
signal[:,-3] = numpy.sum(signal[:,0:cols], axis=1)
signal[:,-2] = makeNoise(noisestd, signal[:,0])
signal[:,-1] = numpy.sign(numpy.rint(numpy.sum(signal[:,-3:-1], axis=1)))
# signal[:,-1] = signal[:,-3]
return signal, prn_matrix
Create a function to examine the raw signal:
def viewSignaltoPRN(N, NPRN, numlong, numstacked):
integers = makeRandomSignal(NPRN, numlong, numstacked)
test, mat = signalToPRN(N, NPRN, integers)
plt.plot(test[:,-1], label='Received Signal+Noise')
plt.plot(test[:,-3], label='Transmitted Signal')
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
print(integers)
plt.show()
viewSignaltoPRN(32, 5, 10, 1)
Write a function to perform the matched filtering, returning the dot product with each PRN code:
def matchedFilter(N, NPRN, signal, prn_matrix):
data = numpy.array([signal[:,-1]]).T
dot_prods = numpy.zeros((NPRN, len(data.T[0])-N))
for i in range(0, len(data.T[0])-N):
dot_prods[:,i] = numpy.dot(prn_matrix, data[i:i+N]).T[0]
return dot_prods, signal[:,-4]
Write a function that takes these dot products and returns the conditional probability:
vectorizedProbability = numpy.vectorize(pPRNxgivenJ, otypes=[numpy.float])
Put the pieces together:
def runSimulation(N, NPRN, numsigs, overlays, B):
integer_signal = makeRandomSignal(NPRN, numsigs, overlays)
signal, prn_matrix = signalToPRN(N, NPRN, integer_signal)
dot_prods, truths = matchedFilter(N, NPRN, signal, prn_matrix)
probability = vectorizedProbability(N, dot_prods, B, NPRN)
return probability, truths, integer_signal
test1, test2, test3 = runSimulation(150, 50, 5, 1, 0.18)
probs = []
dexs = []
times = []
for i in range(len(test1[:,0])):
for j in range(len(test1[0,:])):
if test1[i,j] > 0.01:
probs.extend([test1[i][j]])
dexs.extend([i])
times.extend([j])
fig, ax = plt.subplots()
ax.scatter(times, probs, s=.75)
for i, txt in enumerate(dexs):
ax.annotate(txt+1, (times[i], probs[i]))
plt.title('Decoded PRN Codes')
plt.ylabel('$p(PRN|j)$')
plt.xlabel('Time')
plt.ylim([0, 1.1])
plt.show()
print('Sequence of transmitted PRN codes: ')
for i in test3:
print(str(i))
