End sem 7; week 2 winter session 2025

7th-Semester-Fall-2024/ECOMMS/final-project/bfm_demod.py

@@ -0,0 +1,80 @@
+import numpy as np
+import scipy as sp
+import sounddevice as sd
+from threading import Thread
+import time
+import queue
+
+import matplotlib.animation as anim
+import matplotlib.pyplot as plt
+
+from rtlsdr import RtlSdr
+import asyncio
+import dsp
+
+class BfmDemod: 
+    def __init__(self, sdr:RtlSdr, stream:sd.OutputStream, audio_sample_rate:int=48000, num_samples:int=8192):
+        self.sdr = sdr
+        self.stream = stream
+        self.audio_sample_rate = audio_sample_rate
+        self.buffer = queue.Queue(10)
+        self.num_samples = num_samples
+
+        self.fig, self.ax = plt.subplots()
+        self.duration = self.num_samples/self.sdr.sample_rate
+        #lower_extent = self.sdr.center_freq - self.sdr.sample_rate/2
+        #upper_extent = self.sdr.center_freq + self.sdr.sample_rate/2
+
+        #self.fig, self.ax = plt.subplots()
+        #self.ax.xlabel = "Frequency [Hz]"
+        #self.ax.ylabel = "Amplitude [dB]"
+        #self.ax.set_xlim(left=lower_extent, right=upper_extent)
+    
+    # Plot the spectrum of the incoming samples
+    def _spectrum_update(self, frame):
+        self.ax.clear()
+        if self.buffer is not None:
+            spectrum = np.abs(sp.fft.fftshift(sp.fft.fft(self.buffer)))**2
+            spectrum_db = dsp.db(spectrum)
+            self.ax.plot(spectrum_db)
+        return self.ax,
+
+    # Begin reading and demodulation
+    async def start(self):
+        self.stream.start()
+        plt.ion()
+        self.scatter = self.ax.scatter([], [])
+        plt.show()
+
+        async for samples in self.sdr.stream(self.num_samples):
+            message = dsp.quad_demod(samples)
+            message = dsp.lowpass(message, 16E3, self.sdr.sample_rate)
+            message /= np.max(np.abs(message))
+            
+            time = len(message)/self.sdr.sample_rate
+            num_samples = int(time*self.audio_sample_rate)
+            message = sp.signal.resample(message, num_samples)
+
+            message = message.astype(np.float32)
+            #self.stream.write(message)
+
+            spectrum = np.abs(sp.fft.fftshift(sp.fft.fft(samples)))**2
+            spectrum_db = dsp.db(spectrum)
+            self.ax.clear()
+            self.ax.plot(spectrum_db)
+
+            self.ax.relim()
+            self.ax.autoscale_view()
+
+            self.fig.canvas.draw()
+            self.fig.canvas.flush_events()
+
+        await sdr.stop()
+        
+
+    def stop(self):
+        self.sdr.cancel_read_async()
+        self.reader.join()
+        self.stream.stop()
+
+

7th-Semester-Fall-2024/ECOMMS/final-project/dsp.py

@@ -0,0 +1,37 @@
+import scipy as sp
+import numpy as np
+
+
+
+def chunk_samples(samples, num_rows):
+    return samples.reshape((samples.shape[0]//num_rows, num_rows))
+
+
+def db(x):
+    return 10*np.log10(x)
+
+
+def quad_demod(samples):
+    return 0.5*np.angle(samples[:-1] * np.conj(samples[1:]))
+
+
+def lowpass(data, f_cutoff, f_s):
+    nyq = 0.5*f_s
+    normal_cutoff = f_cutoff/nyq
+    b, a = sp.signal.butter(2, normal_cutoff, btype="low", analog=False)
+    y = sp.signal.lfilter(b, a, data)
+    return y
+
+
+def sdr_bandpass(data, f_pass, f_cutoff, extent_MHz):
+
+    lower_extent = extent_MHz[0]*1E6
+    upper_extent = extent_MHz[1]*1E6
+    extent_range = upper_extent - lower_extent
+
+    normal_pass = (f_pass - lower_extent)/extent_range
+    normal_cutoff = (f_cutoff - lower_extent)/extent_range
+
+    b, a = sp.signal.butter(2, (normal_pass, normal_cutoff), btype="bandpass", analog=False)
+    y = sp.signal.lfilter(b, a, data)
+    return y

7th-Semester-Fall-2024/ECOMMS/final-project/gfsk_demod.py

@@ -0,0 +1,58 @@
+import numpy as np
+import scipy as sp
+from threading import Thread
+import time
+import queue
+
+import matplotlib.animation as animation
+import matplotlib.pyplot as plt
+
+from rtlsdr import RtlSdr
+
+import dsp
+
+
+class GFSKDemod: 
+    def __init__(self, sdr:RtlSdr, baud_rate:int=9600, num_samples:int=512):
+        self.sdr = sdr
+        self.baud_rate = baud_rate
+        self.num_samples = num_samples
+        self.fig, self.ax = plt.subplots()
+
+        self.duration = self.num_samples/self.sdr.sample_rate
+
+    # Demodulate an update the sound buffer
+    def _demod(self, samples):
+        message = dsp.quad_demod(samples)
+        message = dsp.lowpass(message, 16E3, self.sdr.sample_rate)
+        #message /= np.max(np.abs(message))
+        #message = message.astype(np.float32)
+        time = len(message)*self.sdr.sample_rate
+        t = np.linspace(0,time,len(message))
+        return t, message
+
+
+    # Begin reading and demodulation
+    async def start(self):
+        plt.ion()
+        plt.show()
+
+        async for samples in self.sdr.stream(self.num_samples):
+            self.ax.clear()
+            t, message = self._demod(samples)
+
+            #num_samples = int(len(message)*self.resample_rate/self.sdr.sample_rate)
+            #message = sp.signal.resample(message, num_samples)
+
+            #t = np.linspace(0,1,len(message))*t[:-1]
+            #self.ax.plot(t,message)
+            #spectrum = np.abs(sp.fft.fftshift(sp.fft.fft(samples)))**2
+            #spectrum_db = dsp.db(spectrum)
+            #self.ax.plot(spectrum_db)
+
+            self.fig.canvas.draw()
+            self.fig.canvas.flush_events()
+
+
+    def stop(self):
+        self.sdr.cancel_read_async()

7th-Semester-Fall-2024/ECOMMS/final-project/sd_test.py

@@ -0,0 +1,20 @@
+import numpy as np
+import sounddevice as sd
+import scipy as sp
+from scipy.io import wavfile
+import matplotlib.pyplot as plt
+
+sample_rate, wav_data = wavfile.read("fm.wav")
+time = len(wav_data)/sample_rate
+
+new_sample_rate = 48E3
+num_samples = int(time*new_sample_rate)
+plt.subplot(211)
+plt.plot(wav_data)
+wav_data = sp.signal.resample(wav_data, num_samples)
+plt.subplot(212)
+plt.plot(wav_data)
+plt.show()
+
+sd.play(wav_data, samplerate=new_sample_rate)
+sd.wait()

7th-Semester-Fall-2024/ECOMMS/final-project/sdr_stream.py

@@ -1,13 +1,31 @@
-from rtlsdr import RtlSdr
-import numpy as np
-import matplotlib.pyplot as plt
-import scipy as sp
+#from rtlsdr.rtlsdrtcp.client import RtlSdrTcpClient
 import asyncio
+from rtlsdr import RtlSdr
+from bfm_demod import BfmDemod
+from gfsk_demod import GFSKDemod
+import sounddevice as sd
+import numpy as np
 
-async def streaming():
+
+def main():
+    #sdr = RtlSdrTcpClient(hostname='localhost', port=12345)
     sdr = RtlSdr()
-    
-    async for samples in sdr.stream():
 
-    await sdr.stop()
+    sdr.set_sample_rate(2.048E6)
+    sdr.set_center_freq(434.55E6)
+    sdr.set_freq_correction(60)
+    sdr.set_gain("auto")
+
+    stream = sd.OutputStream(samplerate=48E3, dtype=np.float32, channels=1)
+
+    sink = GFSKDemod(sdr=sdr)
+    #sink = BfmDemod(sdr=sdr, stream=stream)
+    
+    loop = asyncio.get_event_loop()
+    loop.run_until_complete(sink.start())
+
     sdr.close()
+
+
+if __name__ == "__main__":
+    main()

7th-Semester-Fall-2024/ECOMMS/final-project/spectrum.py

@@ -0,0 +1,27 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy as sp
+from rtlsdr import RtlSdr
+import asyncio
+import dsp
+
+class Spectrum(Display):
+    def __init__(self):
+        pass
+
+    def put(self, samples):
+        self.buffer.put(samples)
+
+    def update(self):
+        samples = self.buffer.get()
+        spectrum = np.abs(sp.fft.fftshift(sp.fft.fft(samples)))**2
+        spectrum_db = dsp.db(spectrum)
+        self.ax.clear()
+        self.ax.plot(spectrum_db)
+
+        self.ax.relim()
+        self.ax.autoscale_view()
+
+        self.fig.canvas.draw()
+        self.fig.canvas.flush_events()
+

7th-Semester-Fall-2024/ECOMMS/homework/homework-4/homework-4.md

@@ -0,0 +1,28 @@
+# Homework 4 - Aidan Sharpe
+
+## Problem 1
+A multilevel digital communication system sends one of 16 possible levels over the channel every 0.8[ms].
+
+### a
+What is the number of bits corresponding to each level? Or what is the number of bits per symbol?
+
+16 levels $= \log_2(16) = 4$ bits per symbol.
+
+### b
+What is the baud rate?
+
+1 symbol per 0.8[ms] $= \frac{1}{0.8 \times 10^{-3}} = 1250$ symbols per second.
+
+### c
+What is the bit rate?
+
+4 bits per symbol at 1250 baud is 5000[bps].
+
+## Problem 2
+The input to a differential coding signal is 10110010. Begin with reference digit 1. What is the encoded sequence?
+
+Apply XOR operation to each element, diagonally down and left, and store the result below.
+
+|   | 1 | 0 | 1 | 1 | 0 | 0 | 1 | 0 |
+|---|---|---|---|---|---|---|---|---|
+| 1 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 1 |

7th-Semester-Fall-2024/ECOMMS/homework/homework-5/homework-5.md

@@ -0,0 +1,31 @@
+# Homework 5 - Aidan Sharpe
+
+## Problem 1
+In a binary communication system, let the receiver test statistic of the received signal be $r_0$. The received signal consists of a polar digital signal plus noise. The polar signal has values $s_{01}=A$ and $s_{02}=-A$. Assume that the noise has a Laplacian distribution:
+$$f(n_0) = \frac{1}{\sqrt{2} \sigma_0} e^{-\sqrt{2}|n_0|/\sigma_0}$$
+where $\sigma_0$ is the RMS value of the noise, $f(n_0)$ is the probability density function (PDF), and $n_0$ is the signal. In the case of a PDF of $s_{01}$ and $s_{02}$, replace $n_0$ by $r_0-A$ and $r_0+A$. The shape of the PDF for $s_{01}$ and $s_{02}$ is the same. Find the probability of error $P_e$ as a function of $A/\sigma_0$ for the case of equally likely signaling and $V_T$ having the optimum value.
+
+
+Given the two PDFs corresponding to $s_1$ and $s_2$, the probability of a bit error is the same as the area of the intersection of the two PDFs, as seen in figure \ref{error_pdfs}.
+
+![PDFs corresponding to bit error probability \label{error_pdfs}](prob-err.png)
+
+The curve traced out by this area is
+$$p(r_0) = \frac{1}{2}\left[f(r_0 | s_1) + f(r_0 | s_2) - \|f(r_0 | s_1) - f(r_0 | s_2) \|\right]$$
+
+to find the area under the curve (the probability of a bit error $P_e$), we integrate $p(r_0)$ for all values of $r_0$. Since $f(r_0 | s_1)$ and $f(r_0 | s_2)$ are PDFs, the area under each must be unity. Therefore, distributing the $\frac{1}{2}$ term, we are left with:
+$$P_e = 1 - \frac{1}{2} \int\limits_{-\infty}^\infty  \left|f(r_0 | s_1) - f(r_0 | s_2)\right| dr_0$$
+
+Finally, we expand and simplify the integral to:
+
+$$P_e = 1 - \frac{\sqrt{2}}{4\sigma_0} \int\limits_{-\infty}^\infty  \left|e^{-\sqrt{2} |r_0-A| /\sigma_0} - e^{-\sqrt{2} |r_0 + A|/\sigma_0}\right| dr_0$$
+
+After evaluating the intergral, we find that $P_e = e^{-\sqrt{2}A/\sigma_0}$.
+
+
+## Problem 2
+A digital signal with white Gaussian noise is received by a receiver with a matched filter. The signal is a unipolar non-return to zero signal with $s_{01}=1$[V] and $s_{02}=0$[V]. The bit rate is 1 Mbps. The power spectral density of the noise is $N_0/2=10^{-8}$[W/Hz]. What is the probability of error $P_e$? Assume the white Gaussian noise is thermal noise. 
+
+For a unipolar signal received by a receiver with a matched filter, the probability of error is given by:
+$$P_e = Q\left(\sqrt{\frac{A^2 T}{4 N_0}}\right)$$
+where $A = 1 - 0 = 1$ is the amplitude and $T= 1[\mu$s]. Therefore, $P_e = 2.03 \times 10^{-4}$.

7th-Semester-Fall-2024/ECOMMS/homework/homework-5/problem-1.py

@@ -0,0 +1,30 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy as sp
+
+def laplacian_distribution(n_0, std_dev):
+    return np.exp(-np.sqrt(2)*np.abs(n_0)/std_dev) / (np.sqrt(2)*std_dev)
+
+def main():
+    r_0 = np.linspace(-5,5,1000)
+    A = 1
+    std_dev = 1
+
+    low_dist = laplacian_distribution(r_0 - A, std_dev)
+    up_dist = laplacian_distribution(r_0 + A, std_dev)
+
+    plt.plot(r_0, low_dist, label="$s_1$")
+    plt.plot(r_0, up_dist, label="$s_2$")
+
+    P_error = (up_dist + low_dist - np.abs(up_dist-low_dist))/2
+    plt.fill_between(r_0, P_error, color=(0,0,0,0), hatch="//", edgecolor='k', label="$P_e$")
+
+    plt.legend()
+
+    plt.savefig("Probability of error")
+
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()

7th-Semester-Fall-2024/ECOMMS/labs/lab2/am.py

@@ -2,7 +2,7 @@ import numpy as np
 import matplotlib.pyplot as plt
 import scipy as sp
 
-def add_noise_snr_db(s, SNR):
+def add_noise(s, SNR):
     var_s = np.cov(s)
     var_noise = var_s/(10**(SNR/10))
     noise = var_noise**0.5 * np.random.randn(len(s))
@@ -51,14 +51,6 @@ def main():
     s_am = dsb_am(m, f_c, t)
     s_sc = dsb_sc(m, f_c, t)
 
-    plt.subplot(311)
-    plt.plot(t, m, label="Message Signal")
-    plt.subplot(312)
-    plt.plot(t, s_am, label="DSB-AM Modulated Signal")
-    plt.subplot(313)
-    plt.plot(t, s_sc, label="DSB-SC Modulated Signal")
-    plt.show()
-
     m_am_demod = product_demod(s_am, f_m, f_c, f_s, t)
     m_sc_demod = product_demod(s_sc, f_m, f_c, f_s, t)
 
@@ -66,14 +58,17 @@ def main():
     plt.plot(t, m_am_demod, label="Demodulated DSB-AM Signal")
     plt.plot(t, m_sc_demod, label="Demodulated DSB-SC Signal")
     plt.legend(loc="upper right")
+    plt.savefig("demodulation-clean.png")
     plt.show()
 
     m_am_noisy_demod = product_demod(add_noise_snr_db(s_am, 3), f_m, f_c, f_s, t)
     m_sc_noisy_demod = product_demod(add_noise_snr_db(s_sc, 3), f_m, f_c, f_s, t)
 
-    plt.plot(t, m)
-    plt.plot(t, m_am_noisy_demod)
-    plt.plot(t, m_sc_noisy_demod)
+    plt.plot(t, m, label="Original Message Signal")
+    plt.plot(t, m_am_noisy_demod, label="Demodulated DSB-AM Signal With Noise")
+    plt.plot(t, m_sc_noisy_demod, label="Demodulated DSB-SC Signal With Noise")
+    plt.legend(loc="upper right")
+    plt.savefig("demodulation-noisy.png")
     plt.show()
 
 if __name__ == "__main__":

7th-Semester-Fall-2024/ECOMMS/labs/lab3/lab-3.py

@@ -0,0 +1,140 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy as sp
+
+def add_noise(s, SNR):
+    var_s = np.cov(s)
+    var_noise = var_s/(10**(SNR/10))
+    noise = var_noise**0.5 * np.random.randn(len(s))
+    return s + noise
+
+def add_noise_db(s, noise_db):
+    var_noise = 10**(noise_db/10)
+    noise = var_noise**0.5 * np.random.randn(len(s))
+    return s + noise
+
+
+def lowpass(data, f_cutoff, f_s):
+    nyq = 0.5*f_s
+    normal_cutoff = f_cutoff/nyq
+    b, a = sp.signal.butter(2, normal_cutoff, btype="low", analog=False)
+    y = sp.signal.lfilter(b, a, data)
+    return y
+
+def bandpass(data, f_pass, f_cutoff, f_s):
+    nyq = 0.5*f_s
+
+    normal_pass = f_pass/nyq
+    normal_cutoff = f_cutoff/nyq
+
+    b, a = sp.signal.butter(2, (normal_pass, normal_cutoff), btype="bandpass", analog=False)
+    y = sp.signal.lfilter(b, a, data)
+    return y
+
+def fm(m, f_c, beta_f, t):
+    omega_c = 2*np.pi*f_c
+    return np.cos(omega_c*t + beta_f*m)
+
+def product_demod(s, f_baseband, f_c, f_s, t):
+    omega_c = 2*np.pi*f_c
+    mixed = s*np.cos(omega_c*t)
+    return lowpass(mixed, f_baseband, f_s)
+
+
+def fm_demod(s_fm, f_c, f_s, beta_f, t):
+    diff_t = np.gradient(s_fm, t)
+    envelope = np.abs(sp.signal.hilbert(diff_t))
+
+    omega_c = 2*np.pi*f_c
+    dt = 1/f_s
+    m_demod = np.empty_like(envelope)
+
+    #err1 = dt*(envelope[0] + envelope[1] - 2*omega_c)/beta_f**2
+    #err2 = dt*(dm_dt[0] + dm_dt[1])/beta_f
+
+    for i in range(len(envelope)):
+        m_demod[i] = np.sum((envelope[:i+1] - omega_c)/(beta_f*f_s))
+
+    return m_demod
+
+def quad_demod(s_fm, bandwidth, f_s, f_c, t):
+    iq = real_to_iq(s_fm, f_c, bandwidth, f_s, t)
+    iq = np.conj(iq)
+    return 0.5*np.angle(iq)
+
+def real_to_iq(s_fm, f_c, bandwidth, f_s, t):
+    i = lowpass(s_fm * np.cos(2*np.pi*f_c*t), bandwidth, f_s)
+    q = lowpass(s_fm * np.sin(2*np.pi*f_c*t), bandwidth, f_s)
+    return i + 1j*q
+
+
+def m(t):
+    f_m = 10
+    omega_m = 2*np.pi*f_m
+    return 0.1*np.sin(omega_m*t)
+
+def dm_dt(t):
+    f_m = 10
+    omega_m = 2*np.pi*f_m
+    return 0.1*omega_m*np.cos(omega_m*t)
+
+def db(x):
+    return 10*np.log10(x)
+
+def main():
+    T_s = 0.0001
+    f_s = 1/T_s
+    f_m = 20
+    omega_m = 2*np.pi*f_m
+
+    f_c = 500
+
+    beta_f = 2.40483
+    
+    t = np.arange(0,1,T_s)
+    f = np.linspace(-f_s/2, f_s/2, len(t))
+
+    #m = 0.5*(np.sin(omega_m*t) + np.sin(omega_m/2*t))
+    m = np.sin(omega_m*t)
+    s_fm = fm(m, f_c, beta_f, t)
+
+    omega_c = 2*np.pi*f_c
+    
+    #plt.plot(t, m, label="Original Message Signal")
+    #plt.plot(t, quad_demod(s_fm, 2*f_m, f_s, f_c, t), label="Demodulated Message")
+    #plt.legend()
+    #plt.savefig("message-signals")
+    #plt.show()
+
+    #plt.plot(t, s_fm, label="Frequency Modulated Message")
+    #plt.legend()
+    #plt.savefig("fm-signal")
+    #plt.show()
+
+
+    spectrum = np.abs(sp.fft.fftshift(sp.fft.fft(s_fm)))**2
+    #spectrum_db = db(spectrum)
+
+    plt.subplot(121)
+    plt.plot(f, spectrum)
+    #plt.vlines(f_c, -250, 100, color='r')
+    plt.xlim(0,2*f_c)
+    plt.title("Null Carrier Spectrum")
+    plt.xlabel("Frequency [Hz]")
+
+    beta_f = 5
+    m = np.sin(omega_m*t)
+    s_fm = fm(m, f_c, beta_f, t)
+    spectrum = np.abs(sp.fft.fftshift(sp.fft.fft(s_fm)))**2
+
+    plt.subplot(122)
+    plt.title("Non-null Carrier Spectrum")
+    plt.plot(f, spectrum)
+    plt.xlim(0,2*f_c)
+    plt.xlabel("Frequency [Hz]")
+    plt.savefig("null-carrier-comp")
+
+    plt.show()
+
+if __name__ == "__main__":
+    main()