LTspice & NumPy – Part 2: Fast Convolution Filter

Motivation

In the previous post we discussed the possibility to use LTspice as a “plug in” into a Python/Numpy signal processing project. It works quite well: you send a numpy data vector to LTspice, let it run through the simulation and get back a numpy vector again. Everything is abstracted away nicely by the “apply_ltspice_filter.py” module. So far so good.

There is just one problem: It is slow. Every time you process a new signal it takes a few seconds to call up Spice again and to funnel the data through a csv file. Not good for looooong signals or many signals that you want to process with the same filter.

If there only was a way to speed things up …

Appetizer

How about I show you a youtube video with the end results right here, so you stay interested:

Maybe you designed an equalizer for your HiFi system in LTSpice? This is how a wav file would sound if it was played back through the circuit.

How it works

However, there is one trick that we can employ if our filter/circuit is an LTI system (linear time-invariant system). In short, an LTI system is everything that is built with just passive components and ideal (non-saturating) operational amplifiers. Then our filter can be entirely characterized by its impulse response function h(t).

Illustration of an LTI system reacting to various input signals.

The impulse response is obtained by sending an infinitely short pulse with a defined area (~charge) through the circuit and recording its reaction. We can calculate the response of the filter to arbitrary input signals by performing a mathematical procedure called “convolution” of the input signal with the filter’s impulse response.

Yes, that is correct: Doing this trick allows us to calculate what the filter would have done if the input signal were actually processed by it. But it wasn’t 🙂

So we want to use LTspice to get us the impulse response of the filter we designed. Once we have it, we can do all our actual filtering/processing in Python/NumPy by convolution. And performing a convolution can be done quite fast by using Python/SciPy’s on-board signal processing tools (scipy.signal.fftconvolve). As so often, the computational magic lies in the fast fourier transform.

Example I (fast_convolution_filter_demo.py)

Let’s just take the same simple filter again which we used in the original example:

filter_circuit.asc

So far we have just resistors, capacitors, inductors … and wires. Nothing non-linear, nothing that saturates, nothing that changes properties during the simulation. A perfect LTI system. Good.

First, let us repeat the direct, straightforward processing we did last time:

#!/usr/bin/env python3

import numpy as np
from apply_ltspice_filter import \
  apply_ltspice_filter,\
  convolution_filter,\
  get_impulse_response
import matplotlib.pyplot as plt

##################################################
##             generate test signal             ##
##################################################

# our samples shall be 100 ms wide
sample_width=100e-3
# time step between samples: 0.1 ms
delta_t=0.1e-3
samples = int(sample_width/delta_t)

time = np.linspace(0,sample_width,samples)

# we want 1 V between 10 ms and 30 ms, and 2.5 V between 40 and 70 ms
signal_a = 0 + 1*((time > 10e-3) * (time < 30e-3)) + 2.5*((time > 40e-3) * (time < 70e-3))

##################################################
##             apply filter direct              ##
##################################################

# all values in SI units
filter_configuration = {
  "C":100e-6, # 100 uF
  "L":200e-3  # 200 mH
}

dummy, signal_b = apply_ltspice_filter(
      "filter_circuit.asc",
      time, signal_a,
      params = filter_configuration
      )

##################################################
##                     plot                     ##
##################################################

plt.plot(time,signal_a, label="signal_a (LTSpice input)")
plt.plot(time,signal_b, label="signal_b (LTSpice output)")
plt.xlabel("time (s)")
plt.ylabel("voltage (V)")
plt.title("direct processing")
plt.ylim((-1,3.5))

plt.legend()
plt.show()
 

Okay. That looks just as we expected. Also, note that the processing took a few seconds.

Now the new approach. We acquire the impulse response function. As you can see, I wrote a convenient helper function (“get_impulse_response”) to do just that. Inside this function, a very narrow Gaussian pulse is generated with a “charge” of 1 Vs, which is then sent to LTspice and the output is returned. What is “narrow”? Well narrow, compared to the time quantization of our signals, as defined by the above variable “delta_t”.

##################################################
##             get impulse response             ##
##################################################

kernel_delay = 10e-3
kernel_sample_width = 100e-3

kernel_time, kernel = get_impulse_response(
        "filter_circuit.asc",
        params = filter_configuration,
        sample_width = kernel_sample_width,
        delta_t = delta_t,
        kernel_delay = kernel_delay
        )

##################################################
##            plot impulse response             ##
##################################################

plt.plot(kernel_time, kernel, label="impulse response of filter_circuit.asc")
plt.xlabel("time (s)")
plt.ylabel("voltage (V)")
plt.title("impulse response")

plt.legend()
plt.show()
 

Okay. Neat. Also acquiring the response required LTspice to run for a few seconds.

Now that we have the impulse response, we can do the fast convolution filtering:

##################################################
##           apply convolution filter           ##
##################################################

signal_b_conv = convolution_filter(
  signal_a,
  kernel,
  delta_t = delta_t,
  kernel_delay = kernel_delay
  )

##################################################
##      plot direct vs convolution filter       ##
##################################################
  
plt.plot(time,signal_a, label="signal_a (LTSpice input)")
plt.plot(time,signal_b, label="signal_b (LTSpice output)")
plt.plot(time,signal_b_conv, label="signal_b via convolution with IR", linestyle="-.")
plt.xlabel("time (s)")
plt.ylabel("voltage (V)")
plt.ylim((-1,3.5))
plt.title("direct processing vs convolution filter")

plt.legend()
plt.show()
 

Woah. That was fast. And the result looks just like the signal that we actually put through LTspice before.

Now to demonstrate the processing power: Let’s just for fun generate ten different input pulses and process them with the convolution filter:

##################################################
##        demonstration of speed advantage      ##
##################################################

for i in range(0,10):
  ts = i * 6e-3
  signal = 0 + (5+i)*(((time-ts) > 10e-3) * ((time-ts) < (10+2*i)*1e-3)) 
  if i == 0:
    plt.plot(time,signal,color="b",label="in",alpha= 0.9)
  else:
    plt.plot(time,signal,color="b",alpha= 0.9)
  
  signal_conv = convolution_filter(
    signal,
    kernel,
    delta_t = delta_t,
    kernel_delay = kernel_delay
    )
  
  if i == 0:
    plt.plot(time,signal_conv,color="orange",label="out",alpha=0.9)
  else:
    plt.plot(time,signal_conv,color="orange",alpha=0.9)
  
plt.xlabel("time (s)")
plt.ylabel("voltage (V)")
plt.title("demonstration of convolution speed advantage")
plt.ylim((-5,20))
plt.legend()
plt.show()
  

Now you believe me? Ten different input pulses, all processed by the convolution filter in the blink of an eye. Without calling up LTspice at all. 🙂

Download

Find the example on GitHub or download .zip

Example II (fast_conv_wav_filter_demo.py)

Okay. Now we saw that the convolution filter approach is useful for processing a large number of short input signals. But what about one looong input signal, e.g. an audio file? I downloaded a nice example .wav file from freesound.org to test our filter on. Thank you bigmanjoe for creating this sample.

This time we’re going to use a slightly different filter circuit. Something with a more oscillating impulse response:

filter_circuit_2.asc

And here is our example code:

#!/usr/bin/env python3

import numpy as np
import matplotlib.pyplot as plt
from apply_ltspice_filter import \
  apply_ltspice_filter,\
  convolution_filter,\
  get_impulse_response
import scipy.io.wavfile

##################################################
##             read input waveform              ##
##################################################

rate, data = scipy.io.wavfile.read("348275__bigmanjoe__fantasy-orchestra.wav")

# bring to desired numerical format and normalize 
data = data.astype('float32')
# just take left channel
data = data[:,0]
data = data/np.max(np.abs(data))
print("rate: {:d}".format(rate))


##################################################
##             get impulse response             ##
##################################################

delta_t = 1./rate

kernel_delay = 5e-3
kernel_sample_width = 25e-3

# all values in SI units
filter_configuration = {
  "C": 200e-9,  # 200 nF
  "L": 200e-3,  # 200 mH
  "R": 1e3     # 1 kOhm
}

kernel_time, kernel = get_impulse_response(
        "filter_circuit_2.asc",
        params = filter_configuration,
        sample_width = kernel_sample_width,
        delta_t = delta_t,
        kernel_delay = kernel_delay
        )

##################################################
##            plot impulse response             ##
##################################################

plt.plot(kernel_time, kernel, label="impulse response of filter_circuit_2.asc")
plt.xlabel("time (s)")
plt.ylabel("voltage (V)")
plt.title("impulse response")

plt.legend()
plt.show()


##################################################
##      apply IR filter to input waveform       ##
##################################################

filtered = convolution_filter(
  data,
  kernel,
  delta_t = delta_t,
  kernel_delay = kernel_delay
)

##################################################
##        write filtered signal to file         ##
##################################################

scipy.io.wavfile.write("sound_filtered.wav",rate,filtered.astype('float32') )


##################################################
##      plot original and filtered signal       ##
##################################################

time = np.linspace(0,1./rate * len(data) ,len(data) ) # time vector for x axis of plot
plt.plot(time,data,label="original signal",alpha=0.6)
plt.plot(time,filtered, label="filtered signal",alpha=0.6)
plt.xlabel("time (s)")
plt.ylabel("normalized amplitude")
plt.title("waveform before and after filter")
plt.legend()
plt.show()
 
impulse response of filter_circuit_2.asc

So what we did was, process a short sample of 25 ms in LTspice and use the results to filter circa 12 s of audio. Cool eh? 🙂 Now you can scroll to the top and watch the demo video again with the good feeling that you understand how it works.

Download

Find the second example at GitHub or download .zip

Thank You

A warm thank you note to Nuno Brum, who wrote the amazing LTspice RawReader python module and to Alex Stallman and Henk who tested “apply_ltspice_filter” for windows compatibility.

2 responses to “LTspice & NumPy – Part 2: Fast Convolution Filter

  1. Pingback: seamless integration of LTSpice in python/numpy signal processing | a c i d b o u r b o n·

  2. Pingback: LTs[ice and Numpy: a Fast convolution filter #Python #EE « Adafruit Industries – Makers, hackers, artists, designers and engineers!·

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.