In [1]:

    

%matplotlib inline
from matplotlib import pyplot as plt

import numpy as np

n=5
I = np.linspace(-2*np.pi,2*np.pi,n)
odd = np.arange(1,14,2)
h, J = np.meshgrid(odd, I)
h, J


    

    
Out[1]:





(array([[ 1,  3,  5,  7,  9, 11, 13],
        [ 1,  3,  5,  7,  9, 11, 13],
        [ 1,  3,  5,  7,  9, 11, 13],
        [ 1,  3,  5,  7,  9, 11, 13],
        [ 1,  3,  5,  7,  9, 11, 13]]),
 array([[-6.28318531, -6.28318531, -6.28318531, -6.28318531, -6.28318531,
         -6.28318531, -6.28318531],
        [-3.14159265, -3.14159265, -3.14159265, -3.14159265, -3.14159265,
         -3.14159265, -3.14159265],
        [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
          0.        ,  0.        ],
        [ 3.14159265,  3.14159265,  3.14159265,  3.14159265,  3.14159265,
          3.14159265,  3.14159265],
        [ 6.28318531,  6.28318531,  6.28318531,  6.28318531,  6.28318531,
          6.28318531,  6.28318531]]))

In [2]:

    

amp = 1. / odd
plt.bar(odd, amp)
odd


    

    
Out[2]:





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






    

In [3]:

    

n=500
I = np.linspace(-2*np.pi,2*np.pi,n)
h, J = np.meshgrid(odd, I)
harmonics = np.sin(J * h) * amp
sq = np.sum(harmonics, axis=1)

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

for i in range(odd.size):
    ax1.plot(I, harmonics[:,i])
    
ax2.plot(I, sq)


    

    
Out[3]:





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






    

In [7]:

    

def square_wave(I,n_harmonics=7):
    odd = np.arange(1,n_harmonics*2,2)
    h, J = np.meshgrid(odd, I)
    amp = 1. / odd
    harmonics = np.sin(J * h) * amp
    return np.sum(harmonics, axis=1)

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

I = np.linspace(-2*np.pi,2*np.pi,n)

ax1.plot(I, square_wave(I,3))
ax2.plot(I, square_wave(I,20))


    

    
Out[7]:





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