538 Day -3: Trigonometric interpolation and spectral differentiation¶

Visualization of the matrix of the discrete Fourier transform¶

In [1]:
from numpy import *
from matplotlib import rcdefaults
rcdefaults()  # restore default matplotlib rc parameters
%config InlineBackend.figure_format='retina'   
import seaborn as sns  # wrapper for matplotlib that provides prettier styles and more
import matplotlib.pyplot as plt  # use matplotlib functionality directly
%matplotlib inline  
sns.set()
In [2]:
from numpy import linspace,ones_like,real,imag,cos,sin,pi,exp,arange
import matplotlib.pyplot as plt  
from matplotlib import cm

def fft_matrix_viz(n):
    L = 1.
    vsep = -1.
    markersize = 150/float(n)# 75./float(n)
    i = complex(0.,1.)
    t = linspace(0,L,501,endpoint=True)
    kvals = range(0,n)#-n//2,n)
    colormap = cm.get_cmap('hsv') # a cyclic colormap
    nc = len(kvals)
    colors = colormap( linspace(0,1,nc,endpoint=False) )  # get nc colors from the colormap
    # pastelize colors
    p = 0.5
    colors = (1-p)*colors + p*ones_like(colors)
    # darken colors
    d = 0.9
    colors *= [d,d,d,1]
    fontsize = 20*4/n
    plt.figure(figsize=(12, 12))

    theta = linspace(0,2*pi,100)
    c = cos(theta)
    s = sin(theta)
    def plotZ(x,y,Z):
        radius = 0.3
        plt.plot(x+radius*c, y+radius*s,'k',alpha=0.35)  # draw unit circle
        plt.plot([x,x+radius*real(Z)],[y,y+radius*imag(Z)],'k',alpha=0.35) # draw radius
        plt.plot(x+radius*real(Z),y+radius*imag(Z),'o',markersize=markersize,color=colors[(j*k)%len(kvals)]) # mark Z

    plt.subplot(1,1,1,aspect=1,frameon=False)
    for k in kvals:
        zd = exp(-2*pi*i*k*arange(n)/float(n))
        for j,Z in enumerate(zd): plotZ(j,vsep*k,Z)
        plt.text(-0.5,vsep*k,str(k),va='center',ha='right',fontsize=fontsize)

    for j in range(n):
        plt.text(j,vsep*(min(kvals)-0.45),str(j),ha='center',fontsize=fontsize)
    plt.xlim(-.5,n+.5)
    plt.ylim(vsep*(max(kvals)+1),vsep*(min(kvals)-1))
    plt.xticks([])
    plt.yticks([])

    plt.savefig('dft_visualization_'+'n'+str(n).zfill(4)+'.svg',bbox_inches='tight')

fft_matrix_viz(4)
In [5]:
fft_matrix_viz(8)
In [6]:
fft_matrix_viz(20)

Small example of DFT and IDFT¶

In [7]:
import numpy as np
from scipy.fftpack import fft, ifft
x = [1,2,3,5]
y = fft(x)
y
Out[7]:
array([11.-0.j, -2.+3.j, -3.-0.j, -2.-3.j])
In [8]:
ifft(y)
Out[8]:
array([1.+0.j, 2.+0.j, 3.-0.j, 5.+0.j])
In [9]:
## the matrix?
fft(np.eye(4))
Out[9]:
array([[ 1.-0.j,  1.+0.j,  1.-0.j,  1.-0.j],
       [ 1.-0.j,  0.-1.j, -1.-0.j,  0.+1.j],
       [ 1.-0.j, -1.+0.j,  1.-0.j, -1.-0.j],
       [ 1.-0.j,  0.+1.j, -1.-0.j,  0.-1.j]])

Another visualization of the DFT matrix¶

In [3]:
n=16
def fft_matrix_viz2(n):
    L = 1.
    vsep = -4.
    markersize = 150./float(n)
    i = complex(0.,1.)
    t = linspace(0,L,501,endpoint=True)
    kvals = range(n)#-n//2,n)

    plt.figure(figsize=(12, 12))

    recolor,imcolor = '#4444ff','#ff4444'
    plt.subplot(1,1,1)
    for k in kvals:
        zd = exp(-2*pi*i*k*arange(n)/float(n))
        zc = exp(-2*pi*i*k*t)
        plt.plot([0,n],[vsep*k,vsep*k],'#dddddd')
        plt.plot(t*n,vsep*k+real(zc),recolor,alpha=0.25)
        plt.plot(t*n,vsep*k+imag(zc),imcolor,alpha=0.25)
        plt.plot(arange(n),vsep*k+real(zd),'.',color=recolor,markersize=markersize)
        plt.plot(arange(n),vsep*k+imag(zd),'.',color=imcolor,markersize=markersize)
        plt.text(-0.75,vsep*k,str(k),va='center')

    plt.title('$F_{'+str(n)+'}$ : real part blue, imag part red')
    plt.xlim(-.5,n+.5)
    plt.ylim(vsep*(max(kvals)+1),vsep*(min(kvals)-1))
    plt.xticks(range(n))
    plt.yticks([])

    plt.savefig('dft_visualization2_'+'n'+str(n).zfill(4)+'.png')

fft_matrix_viz2(4)
In [11]:
fft_matrix_viz2(8)
In [4]:
fft_matrix_viz2(20)

Example of DFT of a bigger vector¶

The vector is a recording of a piano note piano_low_f.wav (link fixed)

In [12]:
from scipy.io import wavfile
from numpy import *

rate,x = wavfile.read('piano_low_f.wav')
print('sampling rate',rate)
print('x.shape',x.shape)

n = x.shape[0]
d = n/float(rate)  # duration of recording
print('duration',d,'secs')
#assert(0)
import matplotlib.pyplot as plt
#plt.plot(x)

from scipy.fftpack import fft

y = fft(x)

plt.figure(figsize=(20,5))
# mark the nominal frequency and its harmonics (integer multiples)
f = 43.6535
for i in range(20):
        #plt.plot([i*f,i*f],[0,1.e9],'r')
        plt.axvline(i*f,color='r')

plt.plot( arange(n)/d , abs(y),'o-' )

plt.xlim(0,3000/d)
plt.xlabel('frequency (Hz) (proportional to index k)');
plt.ylabel('$|y_k|$');

sampling rate 44100
x.shape (162030,)
duration 3.6741496598639456 secs

new recording of E above middle C, nominally 329.628 Hz¶

In [21]:
from scipy.io import wavfile
from numpy import *

rate,x = wavfile.read('250501_0230_cropped.wav')   # E above middle C, nominally 329.628 Hz
print('sampling rate',rate)
print('x.shape',x.shape)

n = x.shape[0]
d = n/float(rate)  # duration of recording
print('duration',d,'secs')
#assert(0)
import matplotlib.pyplot as plt
#plt.plot(x)

from scipy.fftpack import fft
print(x.shape)
y = fft(x[:,0])

plt.figure(figsize=(20,5))
# mark the nominal frequency and its harmonics (integer multiples)
f = 329.6
for i in range(20):
        #plt.plot([i*f,i*f],[0,1.e9],'r')
        plt.axvline(i*f,color='r')

plt.plot( arange(n)/d , abs(y),'o-' )

plt.xlim(0,6000/d)
plt.xlabel('frequency (Hz) (proportional to index k)');
plt.ylabel('$|y_k|$');

sampling rate 96000
x.shape (116735, 2)
duration 1.2159895833333334 secs
(116735, 2)

Trigonometric interpolation¶

Visualization of IDFT matrix¶

In [25]:
n=16
def ifft_matrix_viz2(n):
    L = 1.
    vsep = -4.
    markersize = 150./float(n)
    i = complex(0.,1.)
    t = linspace(0,L,501,endpoint=True)
    kvals = range(-n//2+1,n//2+1) #(2*n,3*n)

    plt.figure(figsize=(12, 12))

    recolor,imcolor = '#4444ff','#ff4444'
    plt.subplot(1,1,1)
    for k in kvals:
        zd = exp(-2*pi*i*k*arange(n)/float(n))
        zc = exp(-2*pi*i*k*t)
        plt.plot([vsep*k,vsep*k],[0,n],'#dddddd')
        plt.plot(vsep*k+real(zc),t*n,recolor,alpha=0.25)
        plt.plot(vsep*k+imag(zc),t*n,imcolor,alpha=0.25)
        plt.plot(vsep*k+real(zd),arange(n),'.',color=recolor,markersize=markersize)
        plt.plot(vsep*k+imag(zd),arange(n),'.',color=imcolor,markersize=markersize)
        plt.text(vsep*k,-0.75,str(k),ha='center')

    plt.title('$F_{'+str(n)+'}$ : real part blue, imag part red')
    plt.ylim(n+.5,-.5)
    plt.xlim(vsep*(min(kvals)-1),vsep*(max(kvals)+1))
    plt.yticks(range(n))
    plt.xticks([])

    plt.savefig('idft_visualization2_'+'n'+str(n).zfill(4)+'.png')

ifft_matrix_viz2(6)
In [29]:
from numpy import *
from matplotlib import rcdefaults
rcdefaults()  # restore default matplotlib rc parameters
%config InlineBackend.figure_format='retina'   
import seaborn as sns  # wrapper for matplotlib that provides prettier styles and more
import matplotlib.pyplot as plt  # use matplotlib functionality directly
%matplotlib inline  
sns.set()
In [30]:
%%javascript
IPython.OutputArea.prototype._should_scroll = function(lines) { return false; }
In [34]:
#triginterp.py
import numpy as np
from scipy.fftpack import fft,ifft
from pylab import plot,ion,ioff,draw,title,figure,axvline,legend

formula = 'R'

for example in [2,3]:

        figure(example,figsize=(15,10))
        if example==1:
                mytitle = 'vector [1,0,-1,0]'
                c = 0.
                d = 1.
                L = d-c
                x = np.array([1.,0.,-1.,0.])+0.5
                n = len(x)
                t = L*np.arange(n)/float(n)
        elif example==2:
                mytitle = 'distorted sine'
                n = 8
                c = 0.
                d = 1.
                L = d-c
                t = np.linspace(0.,L,n,endpoint=False)
                #print( len(t) )
                r = 2
                def myf(t): return np.exp(r*np.sin(2*pi*t/L))#sin(2*pi*t/L)#
                x = myf(t)
        elif example==3:
                mytitle = 'Runge function'
                n = 12#8
                c = -1.
                d =  1.
                L = d-c
                t = np.linspace(c,d,n,endpoint=False)
                #print( len(t) )
                def myf(t): return 1/(1+12.*t*t)
                x = myf(t)

        elif example==4:
                mytitle = 'odd sine'
                n = 8
                c = -1.
                d =  1.
                L = d-c
                t = linspace(c,d,n,endpoint=False)
                print( len(t) )
                def myf(t): return np.sin(3*pi*(t-c)/L)
                x = myf(t)

        elif example==5:
                mytitle = 'vector [0,1,2,3,4,5,6,7,8,9]'
                c = 0.
                d = 1.
                L = d-c
                x = np.array([0,1,2,3,4,5,6,7,8,9],dtype=float)
                n = len(x)
                t = L*np.arange(n)/float(n)

        else:
                print("Which example do you want me to do?")

        plot( t, x, 'o')
        for tj in t: axvline(tj,color='k',alpha=0.25)

        y = fft(x)/np.sqrt(n) # division to get Sauer normalization of DFT

        tt = np.linspace(c,d,1000)
        #plot(tt,myf(tt),'k')
        plot([c,d],[0,0],'k',alpha=0.5)

        p = np.zeros_like(tt,dtype=complex)
        
        assert n%2==0
        kvals = range(-n//2+1,n//2)
            
        #z = y.copy()
        for k in kvals:  
            term = y[k%n]/np.sqrt(n) * np.exp(complex(0.,1.)*2.*np.pi*k*(tt-c)/L)
            p += term
            
        # deal with Nyquist frequency
        for k in [-n//2,n//2]:
            term = y[k%n]/np.sqrt(n) * np.exp(complex(0.,1.)*2.*np.pi*k*(tt-c)/L)
            p += term/2

        plot(tt,myf(tt)   ,'k',label='f',alpha=0.5)
        plot(tt,np.real(p),'r',label='real(Q(t))')
        plot(tt,np.imag(p),'b',label='imag(Q(t))')
        legend()
        title(mytitle)

Spectral differentiation¶

In [35]:
#triginterp.py
import numpy as np
from scipy.fftpack import fft,ifft
from pylab import plot,ion,ioff,draw,title,figure,axvline,legend

for example in [3,2]:

        #figure(example,figsize=(15,10))
        fig, axes = plt.subplots(2,1,figsize=(15,10))
        
        if example==1:
                mytitle = 'vector [1,0,-1,0]'
                c = 0.
                d = 1.
                L = d-c
                x = np.array([1.,0.,-1.,0.])+0.5
                n = len(x)
                t = L*np.arange(n)/float(n)
        elif example==2:
                mytitle = 'distorted sine'
                n = 12#8
                c = 0.
                d = 1.
                L = d-c
                t = np.linspace(0.,L,n,endpoint=False)
                print( len(t) )
                r = 3
                def myf(t): return np.exp(r*np.sin(2*pi*t/L))
                def myfprime(t): return exp(r*sin(2*pi*t/L))*r*cos(2*pi*t/L)*(2*pi/L)
                x = myf(t)
                xp = myfprime(t)
        elif example==3:
                mytitle = 'Runge function'
                n = 8
                c = -1.
                d =  1.
                L = d-c
                t = np.linspace(c,d,n,endpoint=False)
                #print( len(t) )
                def myf(t): return 1/(1+12.*t*t)
                def myfprime(t): return -myf(t)**2 * 24*t
                x = myf(t)
                xp = myfprime(t)

        elif example==4:
                mytitle = 'odd sine'
                n = 8
                c = -1.
                d =  1.
                L = d-c
                t = linspace(c,d,n,endpoint=False)
                #print( len(t) )
                def myf(t): return np.sin(3*pi*(t-c)/L)
                x = myf(t)

        elif example==5:
                mytitle = 'vector [0,1,2,3,4,5,6,7,8,9]'
                c = 0.
                d = 1.
                L = d-c
                x = np.array([0,1,2,3,4,5,6,7,8,9],dtype=float)
                n = len(x)
                t = L*np.arange(n)/float(n)

        else:
                print("Which example do you want me to do?")


        y = fft(x)/sqrt(n) # division to get Sauer def of DFT
        a = real(y)
        b = imag(y)
        print( 'y',y )
        print( 'a',a )
        print( 'b',b )

        tt = np.linspace(c,d,1000)

        p      = np.zeros_like(tt,dtype=complex)  # for interpolant p to f
        pprime = np.zeros_like(tt,dtype=complex)  # for p'

        smart = True #False        
        if not smart:
            kvals = range(0,n)
        else:
            assert n%2==0
            kvals = range(-n//2+1,n//2)
            
        z = y.copy()
        for k in kvals:  
            term = y[k%n]/np.sqrt(n)*np.exp(complex(0.,1.)*2.*np.pi*k*(tt-c)/L)
            p += term
            ################ FORM THE SPECTRAL DERIVATIVE HERE
        
        for k in [-n//2,n//2]:
            term = y[k%n]/np.sqrt(n) * np.exp(complex(0.,1.)*2.*np.pi*k*(tt-c)/L)
            p += term/2

        ax = axes[0]
        ax.plot([c,d],[0,0],'k',alpha=0.5)
        for tj in t: ax.axvline(tj,color='k',alpha=0.25)
        ax.plot( t, x, 'o')
        ax.plot(tt,myf(tt)   ,'k',label='f',alpha=0.5)
        ax.plot(tt,np.real(p),'r',label='real(Q(t))')
        ax.plot(tt,np.imag(p),'b',label='imag(Q(t))')
        ax.legend()
        ax.set_title(mytitle)

        ax = axes[1]
        for tj in t: ax.axvline(tj,color='k',alpha=0.25)
        ax.plot([c,d],[0,0],'k',alpha=0.5)
        ax.plot(tt,myfprime(tt)   ,'c',label="f'",alpha=0.5)
        ################ PLOT THE SPECTRAL DERIVATIVE HERE ######
        ax.legend()
        ax.set_title(mytitle)

y [ 1.0528271 +0.j -0.54755513+0.j  0.20397311-0.j -0.10515882+0.j
  0.0622259 +0.j -0.10515882+0.j  0.20397311+0.j -0.54755513+0.j]
a [ 1.0528271  -0.54755513  0.20397311 -0.10515882  0.0622259  -0.10515882
  0.20397311 -0.54755513]
b [ 0.  0. -0.  0.  0.  0.  0.  0.]
12
y [ 1.69075637e+01+0.00000000e+00j -4.63366963e-16-1.36948853e+01j
 -7.77771148e+00+1.33139821e-15j -5.12790050e-16+3.32514265e+00j
  1.13109459e+00+1.33139821e-15j  1.48894706e-15-3.31440380e-01j
 -1.50228316e-01+0.00000000e+00j  1.48894706e-15+3.31440380e-01j
  1.13109459e+00-1.33139821e-15j -5.12790050e-16-3.32514265e+00j
 -7.77771148e+00-1.33139821e-15j -4.63366963e-16+1.36948853e+01j]
a [ 1.69075637e+01 -4.63366963e-16 -7.77771148e+00 -5.12790050e-16
  1.13109459e+00  1.48894706e-15 -1.50228316e-01  1.48894706e-15
  1.13109459e+00 -5.12790050e-16 -7.77771148e+00 -4.63366963e-16]
b [ 0.00000000e+00 -1.36948853e+01  1.33139821e-15  3.32514265e+00
  1.33139821e-15 -3.31440380e-01  0.00000000e+00  3.31440380e-01
 -1.33139821e-15 -3.32514265e+00 -1.33139821e-15  1.36948853e+01]
In [36]:
#triginterp.py
import numpy as np
from scipy.fftpack import fft,ifft
from pylab import plot,ion,ioff,draw,title,figure,axvline,legend

for example in [3,2,4]:

        #figure(example,figsize=(15,10))
        fig, axes = plt.subplots(2,1,figsize=(15,10))
        
        if example==1:
                mytitle = 'vector [1,0,-1,0]'
                c = 0.
                d = 1.
                L = d-c
                x = np.array([1.,0.,-1.,0.])+0.5
                n = len(x)
                t = L*np.arange(n)/float(n)
        elif example==2:
                mytitle = 'distorted sine'
                n = 12#8
                c = 0.
                d = 1.
                L = d-c
                t = np.linspace(0.,L,n,endpoint=False)
                print( len(t) )
                r = 3
                def myf(t): return np.exp(r*np.sin(2*pi*t/L))
                def myfprime(t): return exp(r*sin(2*pi*t/L))*r*cos(2*pi*t/L)*(2*pi/L)
                x = myf(t)
                xp = myfprime(t)
        elif example==3:
                mytitle = 'Runge function'
                n = 16
                c = -1.
                d =  1.
                L = d-c
                t = np.linspace(c,d,n,endpoint=False)
                print( len(t) )
                def myf(t): return 1/(1+12.*t*t)
                def myfprime(t): return -myf(t)**2 * 24*t
                x = myf(t)
                xp = myfprime(t)

        elif example==4:
                mytitle = 'odd sine'
                n = 8
                c = -1.
                d =  1.
                L = d-c
                t = linspace(c,d,n,endpoint=False)
                print( len(t) )
                def myf(t): return np.sin(3*pi*(t-c)/L)
                def myfprime(t): return 3*pi/L*np.cos(3*pi*(t-c)/L)
                x = myf(t)

        elif example==5:
                mytitle = 'vector [0,1,2,3,4,5,6,7,8,9]'
                c = 0.
                d = 1.
                L = d-c
                x = np.array([0,1,2,3,4,5,6,7,8,9],dtype=float)
                n = len(x)
                t = L*np.arange(n)/float(n)

        else:
                print("Which example do you want me to do?")


        y = fft(x)/sqrt(n) # division to get Sauer def of DFT
        a = real(y)
        b = imag(y)
        print( 'y',y )
        print( 'a',a )
        print( 'b',b )

        tt = np.linspace(c,d,1000)

        p      = np.zeros_like(tt,dtype=complex)  # for interpolant p to f
        pprime = np.zeros_like(tt,dtype=complex)  # for p'

        smart = True #False        
        if not smart:
            kvals = range(0,n)
        else:
            assert n%2==0
            kvals = range(-n//2+1,n//2)
            
        I = complex(0.,1.)
        z = y.copy()
        for k in kvals:  
            term = y[k%n]/np.sqrt(n)*np.exp(I*2.*np.pi*k*(tt-c)/L)
            p += term
            ################ FORM THE SPECTRAL DERIVATIVE HERE
            dterm = y[k%n]/np.sqrt(n)*np.exp(I*2.*np.pi*k*(tt-c)/L) * I*2.*np.pi*k/L
            pprime += dterm
            
        for k in [-n//2,n//2]:
            term = y[k%n]/np.sqrt(n) * np.exp(I*2.*np.pi*k*(tt-c)/L)
            p += term/2

        ax = axes[0]
        ax.plot([c,d],[0,0],'k',alpha=0.5)
        for tj in t: ax.axvline(tj,color='k',alpha=0.25)
        ax.plot( t, x, 'o')
        ax.plot(tt,myf(tt)   ,'k',label='f',alpha=0.5)
        ax.plot(tt,np.real(p),'r',label='real(Q(t))')
        ax.plot(tt,np.imag(p),'b',label='imag(Q(t))')
        ax.legend()
        ax.set_title(mytitle)

        ax = axes[1]
        for tj in t: ax.axvline(tj,color='k',alpha=0.25)
        ax.plot([c,d],[0,0],'k',alpha=0.5)
        ax.plot(tt,myfprime(tt)   ,'c',label="f'",alpha=0.5)
        ################ PLOT THE SPECTRAL DERIVATIVE HERE ######
        ax.plot(tt,np.real(pprime),'r',label="re(Q'(t))")
        ax.plot(tt,np.imag(pprime),'b',label="im(Q'(t))")
        ax.legend()
        ax.set_title(mytitle)

16
y [ 1.48855217e+00+0.j -7.68392493e-01+0.j  2.82946584e-01-0.j
 -1.26041986e-01+0.j  4.40003545e-02+0.j -2.26750424e-02+0.j
  5.51495441e-03+0.j -5.96740156e-03+0.j  3.70198257e-04+0.j
 -5.96740156e-03+0.j  5.51495441e-03-0.j -2.26750424e-02+0.j
  4.40003545e-02-0.j -1.26041986e-01+0.j  2.82946584e-01+0.j
 -7.68392493e-01+0.j]
a [ 1.48855217e+00 -7.68392493e-01  2.82946584e-01 -1.26041986e-01
  4.40003545e-02 -2.26750424e-02  5.51495441e-03 -5.96740156e-03
  3.70198257e-04 -5.96740156e-03  5.51495441e-03 -2.26750424e-02
  4.40003545e-02 -1.26041986e-01  2.82946584e-01 -7.68392493e-01]
b [ 0.  0. -0.  0.  0.  0.  0.  0.  0.  0. -0.  0. -0.  0.  0.  0.]
12
y [ 1.69075637e+01+0.00000000e+00j -4.63366963e-16-1.36948853e+01j
 -7.77771148e+00+1.33139821e-15j -5.12790050e-16+3.32514265e+00j
  1.13109459e+00+1.33139821e-15j  1.48894706e-15-3.31440380e-01j
 -1.50228316e-01+0.00000000e+00j  1.48894706e-15+3.31440380e-01j
  1.13109459e+00-1.33139821e-15j -5.12790050e-16-3.32514265e+00j
 -7.77771148e+00-1.33139821e-15j -4.63366963e-16+1.36948853e+01j]
a [ 1.69075637e+01 -4.63366963e-16 -7.77771148e+00 -5.12790050e-16
  1.13109459e+00  1.48894706e-15 -1.50228316e-01  1.48894706e-15
  1.13109459e+00 -5.12790050e-16 -7.77771148e+00 -4.63366963e-16]
b [ 0.00000000e+00 -1.36948853e+01  1.33139821e-15  3.32514265e+00
  1.33139821e-15 -3.31440380e-01  0.00000000e+00  3.31440380e-01
 -1.33139821e-15 -3.32514265e+00 -1.33139821e-15  1.36948853e+01]
8
y [ 0.52913004+0.00000000e+00j  1.00683487-1.96261557e-16j
 -0.85355339+2.74766180e-16j -0.29972809-3.92523115e-17j
 -0.23623682+0.00000000e+00j -0.29972809+3.92523115e-17j
 -0.85355339-2.74766180e-16j  1.00683487+1.96261557e-16j]
a [ 0.52913004  1.00683487 -0.85355339 -0.29972809 -0.23623682 -0.29972809
 -0.85355339  1.00683487]
b [ 0.00000000e+00 -1.96261557e-16  2.74766180e-16 -3.92523115e-17
  0.00000000e+00  3.92523115e-17 -2.74766180e-16  1.96261557e-16]
In [ ]: