EE 120 Lab 4: Orthogonal Signaling

v1 - Spring 2019: Dominic Carrano, Sukrit Arora

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

Background

Digital Communication

In communication systems, our model consists of a sender and receiver communicating over a channel. Broadly speaking, there are two types of communication: analog and digital. Orthogonal Signaling is a method for digital communication. Unliked analog communication, where the data we actually want to send is a continuous-time signal, in digital communication, we want to send bits, which are discrete. For example, we may want to send a text file by first encoding it into its binary ASCII values. In the case of Amplitude Modulation (AM), an analog communication method, we saw that we needed to modulate the signal we wanted to send up to a higher frequency range before actually sending it. Similarly with digital communication, we cannot send raw bits themselves, and ultimately must convert our bits into a continuous-time signal before we can send it. Orthogonal signaling is one technique for achieving this.

Modulation

The sender communicates the bits by splitting up the data into fixed-sized chunks (for example, 3-bit chunks), modulating each, and sending the entire modulated signal. For an individual chunk, the sender will choose the signal $x(t)$ to send. If we use $k$ bit chunks, we'll need $2^k$ different signals to choose from, one per possible bit pattern. In orthogonal signaling, the choice of $x(t)$ is made from a family of orthogonal functions, which have nice properties for decoding and noise analysis as you saw in the homework. An orthogonal signaling scheme relies on two choices: $k$, the number of bits ber chunk, and $T$ (which corresponds to frequency $f_c = 1/T$), the duration of a continuous-time signal in the basis.

Specifically, if we use $k$ bits and a signal duration (per chunk) of $T$, then the basis for our Orthogonal Signaling scheme is:

\begin{align*} x_1(t) &= \cos\left(\frac{2\pi}{T} t\right) \\ x_2(t) &= \cos\left(\frac{2\pi}{T} 2t\right) \\ x_3(t) &= \cos\left(\frac{2\pi}{T} 3t\right) \\ &\;\vdots \\ x_m(t) &= \cos\left(\frac{2\pi}{T} mt\right) \\ &\;\vdots \\ x_{2^k}(t)&= \cos\left(\frac{2\pi}{T} 2^k t\right) \end{align*}

We choose the $m^\text{th}$ signal from this basis to represent the $m^\text{th}$ bit pattern.

For example, if $k=2$, then we have four possible bit patterns, $00, 01, 10, 11$, which are mapped to $T$ seconds of the signals $x_1(t), x_2(t), x_3(t), x_4(t)$. For $T = 1$, this basis of four signals looks like:

Just as you can see from the analytical expressions above, the plots suggest that the higher the binary value of the chunk of bits, the higher the frequency of its cosine - $x_4(t)$ is has four times the frequency of $x_1(t)$. Note that each signal completes an integer number of periods over $[0, T]$ and that the $m^\text{th}$ signal completes exactly $m$ periods over $[0, T]$.

Orthogonality of the Basis

We could prove the orthogonality of the basis by taking the inner product of our cosines in the time domain, but this method is a little involved. Instead, we can prove the orthogonality of the fourier transform of the cosines, and use the fact that the fourier transform preserves orthogonality.

In the frequency domain, we see that the proof is very straightforward. Each cosine is pair of deltas symmetric about 0. Now, we see that if the cosines have a different frequencies, the deltas will not overlap and the inner product will be 0. If they do overlap, they have the same frequency, and will have a non-zero inner product!

If you worked it out, you will find that:

$$\langle x_m(t), x_n(t) \rangle = \begin{cases}\displaystyle\frac{T}{2} & n = m \\[0.1in] 0 & n \neq m\end{cases}$$

where

$$\langle x(t), y(t) \rangle = \int_0^T x(t)y(t) dt$$

which, conveniently is just the continuous-time analogue of the dot product - the sum/integral of the pointwise product of the signals! In Fourier analysis, we'd be using a slightly different dot product, where $y(t)$ would be replaced by its complex conjugate, but here we restrict ourselves to real signals, hence the distinction.

Q1a: Orthogonal Signaling - Modulation

Implement modulate, which takes in a bitstring (the data), the number of bits per chunk k, and the length T in seconds of each signal (which is also the fundamental period), and returns the result of applying the orthogonal signaling modulation.

Since T is in seconds, you will have to do a little math to figure out the number of samples it corresponds to. Assume a sampling frequency of 44100 Hz, that is 44100 samples per second, in doing your conversion. This is a very common sampling frequency in signal processing, and used by many different applications including audio formats such as MP3.

fs = 44100 # sampling frequency, in Hz

def mth_basis_signal(m, T):
    samp_per_chunk = int(fs * T)
    t = np.linspace(0, T, samp_per_chunk)
    fc = 1 / T
    return np.cos(2 * np.pi * fc * m * t)

def modulate(bits, k, T):
    """
    Returns the modulated version of bits using an Orthogonal Signaling scheme
    with k bits per chunk and T second signals.
    
    Parameters:
    bits - The bitstring of data.
    k    - The number of bits per chunk. You may assume len(bits) % k == 0.
    T    - The duration, in seconds, of each individual modulated chunk.
    """
    ## TODO ##

To test your implementation of modulate, run the block below which compares the output of your function for various inputs. If the cell successfully runs without any AssertionErrors, you are passing all the tests.

# Tests involving "len" check that your sampling math is correct; others
# verify that the actual modulated waveform your function returns is correct.
bits1 = '0000'
T1 = 1
k1 = 4
sig1 = modulate(bits1, k1, T1)
t1 = np.linspace(0, T1, fs)
assert len(sig1) == fs
assert np.allclose(sig1, np.cos(2 * np.pi * t1))

bits2 = '101010'
T2 = .25
k2 = 2
sig2 = modulate(bits2, k2, T2)
t2 = np.linspace(0, T2, int(fs * T2))
chunks2 = [np.cos(2 * np.pi * t2 * 3 / T2) for _ in range(3)]
assert len(sig2) == len(chunks2) * len(chunks2[0])
assert np.allclose(sig2, np.concatenate(chunks2))

bits3 = '000001010011100101110111'
T3 = .001
k3 = 3
sig3 = modulate(bits3, k3, T3)
t3 = np.linspace(0, T3, int(fs * T3))
chunks3 = [np.cos(2 * np.pi * t3 * v / T3) for v in range(1, 9)]
assert len(sig3) == len(chunks3) * len(chunks3[0])
assert np.allclose(sig3, np.concatenate(chunks3))

Demodulation

This orthgonality is what makes the scheme so nice from a decoding perspective: if you just take the inner product (here, this just corresponds to the dot product since we'll have to implement it with discretized data rather than the truly continuous waveforms that would be transmitted in practice) of one chunk of the received signal with each signal from the basis, you'll get zero if the signals differ, or $T / 2$ if they match! (This is, of course, assuming no noise, which never occurs in practice. But even with noise, the inner products will be close to these values, so as long as a large enough value of $T$ is chosen, and there's not too much noise, we can still decode the received signal back to binary.)

This suggests using an array of $N = 2^k$ "correlators" which take the dot product of the received signal (one chunk at a time) with each signal in the basis, and identifying which signal it was based on whichever inner product returns the highest value, and finally converting it back to binary to complete the demodulation process. Here is a block diagram representation:

The correlators return us the coefficients $Y_i,\ i \in \{1, 2, ..., N=2^k\}$, and we decode/demodulate the received signal based on whichever coefficent is the largest.

In a slightly more sophisticated version of Orthogonal Signaling, if all coefficients were below a certain threshold, we would want to declare an erasure and not decode at all rather than just choose the max, as this means that the transmission was corrupted so badly that it's possible that the correlators will get it wrong regardless. We will not implement this feature here, however.

Q1b: Orthogonal Signaling - Demodulation

Implement demodulate, which takes in a signal (the received waveform), the number of bits per chunk k, and the length T in seconds of each signal (which is also the fundamental period), and returns the original bitstream. Note that in real life, the sender/receiver would actually share knowledge of k and T since these are fixed parameters of our scheme that both parties must agree on, so it makes sense for both functions to take them in as arguments.

You may use the function int_to_bitstring if you find it useful; you are not required to use it. It takes in an integer n and a value, k, and returns the binary value of n (as a bitstring), padded to k bits.

def int_to_bitstring(n, k):
    n_binary = bin(n)
    n_binary = n_binary[2:]  # bin prepends a 0b, which we don't want
    return n_binary.zfill(k) # prepend with zeros to k bits

def demodulate(signal, k, T):
    """
    Returns the demodulated version of signal using an Orthogonal Signaling scheme
    with k bits per chunk and T second signals. The return value will be the recovered
    bitstring.
    
    Parameters:
    signal - The received signal.
    k      - The number of bits per chunk. You may assume len(bits) % k == 0, 
             if bits is the original message.
    T      - The duration, in seconds, of each individual modulated chunk.
    """
    ## TODO ##

To test your implementation of demodulate, run the block below which compares the output of your function for various inputs. If the cell successfully runs without any AssertionErrors, you are passing all the tests. These tests are not comprehensive, and are only meant as basic sanity checks.

# These tests should pass even if modulate doesn't work
t1 = np.linspace(0, 1, 44100)
sig1 = np.cos(2 * np.pi * t1)
assert '00' == demodulate(sig1, 2, 1)

t2 = np.linspace(0, .25, int(.25 * 44100))
sig2 = np.concatenate((np.cos(2 * np.pi * t2 * 8 / .25), np.cos(2 * np.pi * t2 / .25)))
assert '111000' == demodulate(sig2, 3, .25)

# These tests assume a working implementation of modulate
assert '1010' == demodulate(modulate('1010', 2, 1), 2, 1)
assert '0101101010011111' == demodulate(modulate('0101101010011111', 4, .1), 4, .1)
assert '0000000000' == demodulate(modulate('0000000000', 1, .001), 1, .001)
assert '1111111111' == demodulate(modulate('1111111111', 1, .001), 1, .001)

Q2: What's all that noise?

One of the most important concepts in statistics, signal processing, machine learning and pretty much any field involving analysis of data is noise. Broadly speaking, noise is anything that corrupts our signal.

Examples of Noise

For example, suppose we're recording someone singing, and someone yells from a nearby room and messes up part of the recording. This kind of noise is called "shot noise", because of its short duration and high intensity, like a gunshot. Thermal noise from nearby devices giving off heat is another type, as electric circuits have temperature-dependent behavior. Yet another is ADC jitter: in EE 120 we consider the model of uniform sampling, where we obtain a discrete-time signal from a continuous-time one by taking samples of it that are uniformly spaced in time. In practice, there is no such thing as uniform sampling - ADCs (Analog to Digital Converters), the devices that actually perform sampling, have some error margin in where they actually sample, known as "ADC jitter". Fortunately, hardware engineers have shrunk this error margin down to several orders of magnitude less than the actual sampling frequency in modern ADCs to the point where it's negligible from our point of view. You can check out some actual specs for a real ADC and its jitter (sometimes also referred to as "aperture uncertainty") in the first reference.

There are many other types of noise, but hopefully you can already see some of the reasons why careful engineering is necessary to design robust systems from both a hardware development point of view touched on in other classes as well as on the modeling side of things we consider in EE 120.

Simulating Noise

We often treat things as ideal in EE 120, abstracting all of these issues away, partially because it's impossible to properly introduce and understand the noise models used in practical signal processing applications without probability theory. That said, we do want to test how robust our scheme is against some form of noise, so we'll introduce some of the basic concepts here.

Additive White Gaussian Noise (AWGN)

Perhaps the most ubiquitous model for noise in all of engineering is Additive White Gaussian Noise, or AWGN. There's three elements to dissect about this type of noise:

1. ADDITIVE: Additive noise is simply noise that is modeled as corruption caused by adding something to our true signal. For example, if the signal sent is $x(t)$, and the noise is $n(t)$, then we receive $y(t) = x(t) + n(t)$. It may seem like all noise should just be modeled as additive - in many cases it is, but there are cases where addition is not the best model for how noise is introduced. For example, if you were modeling the noise from a signal bouncing off of walls, you'd likely model it as a multiplicative attenuation factor rather than the addition of some noise signal.
2. WHITE: White noise is a signal with a uniform power spectral density. The power spectral density of a signal is the Fourier Transform of its autocorrelation. So, when we say that white noise has a uniform power spectral density, what we mean is that if you took a white noise signal, computed its autocorrelation, and took the Fourier Transform of that, it would be a constant. The name comes from the fact that white light was originally (incorrectly) believed to have uniform intensity over the visible band of frequencies. The useful takeaway from this property is that the noise doesn't favor any range of frequencies over another - they're all corrupted the same amount, on average.
3. GAUSSIAN: It turns out that when you sum up a large number of different random quantities, they approach what's known as a Gaussian distribution (meaning that the actual value itself is random, but the probability it lies within some range is given by a distribution function, which here is the Gaussian, aka the bell curve you've likely seen before). This is known as the Central Limit Theorem (you're not responsible for knowing this). Since any given system has a large number of random noise sources contributing to the overall noise, which is ultimately what we care about from a modeling perspective, this model is a good one that's very common in signal processing and other related fields. The Gaussian is parameterized by two quantities: $\mu$, the average value of the noise, and $\sigma$, the standard deviation of the noise. In almost all cases, $\mu$ is taken to be zero (known as "zero-mean noise") - if it wasn't zero, there would be some systematic bias which could simply be removed post-processing by subtracting this mean. $\sigma$ controls how strong the noise is - a higher value of $\sigma$ means the noise is more powerful, on average. Think of $\sigma$ as a knob that you turn, and as you turn it, the noise gets stronger.

The popularity of this noise model has a lot to do with the Central Limit Theorem mentioned above, and it's what we'll use here. Fortunately, it's very easy to implement for simulations: numpy has a built in function for drawing samples from this distribution, np.random.normal, which takes in $\mu$ and $\sigma$ as its first and second parameters. Run the code below to see what the noise looks like for a few different values of $\sigma$; you don't need to write any code for this cell.

You're not responsible for any probability theory content in this class, so don't worry if you didn't fully understand the above explanation. We just want to expose you to the AWGN model and give you a chance to explore it here due to how ubiquitous it is in signal processing and many other areas of science and engineering as well.

noise_len = 500      # number of noise samples to get
mu = 0               # we assume zero mean noise
sigma = [.5, 5, 10]  # noise stddevs

plt.figure(figsize=(20, 3))
for i in range(3):
    noise = np.random.normal(mu, sigma[i], noise_len)
    plt.subplot(1, 3, i+1)
    plt.title("{0} Samples of AWGN, $\mu = 0, \sigma = {1}$".format(noise_len, sigma[i]))
    plt.ylabel("Noise Amplitude")
    plt.ylim([-40, 40])
    plt.plot(noise)
plt.show()

A few comments:

• Notice how the noise is centered around zero, as we would expect from setting $\mu = 0$.
• As expected, increasing $\sigma$ increases how widespread the noise amplitude is. In fact, probability theory tells us that for the Gaussian distribution, roughly 68% of our noise will have an amplitude within $\pm \sigma$ of our mean of zero, around 95% within $\pm 2\sigma$, and around 99.7% within $\pm 3\sigma$. Look at some of the values on your plots and see if this matches them!
• Additive noise can be negative, as you see here. You may be confused by this - can't noise only be added on top of our signal, not take away from it? Well, sort of - if you're thinking about literal extra sound that's picked up by, say, a receiver listening for a signal, than perhaps (keep in mind we use noise as a more general term in signal processing - it's not just auditory "noise", but anything that corrupts our signal). But remember, we're considering this noise as the contribution of many different sources, such as the fact that the signal will lose some energy as it travels through some wire with resistance or otherwise travel across some distance to the receiver, in which case the signal's actual amplitude will go down from attenuation, hence "negative" noise.

Signal-to-Noise Ratio

Now that we have an idea of what noise is, in particular AWGN, we can study one of the most commonly cited metrics for signal fidelity: Signal-to-Noise Ratio, or SNR for short (sometimes also abbreviated as S/N). SNR allows us to quantify how reliable a particular signal's information is. A higher SNR is always better, as the higher the SNR, the more powerful your signal is relative to the noise corrupting it.

Intuition suggests that since an SNR of 1 is where the "amount" of signal is equal to the "amount" of noise, any SNR above 1 should be good. In most applications, however, the noise will cause noticeable degradation to the signal fidelity for any SNR below around 1000 (this is a coarse estimate, not some exact end-all threshold). An SNR of zero means we have no signal, only noise; an SNR of infinity means we have no noise (or at least noise so weak it can't be measured).

If you take any signal processing classes beyond EE 120, or work on any signal processing projects, you're guaranteed to see this concept come up - you have to deal with noise in all real-world applications.

Formally, SNR has a very intuitively satisfying defintion: it is the ratio of signal energy to noise energy:

$$SNR = \frac{\sum_{n=0}^{N-1}|\text{Signal}[n]|^2}{\sum_{n=0}^{N-1}|\text{Noise}[n]|^2}$$

where the above formula assumes finite duration signals of length $N$ (in general, the energy is computed with a sum over all integers, but on a computer, things always have finite length).

So, SNR measures how many times stronger the true signal is than the noise, in terms of energy. Remember the interpretation of signal energy as the "L2 norm of a signal" from Lab 3? Based on that, SNR can be thought of as (the square of) "how many times longer" your signal is in $N$-dimensional Euclidean space than the noise is!

There's two important things to note about this definition:

1. SNR is a single number describing the noise level, rather than a function at each point in time. Since noise is random, it's only meaningful to talk about its properties in terms of large sample sizes, hence why we use noise energy.
2. SNR relies on knowledge of the ground truth signal as well as the noise, which in practice we don't technically know. Fortunately, in real life, we often have meta-information about the system we're using that we can use to get a good estimate of the SNR. For example, in audio engineering, sinusoids at certain frequencies are the relevant signals, which are often used as the "signal" part of the calculation, and the noise of a system is typically a property of phyiscal hardware that can be measured and/or will be supplied by the hardware vendor.

SNR can also be measured in decibels by taking 20 times the log base 10 of the above expression (just as you saw in EE 16B), although we'll keep it as just the energy ratio for now.

Based on this definition, implement the function SNR below to calculate the signal-to-noise ratio.

def SNR(signal, noise):
    """
    Returns the Signal-to-Noise Ratio (SNR) of signal assuming corruption by noise, where
    SNR is defined as the ratio of the energy of signal to the energy of noise.
    """
    ## TODO ##

Run the block below to test your SNR function and also see some example values.

The plots (going top left to top right to bottom left to bottom right) should have SNRs of roughly 20000, 200, 2, and .02 respectively. If you're close to these, you're fine - keep in mind the noise generated is random, and that the whole point of calculating SNR via an average is that the noise is random, and so we can eliminate any bias that outliers cause. As we increase the number of samples for the SNR calculation, it gets more and more accurate. We only use 10000 samples here, which is pretty good for the noise levels we use, but there will still be some variance run to run, which is okay.

t = np.linspace(0, 1, 10000)
signal = np.cos(2 * np.pi * 4 * t)
sigmas = [.005, .05, .5, 5]
noises = [np.random.normal(0, sigma, len(signal)) for sigma in sigmas]
plt.figure(figsize=(20, 10))
for i in range(4):
    plt.subplot(2, 2, i+1)
    plt.plot(t, signal + noises[i])
    plt.xlabel("Time (s)")
    plt.ylabel("Amplitude of Noised Signal")
    plt.title("4 Hz, Unit Amplitude Cosine with AWGN with $\sigma = {0}$; SNR = {1}" \
              .format(sigmas[i], round(SNR(signal, noises[i]), 3)))

Q: Comment on the plots: how does SNR change with $\sigma$ (as $\sigma$ increases, does the SNR increase or decrease)? Why?

Answer: (TODO)

Q3: Testing Our Scheme!

Let's actually simulate some communications now that we've implemented Orthogonal Signaling and have an understanding of noise! We'll use the AWGN model from problem 2, assuming zero-mean noise, as is often done in both theory and practice due to Central Limit Theorem we discussed in problem 2.

One of the most commonly used metrics for how robust a digital communication system is Bit Error Rate, or BER, which is defined as the the number of incorrect bits (on the receiver end) over the number of total bits. Note that $0 \leq BER \leq 1$. A BER of 0 corresponds to no errors, and a BER of 1 corresponds to all bits being incorrect. In practice, design specifications for digital communication systems routinely require BERs of less than $10^{-12}$.

In this question, we'll simulate actual transmissions using our AWGN model, and gradually crank up the noise level by tuning $\sigma$ (which in turn will decrease our SNR), and see how BER changes as a function of SNR.

Q3a: Bit Error Rate

Implement the function BER below to calculate the Bit Error Rate of a received bitstring against a true bitstring.

def BER(received, actual):
    if len(received) != len(actual):
        raise ValueError("Inputs must have same length.")
    ### TODO ###

Test out your BER function by running this cell as a sanity check. If you don't get any AssertionErrors, then your BER function passes the tests.

def random_bitstring(L):
    """
    Returns a random bitstring of length L, in the sense that each bit 
    is equally likely to be a 1 or 0.
    """
    return int_to_bitstring(np.random.randint(0, 2**L - 1), L)

# use np.isclose rather than == to compare floating point calculation results
bits = '0000'
assert np.isclose(BER(bits, bits), 0)
assert np.isclose(BER('0001', bits), .25)
assert np.isclose(BER('0101', bits), .5)
assert np.isclose(BER('1111', bits), 1)

bits = '1010101'
assert np.isclose(BER(bits, bits), 0)
assert np.isclose(BER('0000000', bits), 4 / 7)
assert np.isclose(BER('1111111', bits), 3 / 7)
assert np.isclose(BER(bits[1:] + '0', bits), 1)

bits       = '010101010101010111100010101010101'
rcvd_bits1 = '010101001001010001010001001001101' 
rcvd_bits2 = '111111111111111111111111111111111'
assert np.isclose(BER(bits, bits), 0)
assert np.isclose(BER(rcvd_bits1, bits), 12 / len(bits))
assert np.isclose(BER(rcvd_bits2, bits), 16 / len(bits))

# Generate 1000 random 60-length bitstrings, and make sure we never go outside
# the range of [0, 1] for the BER - we don't actually care if the exact value is 
# correct here; in fact it most likely won't ever be
L = 60
N = 1000
bits = random_bitstring(L)
rcvd_bitstrs = [random_bitstring(L) for _ in range(N)]
for bitstr in rcvd_bitstrs:
    assert 0 <= BER(bitstr, bits) <= 1 # yes, you can do this in Python 3.. hehe

Q3b: Simulating Communication

Now, we'll simulate communication to see our Orthogonal Signaling in action! We'll do this by taking the bits we want to send, modulating, adding AWGN (to simulate the actual transmission), and demoduating to recover the sent bitstream.

Implement the send_and_receive function to simulate end-to-end communication of a bitstream. To generate the noise, use np.random.normal with a mean of zero and standard deviation of sigma. This should be fairly simple, given your work from question 1.

def send_and_receive(bits, k, T, sigma):
    """
    Simulates transmission of bits over an AWGN channel using a k bit, T duration
    Orthogonal Signaling scheme.
    
    Parameters:
    bits  - The bitstring to send.
    k     - The number of bits per chunk for Orthogonal Signaling.
    T     - Duration of each modulated chunk, in seconds.
    sigma - Standard deviation of the zero-mean AWGN.
    
    Returns:
    rcvd     - The received bitstring.
    sig_sent - The sent analog signal.
    sig_rcvd - The received analog signal.
    snr      - The SNR of the received signal.
    ber      - The BER of the received bitstring versus the transmitted one.
    """
    ### TODO ###
    return rcvd, sig_sent, sig_rcvd, snr, ber

Now, run the cell below to simulate some transmissions for different noise levels, and answer the questions below the plots. We'll plot the received signals, and print the corresponding SNR and BER for each transmission to see just how robust Orthogonal Signaling is!

We'll use k=2 and T=.001 here. We've typically been using T=1 so far, but this corresponds to taking a whole second to send 2 bits, which is unrealistic for practical transmissions. This set of parameters corresponds to 2 bits every millisecond, which is still only 2000 bits per second (try comparing that to the download speed of your home Internet download speed!) but will do for now.

Note that re-running the cell will produce slightly different SNR and BER numbers due to the random nature of the simulation we're performing; ballpark numbers are good enough for our purposes here.

## TODO run me

sigmas = [.01, .1, 1, 4, 10] # noise levels
bits = random_bitstring(60) # generate and use same bitstring across all noise levels
plt.figure(figsize=(18, 25))
for i in range(len(sigmas)):
    # Average the SNR/BER numbers over 100 trials to reduce outlier noise and get reasonably
    # consistent estimates
    k = 2
    T = .001
    num_trials = 100
    bits_rcvd, sig_sent, sig_rcvd, avg_snr, avg_ber = send_and_receive(bits, k, T, sigmas[i])
    for _ in range(num_trials - 1):
        _, _, _, snr, ber = send_and_receive(bits, k, T, sigmas[i])
        avg_snr += snr
        avg_ber += ber
    avg_snr = avg_snr / num_trials
    avg_ber = avg_ber / num_trials
    
    # Plot received signal (blue) and sent signal (bright green) over each other, and include
    # the average SNR, BER for the given noise level
    t = np.linspace(0, len(signal) / fs, len(sig_sent)) * 1000 # x-axis in milliseconds
    plt.subplot(len(sigmas), 1, i+1)
    plt.plot(t, sig_rcvd)
    plt.plot(t, sig_sent, '#54ff44', linewidth=2)
    plt.xlabel("Time (ms)")
    plt.ylabel("Signal Amplitude")
    plt.title("Received Signal; SNR: {0}, BER of Decoded Result: {1}".format(round(avg_snr, 5), round(avg_ber, 2)))
    plt.legend(('Received Signal', 'Sent Signal'))
plt.show()

Looking at the plots, we can see that Orthogonal Signaling is very robust to AWGN - even at low SNR, we can achieve perfect recovery of the sent bitstream for our given choices of k and T!

Eventually, as we increase SNR higher and higher, errors will start to come, and even the best algorithms and digital communication techniques can only handle so much noise.

Q3c: Parameter Tuning

Now that we've fully implemented Orthogonal Signaling, we can actually consider how varying our choices of the two parameters k (chunk size) and T (signal duration) leads to different tradeoffs.

Answer the following questions, assuming a fixed message length (in bits).

Q: If we increase k, keeping T fixed, does the message length (in terms of actual number of seconds to transmit) increase, decrease, or stay the same? Explain why.

A: (TODO)

Q: If we increase T, keeping k fixed, does the message length (in terms of actual number of seconds to transmit) increase, decrease, or stay the same? Explain why.

A: (TODO)

Q4: Play That Funky Music

Now, you're going to take the transmitter and reciever that you built in the earlier parts of the lab and use them to "send" some audio! In the lab, we have included a short violin excerpt of Partita E major, Gavotte en rondeau by Bach, and we will see how Orthogonal signaling lets us "transmit" and "recieve" this audio clip with some noise.

You do not have to write any code for this question - just run all the cells to reap the benefits of your work in previous parts!

The last two cells will take a while to run, about 2-5 minutes each. Don't worry too much if it seems like it is taking a while. It can take around 15 seconds to generate the bitstring and between 3-4 minutes to transmit and receive the data. After you implement this part, we will discuss why it takes so long, and what can be done to speed it up, which we won't do in this lab.

from bitstring import BitArray
def gen_bitstring(data):
    '''
    input: data, an np array of 16 bit integers that represents the music data
    output: bitstring, a large bitstring that represents our data
    '''
    return sum([BitArray(int=i, length=16) for i in data]).bin
    
def decode_bitstring(bits):
    '''
    input: bitstring (data type str), a large bitstring that represents our data
    output: data, an np array of 16 bit integers that represents the music data
    '''
    bit_arr = BitArray(bin=bits)
    decode_dat = []
    for i in range(len(bit_arr)//16):
        curr_int = bit_arr[16*i:16*(i+1)]
        decode_dat.append(curr_int.int)
    return np.array(decode_dat)

# Small sanity check that your code above works
arr = np.random.randint(-2**15, 2**15-1, 60, 'int16')
assert (arr == decode_bitstring(gen_bitstring(arr))).all()

from IPython.display import Audio
from scipy.io import wavfile
import time

fs, data = wavfile.read('a2002011001-e02.wav')
data = data[:,0] # Selecting the left channel from stereo data
n_sec = 3.5
data = data[:int(fs*n_sec)]

Audio(data=data, rate=fs)

# Plot audio signal
plt.figure(figsize=(15,10))
plt.title("Audio Signal")
plt.xlabel("Time (sec)")
plt.ylabel("Amplitude")
plt.plot(np.linspace(0, len(data) / fs, len(data)), data)
plt.show()

# Generate bitstring from data
t0 = time.time()
bitstring = gen_bitstring(data)
print("Time taken: {}".format(time.time()-t0))

# Low error, perfect reconstruction
k = 4
T = .001
sigma = 0.01 # Small!

t0 = time.time()
bits_rcvd, sig_sent, sig_rcvd, avg_snr, avg_ber = send_and_receive(bitstring, k, T, sigma)
print("Time taken: {}".format(time.time()-t0))

rcvd_data = decode_bitstring(bits_rcvd)
Audio(data=rcvd_data, rate=fs)

print("Average SNR: {}".format(avg_snr))
print("Average BER: {}".format(avg_ber))

# Higher error, nonideal reconstruction
k = 4
T = .001
sigma = 1 # Big!

t0 = time.time()
bits_rcvd, sig_sent, sig_rcvd, avg_snr, avg_ber = send_and_receive(bitstring, k, T, sigma)
print("Time taken: {}".format(time.time()-t0))

rcvd_data = decode_bitstring(bits_rcvd)
Audio(data=rcvd_data, rate=fs)

print("Average SNR: {}".format(avg_snr))
print("Average BER: {}".format(avg_ber))

For the second recovery, you should hear a "scratching" sound - this is the noise!

Improving Runtime

As we mentioned above (and as you saw while working through the question), encoding and decoding the data took a really long amount of time.

There are several factors contributing to this, but the main one is the sheer amount of data. Here, even wtih just 3.5 seconds of music, we have around 2.5 million bits of data to transmit. This is due to the fact that a .wav file is an uncompressed audio file format that uses 16 bits to represent each value.

In practice, one would never transmit uncompressed data. There are many different ways to reduce the amount of data to send, including (but not limited to):

• Lossless compression, such as Huffman (used in MP3) or Arithmetic Codes.
• Lossy Compression, such as Quantization and Thresholding (used by JPEG and many image compression algorithms)
• Transmitting some transform of the data, where the specific transform chosen is based on the desire for the signal to be sparse in the representation that transform uses:
  • Fourier Transform (great for combinations of pure tones, similar to what we have here)
  • Discrete Cosine Transform (used in images)
  • Wavelet Transform (many variants; most popular is the Haar wavelet basis, which is used for signals that are roughly piecewise-constant as the basis consists of shifted box functions)

Often, these methods are combined, and a lot of file formats you are familiar with (mp3, jpeg) are doing exactly that. If you are interested in learning more about compression, you can read this link, and can take further classes such as EE 123 (digital signal processing), EE 229A (information theory), or EE 225B (image processing) to learn more.

The decoding itself also adds time, but with a small $k$, we see that it is the data size that contributes the most to the runtime.

Improving Noise Resiliency

Another thing that could be improved about our scheme is how noise resilient it is. As discussed in the previous section, we modulate the raw data. This does create potential issues with the size of the data, but also with how robust our communication method is. Orthogonal signaling will do its best, but after demodulating, it'd be best if we had another layer of defense, so to speak, that could safeguard against potential errors in the bits themselves. This is where error-correcting codes come into play. We won't discuss them here, but in practice, many different codes such as Reed-Solomon codes, Hamming codes, convolutional codes, and more are examples of extra noise resiliency mechanisms that can be added onto the front and end of our communication system to help correct errors by reducing additional redundancy bits.

Concluding Remarks

As we saw here, Orthogonal Signaling is very resilient to noise due to the orthogonality of the different signals/messages sent, because they're very "far apart" - the noise has to work very hard to "push them together", and as long as we pick a large enough value of T for the SNR we're dealing with, we'll be able to win this metaphorical game of tug-of-war.

In practice, Orthogonal Signaling is seldom used due to bandwidth constraints: to send messages as fast as possible but with a low BER, we need a high T and high k. As we saw in question 1, the number of different signals (and thus the number of different frequencies we send at) grows exponentially with k. In practice, we only have a limited bandwidth that we can actually transmit over, due to both FCC allocation constraints as well as actual limits on our hardware. To borrow the example given in the EE 121 notes, if we used an Orthogonal Signaling scheme with $T = 1$ second and wanted to send 30 bits per second (so $k = 30$), then we would need $2^{30}$ different frequencies, which is a bandwidth of just over 1 GHz, which is something we'll never have in practice. Check out the fifth reference if you're interested in seeing how the actual US radio frequencies are allocated by the FCC.

That said, analyzing and understanding why its orthogonality property makes it so robust aids the design of more practical schemes. The tradeoff off having this strong property of orthogonality is that we have to use up a huge amount of bandwidth - to solve that problem, you can probably guess that we just need something that's similar to orthogonality, but isn't as strong - linear independence! Check out the fourth reference if you're interested in reading more.

References

[1] https://www.analog.com/media/en/technical-documentation/data-sheets/ad9461.pdf - example ADC data sheet, including consideration of various noise sources such as jitter.
[2] http://www-inst.eecs.berkeley.edu/~ee121/fa13/notes/lec05.pdf - Orthogonal Signaling notes from EE 121, UC Berkeley's old Digital Comm class.
[3] https://inst.eecs.berkeley.edu/~ee121/fa13/notes/lec04.pdf - Noise modeling notes from EE 121.
[4] http://www-inst.eecs.berkeley.edu/~ee121/fa13/notes/lec06.pdf - Followup to Orthogonal Signaling where orthogonality is traded for independence to help solve the bandwidth issue.
[5] https://www.ntia.doc.gov/files/ntia/publications/2003-allochrt.pdf - FCC spectrum allocation chart.
[6] http://www.music.helsinki.fi/tmt/opetus/uusmedia/esim/index-e.html - Data for Q4