Fourier Series in Numpy. Question about Previous Answer

Question

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

score 3 · Accepted Answer · edited May 23 '17 at 10:32

3

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)

edited May 23 '17 at 10:32

Community

1
1

answered May 01 '11 at 01:24

Paul

42,322
15
106
123

Good to know! I had jumped into list comprehensions because y = dft.real * np.cos(xn) + dft.imag * np.sin(xn) doesn't work since dft is a list. I wasn't sure how to vectorize it. – David Rinck May 01 '11 at 01:37
also, the reason I used the imap was because when I tried y = func(time) it gave me the following erorr: 'The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()' – David Rinck May 01 '11 at 13:49
1

@svaha Yes, `func`, `cn`, and `f` all need to be vectorized to get the speedup numpy offers. `func`, written for numpy would look something like this: `def func(x):; x=np.abs(x)-1; x[x>=1]==0; x[x<=-1]==0; return x` – Paul May 01 '11 at 14:18

score 2 · Answer 2 · answered Apr 30 '11 at 16:52

2

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?

answered Apr 30 '11 at 16:52

regularfry

3,248
2
21
27

score 1 · Answer 3 · answered Jul 31 '12 at 18:53

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

Fourier Series in Numpy. Question about Previous Answer

3 Answers3