{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Lab2 - Time Frequency\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Import functions and libraries\n", "from __future__ import division\n", "import numpy as np, matplotlib.pyplot as plt\n", "from numpy import *\n", "from numpy.fft import *\n", "import scipy.signal as signal\n", "from matplotlib.pyplot import *\n", "from rtlsdr import RtlSdr\n", "import scipy\n", "import scipy.io\n", "from scipy.io.wavfile import write\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below are some functions from Part I that we will need in this part" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Plot an image of the spectrogram y, with the axis labeled with time tl,\n", "# and frequency fl\n", "#\n", "# t_range -- time axis label, nt samples\n", "# f_range -- frequency axis label, nf samples\n", "# y -- spectrogram, nf by nt array\n", "# dbf -- Dynamic range of the spectrum\n", "\n", "def sg_plot( t_range, f_range, y, dbf = 60, fig = None) :\n", " eps = 10.0**(-dbf/20.0) # minimum signal\n", " \n", " # find maximum\n", " y_max = abs(y).max()\n", " \n", " # compute 20*log magnitude, scaled to the max\n", " y_log = 20.0 * np.log10( (abs( y ) / y_max)*(1-eps) + eps )\n", " \n", " # rescale image intensity to 256\n", " img = 256*(y_log + dbf)/dbf - 1\n", " \n", " \n", " fig=figure(figsize=(16,6))\n", " \n", " plt.imshow( np.flipud( 64.0*(y_log + dbf)/dbf ), extent= t_range + f_range ,cmap=plt.cm.gray, aspect='auto')\n", " plt.xlabel('Time, s')\n", " plt.ylabel('Frequency, Hz')\n", " plt.tight_layout()\n", " \n", " return fig\n", " \n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import pyaudio\n", "import wave\n", "\n", "# function that reads wav file to array\n", "def read_wav( wavname ):\n", " \n", " wf = wave.open(wavname, 'rb')\n", " \n", " CHUNK = 1024\n", " frames = []\n", " data_str = wf.readframes(CHUNK) #read a chunk\n", " \n", " while data_str != '':\n", " data_int = np.fromstring( data_str, 'int16') # convert from string to int (assumes .wav in int16)\n", " data_flt = data_int.astype( np.float32 ) / 32767.0 # convert from int to float32\n", " frames.append( data_flt ) #append to list\n", " data_str = wf.readframes(CHUNK) #read a chunk\n", "\n", " return np.concatenate( frames )\n", "\n", "\n", "def play_audio( data, p, fs):\n", " # data - audio data array\n", " # p - pyAudio object\n", " # fs - sampling rate\n", " \n", " # open output stream\n", " ostream = p.open(format=pyaudio.paFloat32, channels=1, rate=fs,output=True)\n", " # play audio\n", " ostream.write( data.astype(np.float32).tostring() )\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part II: Broadcast FM Radio\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Task 5: Time-Frequency plots of the radio-frequency spectrum with the SDR." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The samples that are obtained by the SDR represent a bandwidth of the spectrum around a center frequency. Hence, when demodulating to base-band (i.e. zero frequency) the signal must be imaginary since it has a non symmetric Fourier transform. \n", "\n", "In this case, we would like to display both sides of the spectrum.\n", "\n", "- Modify the function `myspectrogram_hann_ovlp(x,m,fs,fc)` such that it detects whether the input signal `x` is complex. In that case, it will compute a double sided spectrum with frequency range centered around fc (center frequency). For this, it would be useful to use the commands: `isreal` and `fftshift`.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def myspectrogram_hann_ovlp(x, m, fs, fc, dbf = 60, fig = None):\n", " # Plot the spectrogram of x.\n", " # First take the original signal x and split it into blocks of length m\n", " \n", " # Your code here:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will first look at radio FM spectrum. In the US the broadcast FM radio band is 88.1-107.9Mhz. It is split into 200KHz slots. This is relatively a large bandwidth and therefore it is also called wideband FM as opposed to narrowband FM which can be as low as 5 Khz. In FM radio the information is encoded by modulating the frequency of the carrier, \n", "$$y_c(t) = A\\cos \\left (2\\pi f_c t + 2\\pi \\Delta f \\int_0^t x(\\tau) d\\tau \\right ).$$\n", "\n", "Here, $f_c$ is the carrier frequency, $\\Delta f$ is the frequency deviation and $x(t)$ is a normalized baseband signal, which contains all the information the station wants to broadcast. This includes stereo audio, as well as digital information on the station and sometimes additional narrow bandwidth subcarrier channels. \n", "\n", "\n", "\n", "More specifically, the broadcast FM baseband signal, $x(t)$, consists of mono (Left+Right) channels from 30Hz to 15 KHz, a pilot signal at $f_p=19$ KHz, amplitude modulated Stereo (Left - Right) channels around $2\\cdot f_p = 38$KHz and digital RBDS, which encodes the station information, song name etc. at $3\\cdot f_p =57$KHz. (See for more information). \n", "\n", "The baseband signal is:\n", "\n", "$ \\qquad \\qquad x(t) = \\underbrace{(L+R)}_{\\mathrm{mono}} + \\underbrace{0.1 \\cdot \\cos(2\\pi f_p t)}_\\mathrm{pilot} + \\underbrace{(L-R)\\cos(2\\pi (2f_p) t)}_\\mathrm{stereo} + \\underbrace{0.05\\cdot \\mathrm{RBDS}(t)\\cos(2\\pi (3f_p) t)}_\\mathrm{digital~ RBDS} + \\mathrm{subcarriers~signals}. $\n", "\n", "The signal $\\mathrm{RBDS}(t)$ is a $m(t)\\cos(2\\pi(3(f_p))$ where $m(t)$ is a binary signal consists of $\\pm1$ at constant intervals which encode 0, and 1 bits. The subcarriers are often FM modulated signals of foreign stations or carry other information.\n", "This is the spectrum of $x(t)$:\n", "\n", "
\"FM\"
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will listen to our Berkeley own KPFA 94.1MHz station. KPFA transmitter is located close to Grizzly Peak Road, close to where it meets Claremont Ave.\n", "\n", "\n", "Recall that our SDR does IQ demodulation around a center frequency of choice, $f_d$, followed by a low-pass filter and sampling. This results in a complex digital signal: $$\\qquad \\qquad y_b(t) = Ae^{j2\\pi (f_c-f_d) t + j2\\pi \\Delta f \\int_0^t x(\\tau) d\\tau }.$$\n", "\n", "When $f_d=f_c$ we get:\n", "$$\\qquad \\qquad y_b(t) = Ae^{ j2\\pi \\Delta f \\int_0^t x(\\tau) d\\tau }.$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will implement the following system to demodulate and listen to KPFA:\n", "
\"FM\"
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Signal Capture:\n", "\n", "Like always, you will get the best results if you collect samples outside. \n", "\n", "\n", "Task: \n", "- Set the SDR to sample at a center frequency of 94.1MHz (KPFA) and a rate of 2.4MHz. We will first collect 1*2560000 samples, which is just over 1 seconds of data. We acquire only one second because computing a spectrogram on more data make take a while on some computers. \n", "- Compute and display a spectrogram with a window of 1024 samples. Include the spectrogram in your report. What is the spectral and temporal resolution? Explain what you see. \n", "\n", "Tips:\n", "- Don’t forget to play with different dynamic ranges of the spectrogram for best visualization. (30 worked well for me)\n", "- Make sure to set the gain of the SDR so the signal is not under/overloaded! This will make a HUGE difference when demodulating the carriers." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fs = 2400000\n", "fc = 94.1e6\n", "sdr = RtlSdr()\n", "sdr.sample_rate = fs # sampling rate\n", "sdr.center_freq = fc # center frequency\n", "sdr.gain = 24.0\n", "N_Samples = 2560000\n", "data = sdr.read_samples(N_Samples) # get samples\n", "sdr.close()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Plot. Plotting may take a while!\n", "fig = myspectrogram_hann_ovlp(data, 1024, fs,fc,dbf=30)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Downsampling:\n", "\n", "The bandwidth of an FM station is 200KHz. Since we sampled at 2.4MHz, you should be seeing several stations in the spectrogram: at 93.3MHz, 94.1MHz and 94.9MHz. KPFA, the station we are interested in is in the middle around 94.1MHz. We will need to filter out and select the desired station from the rest of the signals within the spectrum. \n", "\n", "Task:\n", "- Design an FIR filter low-pass filter. The filter should be of length 513 and the cutoff frequency be 100KHz. This will select the frequencies between -100KHz to 100KHz. To design the filter, we will use the function `signal.firwin` (We will talk more about this filter design technique later in class -- but in essence the function creates a windowed sinc function with the appropriate bandwidth for the given length). \n", "Type `signal.firwin?` and execute to get help. \n", "\n", " h = signal.firwin(513,100000.0,nyq=2400000.0/2,window='hanning')\n", "\n", "- Plot the double-sided magnitude frequency response (log scale using `plt.semilogy`) by computing a zero-padded FFT to length 1024 and using fftshift. Use KHz as the x-axis units. \n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Your code here:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, filter the acquired SDR signal with the low-pass filter. Use the function `signal.fftconvolve`, which uses the FFT to perform fast convolution. Plot the spectrogram of the filtered signal. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Your code here:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our signal is now bandlimited with a bandwidth of 200KHz out of 2.4MHz. We can now downsample the signal to reduce computation without aliasing. \n", "\n", "- Decimate the signal by a factor of 10 by selecting every 10th sample, to get a signal representing a rate of 240KHz. After that, you have successfully implemented a downsampler through low-pass filtering and decimation! \n", "\n", "- Plot the spectrogram with a window size of 512. Remember to use the new sampling factor `fs/10.0` for the spectrogram.\n", "- Do you see the rolloff of the low-pass fiter? Is there aliasing?\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Your code here:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that you have gone through this processing with 1 second worth of data. We will acquire a longer segment, which will allow you to appreciate the audio that comes out of it once we demodulate it \n", "\n", "- Acquire 2560000*10 samples, which is just over 10 seconds of data\n", "- Filter and decimate (hence downsample) to the rate of 240KHz\n", "- Plot the spectrogram\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Your code here:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### FM Demodulation:\n", "\n", "The plot does not resemble the broadcast FM spectrum figure above, since it is not yet FM demodulated. We can easily see that the signal is frequency modulated -- because its frequency looks like the time-signal of speech or music.\n", "\n", "The next step we are going to perform is to demodulate the FM signal and look at its spectrogram. For this we need to find the instantaneous frequency as a function of time, which is the time derivative of the phase. Computing the phase and then taking the derivative will be sensitive to phase wraps, which we would like to avoid. \n", "\n", "\n", "Instead, we will take the digital implementation version of the classical FM demodulatation using a limiter, followed by a discriminator. The limiter makes the input have constant amplitude, and the discriminator converts frequency deviations into amplitude modulation. Just as a comment, implementation of an accurate analog disciminator is very difficult whereas implementing a digital one is ridiculously easy!\n", "\n", "##### Limiter\n", "Recall that $y_b(t) = A(t)e^{ j2\\pi (f_c-f_d) t + j2\\pi \\Delta f \\int_0^t x(\\tau) d\\tau }.$ The leading coefficient $A(t)$ is some unwanted amplitude modulation, which can be a result of signal fading or other sources of signal variations. The role of the limiter component in an FM radio is to remove this amplitude modulation so the discriminator will only be sensitive to frequency variations. \n", "\n", "In the digital domain, implementing a limiter is done by simply dividing our signal by its amplitude. \n", "\n", "- Plot the magnitude of the downsampled signal. Can you see the amplitude variation?\n", "- Apply a limiter to your signal. To avoid dividing by 0, divide your signal by its amplitude plus epsilon, where epsilon = 1e-6 is adequate. \n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Your code here:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Discriminator\n", "Assuming that $f_c=f_d$, after the limiter our signal is $ y_b(t)=e^{ j2\\pi \\Delta f \\int_0^t x(\\tau) d\\tau }.$ To get $x(t)$ we can compute:\n", "\n", "$$\\left(\\frac{d}{dt}y_b(t)\\right)y_b^*(t) = j2\\pi\\Delta f \\cdot x(t)\\cdot e^{ j2\\pi \\Delta f \\int_0^t x(\\tau) d\\tau }\\cdot e^{ -j2\\pi \\Delta f \\int_0^t x(\\tau) d\\tau } = j2\\pi\\Delta f \\cdot x(t)$$\n", "\n", "Which gives us our desired signal when we take the imaginary part of the result. \n", "\n", "\n", "We will need to design a differentiator filter, which its frequency response approximates the ideal frequency response, $H_{diff}(e^{jw}) = w$, of a differentiator. Since we are only going to demodulate up to 105KHz, we can have the differentiator roll of to zero after that. \n", "\n", "- Design an FIR differentiator filter with length 31 using the function `signal.remez` which implements an equi-ripple min-max optimal based filter design technique. It requires prescribing frequency bands and their corresponding frequency responses. We need to tell the function that it's a differentiator so it knows to match a linear frequency response within the band of interest. We will prescribe a linear frequency extending from 0-110KHz and then tapering to zero at the Nyquist frequency 120KHz. As we will learn later in class, since the filter is of even order (odd number of coefficients) and is antisymmetric, its Nyquis frequency must be zero!\n", "\n", "Type `signal.remez?` for more information.\n", "\n", "Specifically: \n", "\n", " h_diff = signal.remez(31,[0.0,105000.0,120000.0,120000.0],[1.05/1.2,0],Hz = 240000.0, type='differentiator')\n", " \n", "will try to get a linear frequency response from 0-105000 and zero amplitude at 120000.\n", "\n", "- Plot the filter and its two sided magnitude frequency response (linear scale).\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Your code here:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Demodulate the FM signal by first:\n", "- Filter the signal with the differentiator \n", "- Multiply the result with the conjugate of the original signal\n", "- Take the imaginary component\n", "\n", "Note, that the default implementation of `signal.fftconvolve` will have a delay of 16 samples with respect to the original signal. **To avoid that, use the option `mode='same'`.**\n", "\n", "Note, that after FM demodulating the signal should be real (by taking only the imaginary component) and hence only half the spectrum should be displayed. \n", "\n", "- Plot the spectrogram of the frequency demodulated signal and identify the mono audio, the pilot, the stereo and the RBDS signals. Note, that the RBDS signal may be too weak to detect or may need better spectral resolution. Can you identify the subcarriers? KPFA has two subcarriers, one plays French Hatian radio and the other Punjabi. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Your code here:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Play the mono (Left + Right Channel)\n", "\n", "\n", "The mono signal covers the frequency range of 30Hz-16KHz. However, there are many other signals present. There's another problem. Our sampled signal is at a rate of 240KHz. The soundcard on most modern computers can only deal with a sampling rate of 48KHz. Similarly to the downsampling operation we did before, filter our signal and decimate it before being able to play it. \n", "\n", "- Design a 129 length FIR Bandpass filter with a cuttoff frequency of 16KHz using the command:\n", "\n", " h = signal.firwin(129,16000.0,nyq=240000.0/2, window='hanning')\n", "\n", "- Filter the signal and decimate by a factor of 5 to reduce the rate to 48KHz. Store the result in the a variable called `LPR` (Left + Right). Display the spectrogram of the `LPR`. Use a window length of 256. \n", "\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Your code here:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Playing the Mono (Left + Right channel)\n", "\n", "- Normalize `LPR` so it is in the range between -1 and 1. \n", "- Play it. \n", "\n", "Do you hear radio?\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Your code here:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Write to wav file\n", "write('fm_LR.wav', 48000, LPR_sc)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Demodulating the subcarriers\n", "\n", "There are two other signals in the spectrum, centered at 67 kHz, and 92 kHz. These are subsidiary communications authorization (SCA) services. The FCC specifically doesn't regulate the FM band from 57 kHz up to 100 kHz, only requiring that whatever you transmit there doesn't interfere with the principle FM signal. For a long time this extra spectrum was used for Muzak (a generic name for elevator music), and other targeted FM signals. This has become less common with the advent of the internet, but there are still stations that use these channels. KPFA transmits a Hatian-French subchannel at 67 kHz, and a Punjabi subchannel at 92 kHz.\n", "\n", "Oddly enough, the subchannels themselves are FM encoded. This is like the Russian Matrushka dolls, with one inside another. Fortunately, you have all the tools you need to decode these signals! \n", "\n", "First, demodulate the subcarrier to baseband. We will then low-pass to separate it from the other signals. Then, we will then downsample it to 48KHz, and then frequency demodulate again using the technique described above.\n", "\n", "The entire system is described in the following diagram:\n", "\n", "
\"FM\"
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To demodulate the subcarrier to baseband -- or zero frequency we will need to:\n", "\n", "1. create a time index `t` representing the samples of the signal\n", "2. multiply the frequency demodulated signal in Task 5 with $e^{-2\\pi j f_0 t}$ \n", "\n", "Task: Demodulate the frequency demodulated signal in Task 5 by $f_0=67.65$KHz and plot the spectrogram. Since the signal is now complex, you should plot both sides of the spectrum. The subcarrier should be placed at the DC frequency." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Your code here:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Design a 128 tap low-pass filter with a bandwidth of $\\pm 6$KHz in the passband. \n", "- Filter the signal and decimate by a factor of 5 to get a signal with a sampling rate of 48KHz. \n", "- Plot the spectrogram (make sure you adjusted the spectrogram to represent the new sampling rate!). \n", "\n", "do you see aliasing?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Your code here:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To frequency demodulate, you will need to apply a limiter as before. The discriminator should be implemented in the same way as before. However, you need to design a new filter which approximates a differentiator between $\\pm$8KHz and then tapers to zero. Here's an example:\n", "\n", " h_diff = signal.remez(31,[0.0,8000.0,12000.0,24000],[8.0/24.0,0],Hz = 48000, type='differentiator')\n", "\n", "- Perform the limiter operation\n", "- Apply the discriminator using the differentiator filter method. \n", "- Low-pass the result with a 129 length FIR filter with cutoff of 5KHz to eliminate any residual high frequency noise. Don't forget that the sample rate is 48KHz! \n", "- Scale the result to be within $\\pm1$ and play the audio. \n", "- Plot the spectrogram of the resulting signal. \n", "\n", "Can you hear French Hatian? " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Your code here:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Repeat the procedure to play the subcarrier (Punjabi channel) at $f_0=92$ kHz" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Your code here:" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.11" } }, "nbformat": 4, "nbformat_minor": 0 }