```
In [1]:
```# Configure Jupyter so figures appear in the notebook
%matplotlib inline
# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'
# import functions from the modsim.py module
from modsim import *

The following circuit diagram (from Wikipedia) shows a low-pass filter built with one resistor and one capacitor.

A "filter" is a circuit takes a signal, $V_{in}$, as input and produces a signal, $V_{out}$, as output. In this context, a "signal" is a voltage that changes over time.

A filter is "low-pass" if it allows low-frequency signals to pass from $V_{in}$ to $V_{out}$ unchanged, but it reduces the amplitude of high-frequency signals.

By applying the laws of circuit analysis, we can derive a differential equation that describes the behavior of this system. By solving the differential equation, we can predict the effect of this circuit on any input signal.

Suppose we are given $V_{in}$ and $V_{out}$ at a particular instant in time. By Ohm's law, which is a simple model of the behavior of resistors, the instantaneous current through the resistor is:

$ I_R = (V_{in} - V_{out}) / R $

where $R$ is resistance in ohms.

Assuming that no current flows through the output of the circuit, Kirchhoff's current law implies that the current through the capacitor is:

$ I_C = I_R $

According to a simple model of the behavior of capacitors, current through the capacitor causes a change in the voltage across the capacitor:

$ I_C = C \frac{d V_{out}}{dt} $

where $C$ is capacitance in farads (F).

Combining these equations yields a differential equation for $V_{out}$:

$ \frac{d }{dt} V_{out} = \frac{V_{in} - V_{out}}{R C} $

Follow the instructions blow to simulate the low-pass filter for input signals like this:

$ V_{in}(t) = A \cos (2 \pi f t) $

where $A$ is the amplitude of the input signal, say 5 V, and $f$ is the frequency of the signal in Hz.

```
In [2]:
```ohm = UNITS.ohm
farad = UNITS.farad
volt = UNITS.volt
Hz = UNITS.Hz
second = UNITS.second

```
Out[2]:
```

```
In [3]:
```params = Params(
R1 = 1e6 * ohm,
C1 = 1e-9 * farad,
A = 5 * volt,
f = 1000 * Hz)

```
Out[3]:
```

Now we can pass the `Params`

object `make_system`

which computes some additional parameters and defines `init`

.

`omega`

is the frequency of the input signal in radians/second.`tau`

is the time constant for this circuit, which is the time it takes to get from an initial startup phase to`cutoff`

is the cutoff frequency for this circuit (in Hz), which marks the transition from low frequency signals, which pass through the filter unchanged, to high frequency signals, which are attenuated.`t_end`

is chosen so we run the simulation for 4 cycles of the input signal.

```
In [4]:
```def make_system(params):
"""Makes a System object for the given conditions.
params: Params object
returns: System object
"""
f, R1, C1 = params.f, params.R1, params.C1
init = State(V_out = 0)
omega = 2 * np.pi * f
tau = R1 * C1
cutoff = 1 / R1 / C1 / 2 / np.pi
t_end = 4 / f
dt = t_end / 4000
return System(params, init=init, t_end=t_end, dt=dt,
omega=omega, tau=tau, cutoff=cutoff.to(Hz))

Let's make a `System`

```
In [5]:
```system = make_system(params)

```
Out[5]:
```

**Exercise:** Write a slope function that takes as an input a `State`

object that contains `V_out`

, and returns the derivative of `V_out`

.

Note: The ODE solver requires the return value from slope_func to be a sequence, even if there is only one element. The simplest way to do that is to return a list with a single element:

` return [dV_out_dt]`

```
In [44]:
```# Solution
def slope_func(state, t, system):
"""Compute derivatives of the state.
state: V_out
t: time
system: System object with A, omega, R1 and C1
returns: dV_out/dt
"""
[V_out] = state
R1, C1 = system.R1, system.C1
A, omega = system.A, system.omega
V_in = A * np.cos(omega * t)
V_R1 = V_in - V_out
I_R1 = V_R1 / R1
I_C1 = I_R1
dV_out_dt = I_C1 / C1
return [dV_out_dt]

Test the slope function with the initial conditions.

```
In [45]:
```slope_func(system.init, 0*UNITS.s, system)

```
Out[45]:
```

`t_eval=ts`

to make sure we have enough data points to plot and analyze the results.

```
In [46]:
```results, details = run_ode_solver(system, slope_func)
details

```
Out[46]:
```

```
In [47]:
```results.head()

```
Out[47]:
```

Here's a function you can use to plot `V_out`

as a function of time.

```
In [48]:
```def plot_results(results):
xs = results.V_out.index
ys = results.V_out.values
t_end = get_last_label(results)
if t_end < 10:
xs *= 1000
xlabel = 'Time (ms)'
else:
xlabel = 'Time (s)'
plot(xs, ys)
decorate(xlabel=xlabel,
ylabel='$V_{out}$ (volt)',
legend=False)
plot_results(results)

```
```

If things have gone according to plan, the amplitude of the output signal should be about 0.8 V.

Also, you might notice that it takes a few cycles for the signal to get to the full amplitude.

```
In [11]:
```fs = [1, 10, 100, 1000, 10000, 100000] * Hz
for i, f in enumerate(fs):
system = make_system(Params(params, f=f))
results, details = run_ode_solver(system, slope_func)
subplot(3, 2, i+1)
plot_results(results)

```
```

`tau`

, which is 1 ms.

```
In [12]:
```system = make_system(Params(params, f=1000*Hz))
results, details = run_ode_solver(system, slope_func)
V_out = results.V_out
plot_results(results)

```
```

The following function computes `V_in`

as a `TimeSeries`

:

```
In [13]:
```def compute_vin(results, system):
"""Computes V_in as a TimeSeries.
results: TimeFrame with simulation results
system: System object with A and omega
returns: TimeSeries
"""
A, omega = system.A, system.omega
ts = results.index.values * UNITS.second
V_in = A * np.cos(omega * ts)
return TimeSeries(V_in, results.index, name='V_in')

```
In [14]:
```V_in = compute_vin(results, system)
plot(V_out)
plot(V_in)
decorate(xlabel='Time (s)',
ylabel='V (volt)')

```
```

```
In [15]:
```def estimate_A(series):
"""Estimate amplitude.
series: TimeSeries
returns: amplitude in volts
"""
return (series.max() - series.min()) / 2

`V_in`

should be near 5 (but not exact because we evaluated it at a finite number of points).

```
In [16]:
```A_in = estimate_A(V_in)

```
Out[16]:
```

The amplitude of `V_out`

should be lower.

```
In [17]:
```A_out = estimate_A(V_out)

```
Out[17]:
```

And here's the ratio between them.

```
In [18]:
```ratio = A_out / A_in

```
Out[18]:
```

```
In [19]:
```ratio.to_base_units()

```
Out[19]:
```

**Exercise:** Encapsulate the code we have so far in a function that takes two TimeSeries objects and returns the ratio between their amplitudes.

```
In [20]:
```# Solution
def estimate_ratio(V1, V2):
"""Estimate the ratio of amplitudes.
V1: TimeSeries
V2: TimeSeries
returns: amplitude ratio
"""
a1 = estimate_A(V1)
a2 = estimate_A(V2)
return a1 / a2

And test your function.

```
In [21]:
```estimate_ratio(V_out, V_in)

```
Out[21]:
```

The delay between the peak of the input and the peak of the output is call a "phase shift" or "phase offset", usually measured in fractions of a cycle, degrees, or radians.

To estimate the phase offset between two signals, we can use cross-correlation. Here's what the cross-correlation looks like between `V_out`

and `V_in`

:

```
In [22]:
```corr = correlate(V_out, V_in, mode='same')
corr = TimeSeries(corr, V_in.index)
plot(corr, color='C4')
decorate(xlabel='Lag (s)',
ylabel='Correlation')

```
```

```
In [23]:
```peak_time = corr.idxmax() * UNITS.second

```
Out[23]:
```

We can express the phase offset as a multiple of the period of the input signal:

```
In [24]:
```period = 1 / system.f

```
Out[24]:
```

```
In [25]:
```(peak_time / period).to_reduced_units()

```
Out[25]:
```

We don't care about whole period offsets, only the fractional part, which we can get using `modf`

:

```
In [26]:
```frac, whole = np.modf(peak_time / period)
frac = frac.to_reduced_units()

```
Out[26]:
```

Finally, we can convert from a fraction of a cycle to degrees:

```
In [27]:
```frac * 360 * UNITS.degree

```
Out[27]:
```

**Exercise:** Encapsulate this code in a function that takes two `TimeSeries`

objects and a `System`

object, and returns the phase offset in degrees.

Note: by convention, if the output is shifted to the right, the phase offset is negative.

```
In [28]:
```# Solution
def estimate_offset(V1, V2, system):
"""Estimate phase offset.
V1: TimeSeries
V2: TimeSeries
system: System object with f
returns: amplitude ratio
"""
corr = correlate(V1, V2, mode='same')
corr = TimeSeries(corr, V1.index)
peak_time = corr.idxmax() * UNITS.second
period = 1 / system.f
frac, whole = np.modf(peak_time / period)
frac = frac.to_reduced_units()
return -frac * 360 * UNITS.degree

Test your function.

```
In [29]:
```estimate_offset(V_out, V_in, system)

```
Out[29]:
```

**Exercise:** Write a function that takes as parameters an array of input frequencies and a `Params`

object.

For each input frequency it should run a simulation and use the results to estimate the output ratio (dimensionless) and phase offset (in degrees).

It should return two `SweepSeries`

objects, one for the ratios and one for the offsets.

```
In [30]:
```# Solution
def sweep_frequency(fs, params):
ratios = SweepSeries()
offsets = SweepSeries()
for i, f in enumerate(fs):
system = make_system(Params(params, f=f))
results, details = run_ode_solver(system, slope_func)
V_out = results.V_out
V_in = compute_vin(results, system)
f = magnitude(f)
ratios[f] = estimate_ratio(V_out, V_in)
offsets[f] = estimate_offset(V_out, V_in, system)
return ratios, offsets

Run your function with these frequencies.

```
In [31]:
```fs = 10 ** linspace(0, 4, 9) * Hz

```
Out[31]:
```

```
In [32]:
```ratios, offsets = sweep_frequency(fs, params)

We can plot output ratios like this:

```
In [33]:
```plot(ratios, color='C2', label='output ratio')
decorate(xlabel='Frequency (Hz)',
ylabel='$V_{out} / V_{in}$')

```
```

```
In [34]:
```def plot_ratios(ratios, system):
"""Plot output ratios.
"""
# axvline can't handle a Quantity with units
cutoff = magnitude(system.cutoff)
plt.axvline(cutoff, color='gray', alpha=0.4)
plot(ratios, color='C2', label='output ratio')
decorate(xlabel='Frequency (Hz)',
ylabel='$V_{out} / V_{in}$',
xscale='log', yscale='log')

```
In [35]:
```plot_ratios(ratios, system)

```
```

This plot shows the cutoff behavior more clearly. Below the cutoff, the output ratio is close to 1. Above the cutoff, it drops off linearly, on a log scale, which indicates that output ratios for high frequencies are practically 0.

Here's the plot for phase offset, on a log-x scale:

```
In [36]:
```def plot_offsets(offsets, system):
"""Plot phase offsets.
"""
# axvline can't handle a Quantity with units
cutoff = magnitude(system.cutoff)
plt.axvline(cutoff, color='gray', alpha=0.4)
plot(offsets, color='C9')
decorate(xlabel='Frequency (Hz)',
ylabel='Phase offset (degree)',
xscale='log')

```
In [37]:
```plot_offsets(offsets, system)

```
```

For low frequencies, the phase offset is near 0. For high frequencies, it approaches 90 degrees.

By analysis we can show that the output ratio for this signal is

$A = \frac{1}{\sqrt{1 + (R C \omega)^2}}$

where $\omega = 2 \pi f$, and the phase offset is

$ \phi = \arctan (- R C \omega)$

**Exercise:** Write functions that take an array of input frequencies and returns $A(f)$ and $\phi(f)$ as `SweepSeries`

objects. Plot these objects and compare them with the results from the previous section.

```
In [38]:
```# Solution
def output_ratios(fs, system):
R1, C1, omega = system.R1, system.C1, system.omega
omegas = 2 * np.pi * fs
rco = R1 * C1 * omegas
A = 1 / np.sqrt(1 + rco**2)
return SweepSeries(A, magnitude(fs))

Test your function:

```
In [39]:
```A = output_ratios(fs, system)

```
Out[39]:
```

```
In [40]:
```# Solution
def phase_offsets(fs, system):
R1, C1, omega = system.R1, system.C1, system.omega
omegas = 2 * np.pi * fs
rco = R1 * C1 * omegas
phi = np.arctan(-rco).to(UNITS.degree)
return SweepSeries(phi, magnitude(fs))

Test your function:

```
In [41]:
```phi = phase_offsets(fs, system)

```
Out[41]:
```

Plot the theoretical results along with the simulation results and see if they agree.

```
In [42]:
```plot(A, ':', color='gray')
plot_ratios(ratios, system)

```
```

```
In [43]:
```plot(phi, ':', color='gray')
plot_offsets(offsets, system)

```
```

**Exercise:** Consider modifying this notebook to model a first order high-pass filter, a two-stage second-order low-pass filter, or a passive band-pass filter.

```
In [ ]:
```