Phy408 - Lab 1 - Convolution - Problems - Solved

1          A Discrete Convolution Program (5 pts)
Write a discrete convolution function myConv that convolves two arrays {fi,i = 0, . . . , Nf − 1} and {wj, j = 0, . . . , Nw − 1} to obtain an output time series {gn}. For simplicity, assume a fixed sampling interval ∆ = 1, and further, that f and w are 0 outside of their sampled regions.

1.    How long is {gn}? In other words, how many non-zero points can it have? Justify your answer.

2.    Please copy and paste your function g = myConv(f, w) to the PDF report.

3.    Provide a test to convince yourself (and me) that your function agrees with numpy.convolve. For example, generate two random timeseries f, w with Nf = 50, Nw = 100, drawing each element from U[0,1], and plot the difference between your function’s output and numpy’s. Include the code for your test in the PDF report.

4.    Compare the speed of your myConv function to the NumPy function. Provide a plot of the comparison, and include your python code in the PDF report. Is your function faster or slower than the NumPy function? Can you suggest why that is the case?

Hint: For the speed test part, make up your own fi and wj time series, and for simplicity, study the cases of Nf = Nw = 10,100,1000,10000. To accurately time each computation of the convolution function, import the time module and place calls to time.time around your code:

import time t1 = time.time() g = myConv(f, w) t2 = time.time() print(t2-t1)

Alternatively, use the timeit module:

import timeit print(timeit.timeit('g = myConv(f, w)', number=10000))

2          Simple Physical System: RL Circuit Response (8 pts)
Consider a simple physical system consisting of a resistor (with resistance R) and an inductor (with inductance L) in series. We apply an input voltage a(t) across the pair in series, and measure the output voltage b(t) across the inductor alone. For this linear system,

1.    Show analytically that its step response (i.e., the b(t) we obtain when the input voltage a(t) = H(t), the Heaviside function) is given by

S(t) = e−Rt/LH(t),

and its impulse response (i.e., the output voltage b(t) when a(t) = δ(t)) is given by

                                                                                                            − R      Rt/LH(t).

                                                                                  R(t) = δ(t)            e−

L

Hint: Construct and solve the ODE relating the voltages under consideration. Consider the two b(t) choices to derive S(t) and R(t). Formulas dtd H(t) = δ(t) and δ(t)f(t) = δ(t)f(0) may help.

2.    Discretize the impulse response R(t) function, realizing that H(t) should be discretized as

H = [0.5,1,1, . . . ],

and δ(t) should be discretized as

D = [1/dt,0,0, . . . ].

Take advantage of your myConv function, or the NumPy built-in function convolve, and write your own Python function V_out = RLresponse(R,L,V_in,dt) to take an input series Vin sampled at ∆ = dt, and calculate the output series Vout sampled by the same dt. Please paste your Python function here. (Hint: here ∆ may not be 1, so remember to build the multiplication of ∆ into your convolution function.)

3.    Using R = 850Ω, L = 2H, and sampling period dt = 0.20 ms, test your RL-response function with {Hn} series (discretized H(t)) as input, and plot the output time series (as circles) on top of the theoretical curve S(t) given by part 1 (as a solid line). Repeat this for {Dn} (discretized δ(t)) and R(t). Make the time range of the plots 0 to at least 20 ms. Please list your Python code here.

3          Convolution of Synthetic Seismograms (5 pts)
Numerical simulations of seismic wave propagation can now be routinely done for global and regional earthquakes. For a recent southern Pakistan earthquake (Jan 18, 2011, 20:23:33 UTC), raw vertical synthetic seismogram (i.e., displacement field simulated at a seismic station) for station RAYN (Ar Rayn, Saudi Arabia) is provided (RAYN.II.LHZ.sem). A common practice in seismology is to convolve synthetic seismograms with a Gaussian function

g e−(t/tH)2

to reflect either the time duration of the event or the accuracy of the numerical simulation.

1.    Provide two plots. Plot 1: the raw synthetic seismogram for station RAYN between 0 and 800 seconds. Plot 2: Gaussian functions with half duration tH = 10 sec and tH = 20 sec (include a legend). For the gaussians, use the same timestep dt as the seismogram data.

2.    Use numpy’s convolve function to convolve the raw timeseries with a Gaussian function (both tH = 10 and tH = 20 cases). Plot the raw data and the two convolved time series between 0 and 800 seconds on the same graph (include a legend) and comment on your results.

Hints

•      The raw synthetics RAYN.II.LHZ.sem is given as a text file with two columns: time in seconds and displacement in meters.

•      Gaussian functions quickly decay to zero beyond [−3tH,3tH], therefore it is sufficient to sample g(t) within this interval.

•      Use mode='same' when calling numpy convolve to truncate the convolution to the max of the supplied arrays (i.e. length of the raw timeseries in our case). This is convenient, since we want to compare the convolution output to the original timeseries. Alternatively, use the default mode ('full') and truncate the output manually.

•      As a check for part 2, ensure that your convolved timeseries is aligned with (or “overlaps”) the raw data timeseries.