2

This example is taken from a tutorial related to convolution integral.

I would like to export this example as animation in mp4 format. So far, the code looks like this :

import scipy.integrate
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def showConvolution(f1, f2, t0):
    # Calculate the overall convolution result using Simpson integration
    convolution = np.zeros(len(t))
    for n, t_ in enumerate(t):
        prod = lambda tau: f1(tau) * f2(t_-tau)
        convolution[n] = scipy.integrate.simps(prod(t), t)

    # Create the shifted and flipped function
    f_shift = lambda t: f2(t0-t)
    prod = lambda tau: f1(tau) * f2(t0-tau)

    # Plot the curves
    plt.gcf().clear() # il

    plt.subplot(211)
    plt.gca().set_ymargin(0.05) # il
    plt.plot(t, f1(t), label=r'$f_1(\tau)$')
    plt.plot(t, f_shift(t), label=r'$f_2(t_0-\tau)$')
    plt.fill(t, prod(t), color='r', alpha=0.5, edgecolor='black', hatch='//') # il
    plt.plot(t, prod(t), 'r-', label=r'$f_1(\tau)f_2(t_0-\tau)$')
    plt.grid(True); plt.xlabel(r'$\tau$'); plt.ylabel(r'$x(\tau)$') # il
    plt.legend(fontsize=10) # il
    plt.text(-4, 0.6, '$t_0=%.2f$' % t0, bbox=dict(fc='white')) # il

    # plot the convolution curve
    plt.subplot(212)
    plt.gca().set_ymargin(0.05) # il
    plt.plot(t, convolution, label='$(f_1*f_2)(t)$')

    # recalculate the value of the convolution integral at the current time-shift t0
    current_value = scipy.integrate.simps(prod(t), t)
    plt.plot(t0, current_value, 'ro')  # plot the point
    plt.grid(True); plt.xlabel('$t$'); plt.ylabel('$(f_1*f_2)(t)$') # il
    plt.legend(fontsize=10) # il
    plt.show() # il

Fs = 50  # our sampling frequency for the plotting
T = 5    # the time range we are interested in
t = np.arange(-T, T, 1/Fs)  # the time samples
f1 = lambda t: np.maximum(0, 1-abs(t))
f2 = lambda t: (t>0) * np.exp(-2*t)

t0 = np.arange(-2.0,2.0, 0.05)

fig = plt.figure(figsize=(8,3))
anim = animation.FuncAnimation(fig, showConvolution(f1,f2, t0), frames=np.linspace(0, 50, 500), interval=80)

anim.save('animation.mp4', fps=30) # fps = frames per second

plt.show()

As I understood, I should be able to change t0 value between -2.00 and 2.00 with 0.05 steps. At first glance I've tried to use numpy's arange function.

t0 = np.arange(-2.0,2.0, 0.05)

But it gives an error message :

ValueError: operands could not be broadcast together with shapes (80,) (500,)

How should I change t0 value so that I should be able to generate animation video?

Edit : I tried the suggested changes. I run this example with

python convolution.py

Rather than an animation I see the output of see the output of convolution integral at t0 = -0.20.

convolution integral

Is there a way to change t0 so that I'd like to able to save it as animation as in the tutorial In the example t0 decreases from -2.0 to -1.95 etc the green curve is shifted right, and the area between curves, product increases. In the example there is an html animation and I would like to save as mp4 file.

Edit 2 :

Removing the plt.show() calls from inside of the re-draw function allows it to run end-to-end and write out the animation.

2 Answers2

1

It seems that the example is wrong.

The second argument in FuncAnimation takes a callable of which the first argument will get a new value from the 'frames' keyword argument every loop. Please take a look at the matplotlib documentation to find more information about the callable's required signature.

I simply changed the showConvolution() arguments around so that t0 is the first argument. The range t0 is used as the desired frames argument. The lambda functions f1 and f2 are passed in a tuple in 'fargs'.

Hope it helps you,

Cheers,

BdeG

import scipy.integrate
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def showConvolution(t0,f1, f2):
    # Calculate the overall convolution result using Simpson integration
    convolution = np.zeros(len(t))
    for n, t_ in enumerate(t):
        prod = lambda tau: f1(tau) * f2(t_-tau)
        convolution[n] = scipy.integrate.simps(prod(t), t)

    # Create the shifted and flipped function
    f_shift = lambda t: f2(t0-t)
    prod = lambda tau: f1(tau) * f2(t0-tau)

    # Plot the curves
    plt.gcf().clear() # il

    plt.subplot(211)
    plt.gca().set_ymargin(0.05) # il
    plt.plot(t, f1(t), label=r'$f_1(\tau)$')
    plt.plot(t, f_shift(t), label=r'$f_2(t_0-\tau)$')
    plt.fill(t, prod(t), color='r', alpha=0.5, edgecolor='black', hatch='//') # il
    plt.plot(t, prod(t), 'r-', label=r'$f_1(\tau)f_2(t_0-\tau)$')
    plt.grid(True); plt.xlabel(r'$\tau$'); plt.ylabel(r'$x(\tau)$') # il
    plt.legend(fontsize=10) # il
    plt.text(-4, 0.6, '$t_0=%.2f$' % t0, bbox=dict(fc='white')) # il

    # plot the convolution curve
    plt.subplot(212)
    plt.gca().set_ymargin(0.05) # il
    plt.plot(t, convolution, label='$(f_1*f_2)(t)$')

    # recalculate the value of the convolution integral at the current time-shift t0
    current_value = scipy.integrate.simps(prod(t), t)
    plt.plot(t0, current_value, 'ro')  # plot the point
    plt.grid(True); plt.xlabel('$t$'); plt.ylabel('$(f_1*f_2)(t)$') # il
    plt.legend(fontsize=10) # il
    plt.show() # il

Fs = 50  # our sampling frequency for the plotting
T = 5    # the time range we are interested in
t = np.arange(-T, T, 1/Fs)  # the time samples
f1 = lambda t: np.maximum(0, 1-abs(t))
f2 = lambda t: (t>0) * np.exp(-2*t)

t0 = np.arange(-2.0,2.0, 0.05)

fig = plt.figure(figsize=(8,3))
anim = animation.FuncAnimation(fig, showConvolution, frames=t0, fargs=(f1,f2),interval=80)

anim.save('animation.mp4', fps=30) # fps = frames per second

plt.show()
BdeG
  • 76
  • 3
0

I have a quick implementation to visualize the convolution of two real functions defined with Sympy. The code outputs a .GIF though.

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt    
from sympy import lambdify
from matplotlib.animation import FuncAnimation

def symplot(t, F, interval, funLabel):
'''
Create plots of sympy symbolic functions

:param t: sympy variable
:param F: sympy function F(t)
:param interval: array of values of t where F should be evaluated [np.array]
:funLabel: curve label be displayed in the plot [string].    
'''
fig = plt.figure()
if type(F) == list:
    indLabel = 0
    for f in F:
        func  = lambdify(t, f, modules=['numpy', {'Heaviside': lambda t:np.heaviside(t,0)}])
        f_num = func(interval)
        
        plt.plot(interval, f_num, label=funLabel[indLabel])
        plt.legend();
        plt.xlim([min(interval), max(interval)]);
        plt.xlabel('tempo [s]');
        indLabel += 1
else:        
    func  = lambdify(t, F, modules=['numpy', {'Heaviside': lambda t:np.heaviside(t,0)}])
    f_num = func(interval)           
                
    plt.plot(interval, f_num, label=funLabel)
    plt.legend(loc="upper right");
    plt.xlim([min(interval), max(interval)]);
    plt.xlabel('tempo [s]');            

plt.grid();
plt.close();
return fig

def genConvGIF(x, h, t, totalTime, ti, tf, figName, xlabel=[], ylabel=[], fram=200, inter=20, plotConv=False):    
'''
Create and save a convolution plot animation as GIF
 
:param x: x(t) function [sympy expr]
:param h: h(t) function [sympy expr]
:param t: t time variable [sympy variable]
:param totalTime: array of time instants where the functions will be evaluated [nparray]
:param ti: time when animation starts [scalar]
:param tf: time when animation stops [scalar]
:param figName: figure file name w/ folder path [string]
:param xlabel: xlabel [string]
:param ylabel: ylabel [string]
:param fram: number of frames [int]
:param inter: time interval between frames [ms]

'''
x_func = lambdify(t, x, modules=['numpy', {'Heaviside': lambda t:np.heaviside(t,0)}])
h_func = lambdify(t, h, modules=['numpy', {'Heaviside': lambda t:np.heaviside(t,0)}])

x_num  = x_func(totalTime)
h_num  = h_func(totalTime)    
dt = totalTime[1]-totalTime[0]

if plotConv:
    y_num  = np.convolve(h_num, x_num, 'same')*dt
    ymax = np.max([x_num, h_num, y_num])
    ymin = np.min([x_num, h_num, y_num])
else:
    ymax = np.max([x_num, h_num])
    ymin = np.min([x_num, h_num])

figAnim = plt.figure()
ax      = plt.axes(xlim=(totalTime.min(), totalTime.max()),ylim=(ymin-0.1*np.abs(ymax), ymax+0.1*np.abs(ymax)))
line1, line2, line3 = ax.plot([], [], [], [], [], [])
line1.set_label(ylabel[0])
line2.set_label(ylabel[1])

if plotConv:
    line3.set_label(ylabel[2])
    
ax.grid()
ax.legend(loc="upper right");

# plot static function
figh = symplot(t, h, totalTime, 'h(t)')

if len(xlabel): 
       plt.xlabel(xlabel)            
    
def init():        
    line1.set_data(figh.get_axes()[0].lines[0].get_data())
    return line1,

plt.close(figh)

delays = totalTime[::int(len(totalTime)/fram)]
ind    = np.arange(0, len(totalTime), int(len(totalTime)/fram))

ind    = ind[delays > ti]
delays = delays[delays > ti]

ind    = ind[delays < tf]   
delays = delays[delays < tf]

totalFrames = len(delays)

def animate(i):
    figx = symplot(t, x.subs({t:delays[i]-t}), totalTime, 'x(t-τ)')
    line2.set_data(figx.get_axes()[0].lines[0].get_data())
    
    if plotConv:
        line3.set_data(totalTime[0:ind[i]], y_num[0:ind[i]])
        
    plt.close(figx)
    return line2, line3

anim = FuncAnimation(figAnim, animate, init_func=init, frames=totalFrames, interval=inter, blit=True)

anim.save(figName, dpi=200, writer='imagemagick')
plt.close()

example of .GIF generated by the code

You can find a notebook with a quick tutorial of how to use it here:

https://github.com/edsonportosilva/ElectricCircuits/blob/master/Jupyter%20notebooks/Visualize%20convolution.ipynb

David Buck
  • 3,752
  • 35
  • 31
  • 35
epsilva
  • 1
  • 2