{
"metadata": {
"name": ""
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Lab1 - Time Domain Lab\n",
"\n",
"### Written by Miki Lustig and Frank Ong 2014"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This week we will interact with physical time-domain signals. The first task will involve generating and recording sounds on your computer. \n",
"In the second we will use the computer sound system as a simple sonar. The last task will involve sampling the Automatic Dependent Survaillance Broadcast (ADS-B) frequency band (around 1090MHz) using the rtl-sdr and looking at packets sent from airplanes in the bay area. "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Import functions and libraries\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import pyaudio\n",
"from numpy import pi\n",
"from numpy import sin\n",
"from numpy import zeros\n",
"from numpy import r_\n",
"from scipy import signal\n",
"\n",
"# Task II\n",
"import threading,time\n",
"\n",
"# Task IV\n",
"from rtlsdr import RtlSdr\n",
"from numpy import mean\n",
"from numpy import power\n",
"from numpy.fft import fft\n",
"from numpy.fft import fftshift\n",
"from numpy.fft import ifft\n",
"from numpy.fft import ifftshift\n",
"%matplotlib inline"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Part 1: Chirping"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For this assignment you will have to use iPython, and a laptop equipped with a speaker and a microphone. When playing a sound and recording on the computer speaker, the signal goes through several systems. In particular it goes through the response of the speaker, the room we are in and the response of the microphone recording system.\n",
"\n",
"A chirp is a a signal in which the frequency increases linearly with time. In this assignment we will generate a chirp signal and use it to measure the amplitude of the frequency response of our speaker-room-microphone system.\n",
"\n",
"A simultaneous frequency is defined as the derivative of the phase of a signal, $f = \\frac{d\\phi (t)}{2\\pi dt} $. For example, the simultanious frequency of $cos(\\phi(t))=cos(2\\pi f_0 t)$ is $f = \\frac{d\\phi (t)}{2\\pi dt} = f_0$. \n",
"\n",
"\n",
"For a linear chirp, the frequency changes linearly over time. The simultaneous frequency is therefore defined as \n",
"\n",
"$$ f(t) = f_0 + kt. $$\n",
"\n",
"\n",
"So, \n",
"\n",
"$$ x(t) = \\sin(2\\pi\\int_0^t f(t')dt') = \\sin(2\\pi\\int_o^t(f_0+kt')dt') = \\sin(2\\pi(f_0+\\frac{k}{2}t)t) $$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###Task I:\n",
"Generate a 15 seconds long chirp signal, sampled at 44,100[Hz] with a frequency range of 20[Hz] to 20000[Hz]. Set the magnitude of the chirp to 0.5. This will help prevent non-linearities when we play the sound later. \n",
"\n",
"\n",
"1. Set the sample-rate frequency fs = 44100 Hz\n",
"2. Generate a time index from t=0 to t=15 with sampling rate of 44100 Hz\n",
"3. Generate a vector of frequency vs time: f_of_t ($f(t)$) that changes linearly from 20Hz to 20Khz over 15 seconds\n",
"4. Generate a vector of phase vs time: phi_of_t ($\\phi(t) = 2\\pi \\int_0^t f(t)dt$) by numerically integrating f_of_t. You will find the function `np.cumsum` useful.\n",
"5. Generate the chirp function with amplitude of 0.5\n",
" \n",
" \n"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"fs = 44100\n",
"# generate time index\n",
"t = r_[0.0:15*fs:1.0]/fs\n",
"\n",
"# Your code here:"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* Plot the first $\\frac{1}{2}$ a second of the chirp, you will notice that the carrier frequency increases and that the chirp has a constant envelope. To get a nice figure, make sure the aspect ratio of the figure is height/width = 0.2. Label the axis and figure appropriately."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Set the aspect ratio such that the image is wide\n",
"width, height = figaspect(0.2)\n",
"fig = figure(figsize=(width,height))\n",
"\n",
"#Your code below:"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"* Plot the magnitude frequency response of the sequence from 0 to $\\pi$ using the function `signal.freqz`. Note, that the digital frequency range represents a physical frequency range of 0[hz] to 22050[Hz]. To get a nice figure, make sure the aspect ratio of the figure is height/width = 0.2. Label the axis and figure appropriately.\n"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Your code below:\n",
"\n"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* Explain why the chirp is an appropriate signal to measure the magnitude frequency response of a system. \n"
]
},
{
"cell_type": "raw",
"metadata": {},
"source": [
"#### Your answer here:\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###Task II:\n",
"Now, we will play the sound of the chirp on our computer speaker and simultaneously record using the microphone. \n",
"\n",
"* On an Apple computers it is recommended that you turn off the ambient noise reduction by going to system-preferences, selecting sound, choose the input tab and make sure that the \"Use ambient noise reduction\" box is unchecked. In some windows system there's ambient noise reduction as well. Make sure it is also turned off. \n",
"\n",
"* Your laptop most likely has two speakers. It is best if we work only with one. Go to the operating system's sound settings and change the stereo settings such that the speaker that is closest to the microphone is active. Your result will be much better that way. \n",
"\t\t\n",
"* Make sure your output volume is at 70-80% and that the laptop's microphone is on, again to avoid non-linear distorsions. \n",
"\n",
"* We will record 17 seconds just to make sure we capture the entire sequence. \n",
"\n",
"The code below defines some functions to use with pyaudio -- a multi-platform audio python interface. "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"\n",
"import numpy as np\n",
"import threading,time\n",
"import pyaudio\n",
"\n",
"## Define functions that play and record audio\n",
"\n",
"def play_audio( data, p, fs):\n",
" # play_audio plays audio with sampling rate = fs\n",
" # data - audio data array\n",
" # p - pyAudio object\n",
" # fs - sampling rate\n",
" # \n",
" # Example:\n",
" # fs = 44100\n",
" # p = pyaudio.PyAudio() #instantiate PyAudio\n",
" # play_audio( data, p, fs ) # play audio\n",
" # p.terminate() # terminate pyAudio\n",
" \n",
" # open output stream\n",
" ostream = p.open(format=pyaudio.paFloat32, channels=1, rate=int(fs),output=True)\n",
" # play audio\n",
" ostream.write( data.astype(np.float32).tostring() )\n",
" \n",
" \n",
"def record_audio( odata, p, fs, record_seconds ):\n",
" # record_audio records audio with sampling rate = fs\n",
" # odata - output data\n",
" # p - pyAudio object\n",
" # fs - sampling rate\n",
" # record_seconds - record seconds\n",
" #\n",
" # Example:\n",
" # fs = 44100\n",
" # record_seconds = 5\n",
" # odata = zeros( fs * record_seconds ) # initialize odata\n",
" # p = pyaudio.PyAudio() #instantiate PyAudio\n",
" # record_audio( odata, p, fs, record_seconds ) # play audio\n",
" # p.terminate() # terminate pyAudio\n",
" \n",
" # open input stream\n",
" chunk = 1024\n",
" istream = p.open(format=pyaudio.paFloat32, channels=1, rate=int(fs),input=True,frames_per_buffer=chunk)\n",
"\n",
" # record audio in chunks and append to frames\n",
" frames = [];\n",
" for i in range(0, int(fs / chunk * record_seconds)):\n",
" data_str = istream.read(chunk) # read a chunk of data\n",
" data_flt = np.fromstring( data_str, 'float32' ) # convert string to float\n",
" frames.append( data_flt ) # append to list\n",
" \n",
" # flatten list to array\n",
" data = np.concatenate( frames )\n",
" # copy to output\n",
" np.copyto(odata[0:len(data)], data)\n"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### playing and recording audio:\n",
"Run the following code. It is an example how to play and record sound at the same time and uses threading for the play and record threads.\n",
"\n",
"The resulting received sequence will be stored in the variable `rcv_chirp`."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"## Play and record chirp at the same time\n",
"\n",
"RECORD_SECS = 17 # record 10 seconds\n",
"\n",
"# allocate receive chirp response\n",
"rcv_chirp = zeros( RECORD_SECS * fs );\n",
"\n",
"# instantiate PyAudio\n",
"p = pyaudio.PyAudio()\n",
"\n",
"# initialize threads\n",
"t_play = threading.Thread( target = play_audio, args = (s_chirp, p, fs, ))\n",
"t_record = threading.Thread( target = record_audio, args = (rcv_chirp, p, fs, RECORD_SECS , ))\n",
"\n",
"# start recording and playing threads at the same time\n",
"t_record.start()\n",
"t_play.start()\n",
"\n",
"# pause for seconds\n",
"time.sleep( RECORD_SECS +1)\n",
"\n",
"# terminate pyAudio\n",
"p.terminate()\n"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* Plot the frequency response of the received sequence. \n",
"* Also, plot the absolute value of the received signal. Plotting the absolute value (sort of) displays the envelope of the chirp. \n",
"\n",
"Label the figures and use an aspect ration of Height/Width = 0.2\n"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"## Plot chirp response\n",
"# Your code below:\n",
"\n",
"\n"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* Comment on the results you got. In addition, what is the implicit assumption we are making in order to claim that the result is a frequency response? \n",
"(HINT: consider the case when the chirp was very short)"
]
},
{
"cell_type": "raw",
"metadata": {},
"source": [
"#### Answers here:\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Envelope detection with Hilbert transform. \n",
"The absolute value of the of the result \"sort of\" display the envelope, however it still modulated by the (now rectified) frequency sweep carrier. If we write down the response, it can be expressed approximately as $$y[n] = |H[n]|sin(2\\pi (f_0 +k[n*T])nT + \\angle H[n])$$,\n",
"where $|H[n]|$ is the frequency response for the instantaniouse frequency at the nth sample and $\\angle H[n]$ is its phase response. \n",
"The reason that it is only approximate is that there is an inherent assumption that we do not look at transient effects, only steady state effect for each frequency. This is a good approximation because our chirp is very slow compared to the propagation of sound in the room. \n",
"\n",
"One way to get the envelope |H[n]| is to convert it to its analytic signal. The analytic signal $x_a(t)$ of signal $x(t)$ is:\n",
"$$x_a = F^{-1}(F(x)\\cdot 2U) = x + j y$$\n",
"where $F$ is the Fourier transform, $U$ the unit step function,\n",
"and $y$ the Hilbert transform of $x$. In other words, the negative half of the frequency spectrum is zeroed\n",
"out, turning the real-valued signal into a complex signal. This is similar to the question in HW2!\n",
"\n",
"The analytic signal of the received chirp will then be: \n",
"$$ y_a[n] = |H[n]|e^{j2\\pi (f_0 +k[n*T])nT + \\angle H[n]} $$. \n",
"The envelope can be detected by taking the magnitude. \n",
"\n",
"Task: Compute the analytic signal by using the function `signal.hilbert` and plot its absolute value. Note that the discrete hilbert transform is not perfect, since it uses FIR filtering. This will show up as ripple in the envelope.\n",
"\n",
"Label the figures and use an aspect ration of Height/Width = 0.2\n",
"\n"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"## Your lovely code here:\n"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Task III: \n",
"The chirp has a very nice property that its auto-correlation is very narrow. A discrete autocorrelation of a signal is defined as: \n",
"$$ R_{xx}[n] = \\sum_{m=-\\infty}^\\infty x[m]x^*[m-n] = (x[m]*x^*[-m])[n]$$ \n",
"It is basically a cross correlation of the signal with itself. A cross correlation is defined as:\n",
"$$ R_{xy}[n] = \\sum_{m=-\\infty}^\\infty x[m]y^*[m-n] = (x[m]*y^*[-m])[n]$$, \n",
"which is like a convolution, without flipping one of the signals. It can be implemented using a convolution as shown above. \n",
"\n",
"This property is called pulse compression and is used widely in radar. Random noise and some other pseudo-random like sequences also posseses this property. \n",
"\n",
"* Generate a 512 sample chirp pulse with a frequency sweep from 17KHz-19KHz and sampling rate fs=44100. \n",
"* Validate its frequency response by plotting it."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"fs = 44100\n",
"t = r_[0.0:512]/fs\n",
"\n",
"f0 = 17000.0\n",
"f1 = 19000.0\n",
"\n",
"## Your beautiful code here:\n",
"\n"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* Compute the autocorrelation of the chirp discrete convolution, either with `signal.convolve` or `signal.fftconvolve`. Remember that you have to flip the signal since convolution does that already. Use mode=''full''\n",
"* plot the autocorrelation"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"## Your fantastic code here:\n",
"\n"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In a similar way as we did before, it is possible to recover the envelope of the autocorrelation by performing a cross-correlation with the analytic signal and then taking the absolute value. In this case, we know exactly what is the analytic function is!\n",
"\n",
"* Generate s_chirp_a, the analytic function of the chirp by computing: `s_chirp_a = exp(1j* phi_of_t )`. Perform the cross correlation and show the envelope. This is also called a matched filter. \n",
"* Comment on the size of the main lobe of the matched-filter with respect to the size of the cross-correlation. "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"## your amazing code here:\n",
"\n"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "raw",
"metadata": {},
"source": [
"#### Your answer here:\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Repeat Task III for:\n",
" 1. A constant frequency of 19000Hz, 512 samples in length. \n",
" 2. A chirp with a frequency sweep from 15000Hz - 20000Hz, 512 in length. \n",
"\n",
"- Compare the size of the main lobes. How much \"Pulse Compression\" are you getting by using a chirp for detection compared to a single frequency pulse?"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# your nice script to produce beautiful chirps, xcorralations and figures here:\n"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "raw",
"metadata": {},
"source": [
"#### Your answers and interpretations here:\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Dealing with sidelobes\n",
"As you can see, the chirp provides good pulse compression of the main-lobe. However, there are very strong sidelobes. This is because the chirp is multiplied with a rect function, that is abrupt. Instead, we will window the chirp with one of the smooth window functions to taper off the sidelobes. \n",
" \n",
"* Repeat the task above for a chirp with a sweep from 17KHz to 19KHz. This time, multiply the chirp (and its analytic function) with the square-root of the hamming window. You will find the function `np.hamming` useful. The reason for using a square root is that when performing the matched filtering will result in squaring of the frequency response of the result, which gives a multiplication with the window, rathern than squared function of the window. "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# your solution here\n",
"\n",
"\n",
"\n"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### You are now ready to proceed to the Sonar Lab"
]
}
],
"metadata": {}
}
]
}