Merge pull request #3 from AUAS-Pulsar/vzero

modularity for algorithms

Author: maxvhas

examples/visualize_psd.py

@@ -0,0 +1,42 @@
+import os,sys,inspect
+currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
+parentdir = os.path.dirname(currentdir)
+sys.path.insert(0,parentdir) 
+
+from rtlsdr import RtlSdr
+import matplotlib.pyplot as plt
+import numpy as np
+from fourier import fourier
+from plot import psd
+
+# Initiate RtlSdr
+sdr = RtlSdr()
+
+# configure device
+sdr.sample_rate = 2.4e6
+sdr.center_freq = 102.2e6
+
+# Read samples
+samples = sdr.read_samples(1024)
+
+# Close RTLSDR device connection
+sdr.close()
+
+# Number of samples equals the length of samples
+N = samples.shape[0]
+
+# T equals N/Fs 
+T = N/sdr.sample_rate
+
+# Get the powerlevels and the frequencies
+Pxx, freqs, _ = psd(samples, NFFT=1024, Fs=sdr.sample_rate/1e6, scale_by_freq=True, sides='twosided')
+
+# Calculate the powerlevel dB's 
+power_levels = 10*np.log10(Pxx/(sdr.sample_rate/1e6))
+
+# Add the center frequency to the frequencies so it matches the actual frequencies
+freqs = freqs + sdr.center_freq/1e6
+
+# Plot the PSD
+plt.plot(freqs, power_levels)
+plt.show()

examples/visualize_time_domain.py

@@ -0,0 +1,31 @@
+from rtlsdr import RtlSdr
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Initiate RtlSdr
+sdr = RtlSdr()
+
+# configure device
+sdr.sample_rate = 2.4e6
+sdr.center_freq = 102.2e6
+# sdr.gain = 0
+
+# Read samples
+samples = sdr.read_samples(256*1024)
+
+# Close RTLSDR device connection
+sdr.close()
+
+# Number of samples equals the length of samples
+N = samples.shape[0]
+
+# T equals N/Fs 
+T = N/sdr.sample_rate
+
+# Define the x axis
+x = np.linspace(0.0, T, N)
+
+# Generate the plot
+fig, ax = plt.subplots()
+ax.plot(x, samples)
+plt.show()

fourier/__init__.py

@@ -0,0 +1 @@
+from .fourier import *

fourier/fourier.py

@@ -0,0 +1,88 @@
+import numpy as np    
+
+
+def DFT_slow(x):
+    """ Compute the discrete Fourier Transform of the one-dimensional array x"""
+    x = np.asarray(x)
+    N = x.shape[0]
+    n = np.arange(N)
+    k = n.reshape((N, 1))
+    M = np.exp(-2j * np.pi * k * n / N)
+    return np.dot(M, x)
+
+
+def FFT_vectorized(x):
+    """A vectorized, non-recursive version of the Cooley-Tukey FFT"""
+    x = np.asarray(x)
+    N = x.shape[0]
+
+    if np.log2(N) % 1 > 0:
+        raise ValueError("size of x must be a power of 2")
+
+    # N_min here is equivalent to the stopping condition above,
+    # and should be a power of 2
+    N_min = min(N, 32)
+    
+    # Perform an O[N^2] DFT on all length-N_min sub-problems at once
+    n = np.arange(N_min)
+    k = n[:, None]
+    M = np.exp(-2j * np.pi * n * k / N_min)
+    X = np.dot(M, x.reshape((N_min, -1)))
+
+    # build-up each level of the recursive calculation all at once
+    while X.shape[0] < N:
+        X_even = X[:, :X.shape[1] // 2]
+        X_odd = X[:, X.shape[1] // 2:]
+        factor = np.exp(-1j * np.pi * np.arange(X.shape[0])
+                        / X.shape[0])[:, None]
+        X = np.vstack([X_even + factor * X_odd,
+                       X_even - factor * X_odd])
+
+    return X.ravel()
+
+
+def fftfreq(n, d=1.0):
+    """
+    Return the Discrete Fourier Transform sample frequencies.
+
+    The returned float array `f` contains the frequency bin centers in cycles
+    per unit of the sample spacing (with zero at the start).  For instance, if
+    the sample spacing is in seconds, then the frequency unit is cycles/second.
+
+    Given a window length `n` and a sample spacing `d`::
+
+      f = [0, 1, ...,   n/2-1,     -n/2, ..., -1] / (d*n)   if n is even
+      f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (d*n)   if n is odd
+
+    Parameters
+    ----------
+    n : int
+        Window length.
+    d : scalar, optional
+        Sample spacing (inverse of the sampling rate). Defaults to 1.
+
+    Returns
+    -------
+    f : ndarray
+        Array of length `n` containing the sample frequencies.
+
+    Examples
+    --------
+    >>> signal = np.array([-2, 8, 6, 4, 1, 0, 3, 5], dtype=float)
+    >>> fourier = np.fft.fft(signal)
+    >>> n = signal.size
+    >>> timestep = 0.1
+    >>> freq = np.fft.fftfreq(n, d=timestep)
+    >>> freq
+    array([ 0.  ,  1.25,  2.5 ,  3.75, -5.  , -3.75, -2.5 , -1.25])
+
+    """
+    
+    val = 1.0 / (n * d)
+    results = np.empty(n, int)
+    N = (n-1)//2 + 1
+    p1 = np.arange(0, N, dtype=int)
+    results[:N] = p1
+    p2 = np.arange(-(n//2), 0, dtype=int)
+    results[N:] = p2
+    return results * val

plot/__init__.py

@@ -0,0 +1 @@
+from .plot import *