Phy408 Lab 1- Convolution Solved

1 A Discrete Convolution Program

WriteadiscreteconvolutionfunctionmyConvthatconvolvestwoarrays{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. Provideatesttoconvinceyourself(andme)thatyourfunctionagreeswithnumpy.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 thecomparison, 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))

1

2 Simple Physical System: RL Circuit Response

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,

3

Take advantage of your myConv function, or the NumPy built-in function convolve, and writeyourownPythonfunctionV_out = RLresponse(R,L,V_in,dt)totakeaninputseries 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.

Convolution of Synthetic Seismograms (5 pts)

1. Showanalyticallythatitsstepresponse(i.e.,theb(t)weobtainwhentheinputvoltagea(t)= H(t), the Heaviside function) is given byS(t) = e−Rt/L H(t),
   and its impulse response (i.e., the output voltage b(t) when a(t) = δ(t)) is given by

   R(t) = δ(t) − RL e−Rt/L H(t).
   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 d 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

and δ(t) should be discretized as

H = [0.5,1,1,…],

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

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

2

dt

g(t) = √ 1 e−(t/tH )2 πtH

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 sec- onds 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.