CODEWITHSUNDEEP
pythonnumpy
David Rinck
David Rinck

Reputation: 6986

Fourier Series in Numpy. Question about Previous Answer

I'm trying to replicate an answer given in a previous thread: How to calculate a Fourier series in Numpy?

import numpy as np
import matplotlib.pyplot as plt
import itertools

def func(x):
    if x >= 1.0 or x <= -1.0:
        return 0
    else:
        return (abs(x) - 1.0)
a = 1.0
b = -1.0
N = 128.
time = np.linspace( a, b, N )
y = (np.fromiter(itertools.imap(func, time), 
                  dtype=time.dtype, count=time.shape[0]))
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(time,y) 
period = 2.
def cn(n):
    c = y*np.exp(-1j*2*n*np.pi*time/period)
    return c.sum()/c.size
def f(x, Nh):
    f = np.array([2*cn(i)*np.exp(1j*2*i*np.pi*x/period) for i in range(1,Nh+1)])
    return f.sum()
y2 = np.array([f(t,10).real for t in time])
ax.plot(time, y2)
plt.show()

I'm getting a solution that's close to the right answer, but shifted. I wasn't sure what I am doing wrong. shifted fourier series

Upvotes: 2

Views: 4700

Answers (3)

chan gimeno
chan gimeno

Reputation: 111

A vectorized version of your code:

import numpy as np
import matplotlib.pyplot as plt
from optparse import OptionParser

def func(x):
    return np.where(np.abs(x) >= 1, 0., np.abs(x) - 1.0)

def cn(x, y, n, period):
    c = y * np.exp(-1j * 2. * np.pi * n * x / period)
    return c.sum()/c.size

def f(x, y, Nh, period):
    rng = np.arange(.5, Nh+.5)
    coeffs = np.array([cn(x,y,i,period) for i in rng])
    f = np.array([2. * coeffs[i] * np.exp(1j*2*i*np.pi*x/period) for i in rng])
    return f.sum(axis=0)

if __name__=='__main__':

    Version = '0.1'
    usage = "usage: %prog [options]"

    parser = OptionParser(usage = usage,version="%prog "+Version)
    parser.add_option("-a", dest='a', type='float', default=1., help="initial time")
    parser.add_option("-b", dest='b', type='float', default=-1., help="end time")
    parser.add_option("-N", "--Nt", dest='N', type='int', default=128, help="number of time steps")
    parser.add_option("-p", "--period", dest='period', type='float', default=2., help="period [time span]")
    parser.add_option("--Nh", dest='Nh', type='int', default=10, help="number of fourier series terms")

    (options, args) = parser.parse_args()

    for key,value in options.__dict__.iteritems():
        exec key + ' = ' + repr(value)

    time = np.linspace( a, b, N )
    y = func(time)
    period = np.abs(a-b)

    y2 = f(time,y,Nh,period).real

    fig = plt.figure()
    ax = fig.add_subplot(1,1,1)
    ax.plot(time, y)
    ax.plot(time, y2)
    plt.show()

Given that the code has been saved with the name "fourier_series.py", you could try:

python fourier_series.py -N 512 --Nh 128

in a normal terminal or:

%run fourier_series.py -N 512 --Nh 128

in the ipython console

Upvotes: 1

Paul
Paul

Reputation: 43620

The error seems to be related to your Riemann sum method (right/middle/left) - indicated by regularfry. Using the middle method gives:

enter image description here

Code:

import numpy as np
import matplotlib.pyplot as plt
import itertools

def func(x):
    if x >= 1.0 or x <= -1.0:
        return 0
    else:
        return (abs(x) - 1.0)

a = 1.0
b = -1.0
N = 128.
time = np.linspace( a, b, N )
y = (np.fromiter(itertools.imap(func, time), 
                  dtype=time.dtype, count=time.shape[0]))

period = 2.
def cn(n):
    c = y*np.exp(-1j*2*n*np.pi*time/period)
    return c.sum()/c.size
def f(x, Nh):
    rng = np.arange(.5, Nh+.5)
    f = np.array([2*cn(i)*np.exp(1j*2*i*np.pi*x/period) for i in rng])
    return f.sum()

y2 = np.array([f(t,10).real for t in time])

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(time, y)
ax.plot(time, y2)
plt.show()

As noted by Sven in another question your use of list comprehensions (and imap) instead of arrays and ufuncs is quite inefficient (should you run into performance issues)

Upvotes: 3

regularfry
regularfry

Reputation: 3268

Looks to me like your DC term has got lost somewhere. I can't check for myself right now, but are you sure the 1 in range(1, Nh+1) in f() is correct?

Upvotes: 2

Related Questions