In [19]:

    

%matplotlib inline
from matplotlib import pyplot as plt

import numpy as np

# create 100 linearly-spaced points along the interval [-6,6]
I = np.linspace(-6,6,100)

# compute some functions over those points
C = -I**3
S = np.sin(np.pi * I) / (np.pi * I)

fig, (ax1, ax2) = plt.subplots(ncols=2,figsize=(15,5))

ax1.plot(I, C)
ax2.plot(I, S)


    

    
Out[19]:





[<matplotlib.lines.Line2D at 0x7f263df9d390>]






    

In [20]:

    

from scipy import stats

I = np.linspace(-6,6,100)

# compute the normal PDF
N = stats.norm.pdf(I)

# now integrate it using cumulative sum to produce CDF
cdf = np.cumsum(N)

fig, (ax1, ax2) = plt.subplots(ncols=2,figsize=(15,5))

ax1.plot(I, N)
ax2.plot(I, cdf)


    

    
Out[20]:





[<matplotlib.lines.Line2D at 0x7f263db5a990>]






    

In [29]:

    

from scipy.optimize import minimize

# here's how to use optimization to find a minimum

I = np.linspace(-4,4,200)

# define an objective function
def obj_fn(a):
    return a**2 + np.sin(a*6)

# compute it over the range
F = obj_fn(I)

# find the minimum with initial guess of 0
opt = minimize(obj_fn, 0)

# compute numerical 1st derivative
D = np.diff(F)

fig, (ax1, ax2) = plt.subplots(ncols=2,figsize=(15,5))

ax1.plot(I, F)
ax1.plot(opt.x, opt.fun, marker='+', markersize=20, mew=3)
ax1.text(opt.x, opt.fun-2,'x=%f, y=%f' % (opt.x, opt.fun))
ax2.plot(I[:-1], D)


    

    
Out[29]:





[<matplotlib.lines.Line2D at 0x7f263ccbb990>]






    

In [22]:

    

# compute some sine waves
I = np.linspace(-2*np.pi,2*np.pi,400)

S1 = np.sin(I)
S3 = np.sin(I*3)

fig, (ax1, ax2) = plt.subplots(ncols=2,figsize=(15,5))

ax1.plot(I, S1)
ax2.plot(I, S3)


    

    
Out[22]:





[<matplotlib.lines.Line2D at 0x7f263d2e7210>]






    

In [23]:

    

# now make a square wave out of sine waves

# number of partials
n_partials = 5

# frequency ratios are odd; make this a column vector
freq_rat = np.arange(1,n_partials*2+1,2).reshape(-1,1)
freq_rat


    

    
Out[23]:





array([[1],
       [3],
       [5],
       [7],
       [9]])

In [24]:

    

# now create a phase vector for each partial and each grid point
phase = freq_rat * np.tile(I,(n_partials,1))

# result is n_partials rows, and grid size cols
phase.shape


    

    
Out[24]:





(5, 400)

In [25]:

    

# compute sine of all partials in one operation
S = np.sin(phase)

fig, (ax1, ax2) = plt.subplots(ncols=2,figsize=(15,5))

ax1.plot(I, S[0,:])
ax2.plot(I, S[1,:])


    

    
Out[25]:





[<matplotlib.lines.Line2D at 0x7f263d189890>]






    

In [26]:

    

# for a square wave, the amplitude of the ith partial is 1/i,
# for odd partials
amp = 1. / freq_rat
amp


    

    
Out[26]:





array([[ 1.        ],
       [ 0.33333333],
       [ 0.2       ],
       [ 0.14285714],
       [ 0.11111111]])

In [27]:

    

# now multiply sine for each partial by its amplitude
partials = S * amp

# sum along the row axis, that is the partials axis
sq = np.sum(partials,axis=0)

fig, (ax1, ax2) = plt.subplots(ncols=2,figsize=(15,5))

ax1.bar(freq_rat.ravel(), amp)
ax2.plot(I, sq)


    

    
Out[27]:





[<matplotlib.lines.Line2D at 0x7f263d0665d0>]






    

In [28]:

    

# now generalize to a function with any number of partials
def square_wave(I, n_partials):
    freq_rat = np.arange(1,n_partials*2+1,2).reshape(-1,1)
    phase = freq_rat * np.tile(I,(n_partials,1))
    amp = 1. / freq_rat
    partials = np.sin(phase) * amp
    return np.sum(partials,axis=0)

fig, (ax1, ax2) = plt.subplots(ncols=2,figsize=(15,5))

ax1.plot(I, square_wave(I,7))
ax2.plot(I, square_wave(I,13))


    

    
Out[28]:





[<matplotlib.lines.Line2D at 0x7f263ced57d0>]