"""Fourier Series"""
from __future__ import print_function, division
from sympy import pi, oo, Wild, Basic
from sympy.core.expr import Expr
from sympy.core.add import Add
from sympy.core.compatibility import is_sequence
from sympy.core.containers import Tuple
from sympy.core.singleton import S
from sympy.core.symbol import Dummy, Symbol
from sympy.core.sympify import sympify
from sympy.functions.elementary.trigonometric import sin, cos, sinc
from sympy.series.series_class import SeriesBase
from sympy.series.sequences import SeqFormula
from sympy.sets.sets import Interval
from sympy.simplify.fu import TR8, TR2, TR1, TR10, sincos_to_sum
def fourier_cos_seq(func, limits, n):
"""Returns the cos sequence in a Fourier series"""
from sympy.integrals import integrate
x, L = limits[0], limits[2] - limits[1]
cos_term = cos(2*n*pi*x / L)
formula = 2 * cos_term * integrate(func * cos_term, limits) / L
a0 = formula.subs(n, S.Zero) / 2
return a0, SeqFormula(2 * cos_term * integrate(func * cos_term, limits)
/ L, (n, 1, oo))
def fourier_sin_seq(func, limits, n):
"""Returns the sin sequence in a Fourier series"""
from sympy.integrals import integrate
x, L = limits[0], limits[2] - limits[1]
sin_term = sin(2*n*pi*x / L)
return SeqFormula(2 * sin_term * integrate(func * sin_term, limits)
/ L, (n, 1, oo))
def _process_limits(func, limits):
"""
Limits should be of the form (x, start, stop).
x should be a symbol. Both start and stop should be bounded.
* If x is not given, x is determined from func.
* If limits is None. Limit of the form (x, -pi, pi) is returned.
Examples
========
>>> from sympy import pi
>>> from sympy.series.fourier import _process_limits as pari
>>> from sympy.abc import x
>>> pari(x**2, (x, -2, 2))
(x, -2, 2)
>>> pari(x**2, (-2, 2))
(x, -2, 2)
>>> pari(x**2, None)
(x, -pi, pi)
"""
def _find_x(func):
free = func.free_symbols
if len(free) == 1:
return free.pop()
elif not free:
return Dummy('k')
else:
raise ValueError(
" specify dummy variables for %s. If the function contains"
" more than one free symbol, a dummy variable should be"
" supplied explicitly e.g. FourierSeries(m*n**2, (n, -pi, pi))"
% func)
x, start, stop = None, None, None
if limits is None:
x, start, stop = _find_x(func), -pi, pi
if is_sequence(limits, Tuple):
if len(limits) == 3:
x, start, stop = limits
elif len(limits) == 2:
x = _find_x(func)
start, stop = limits
if not isinstance(x, Symbol) or start is None or stop is None:
raise ValueError('Invalid limits given: %s' % str(limits))
unbounded = [S.NegativeInfinity, S.Infinity]
if start in unbounded or stop in unbounded:
raise ValueError("Both the start and end value should be bounded")
return sympify((x, start, stop))
def finite_check(f, x, L):
def check_fx(exprs, x):
return x not in exprs.free_symbols
def check_sincos(_expr, x, L):
if isinstance(_expr, (sin, cos)):
sincos_args = _expr.args[0]
if sincos_args.match(a*(pi/L)*x + b) is not None:
return True
else:
return False
_expr = sincos_to_sum(TR2(TR1(f)))
add_coeff = _expr.as_coeff_add()
a = Wild('a', properties=[lambda k: k.is_Integer, lambda k: k != S.Zero, ])
b = Wild('b', properties=[lambda k: x not in k.free_symbols, ])
for s in add_coeff[1]:
mul_coeffs = s.as_coeff_mul()[1]
for t in mul_coeffs:
if not (check_fx(t, x) or check_sincos(t, x, L)):
return False, f
return True, _expr
[docs]class FourierSeries(SeriesBase):
r"""Represents Fourier sine/cosine series.
This class only represents a fourier series.
No computation is performed.
For how to compute Fourier series, see the :func:`fourier_series`
docstring.
See Also
========
sympy.series.fourier.fourier_series
"""
def __new__(cls, *args):
args = map(sympify, args)
return Expr.__new__(cls, *args)
@property
def function(self):
return self.args[0]
@property
def x(self):
return self.args[1][0]
@property
def period(self):
return (self.args[1][1], self.args[1][2])
@property
def a0(self):
return self.args[2][0]
@property
def an(self):
return self.args[2][1]
@property
def bn(self):
return self.args[2][2]
@property
def interval(self):
return Interval(0, oo)
@property
def start(self):
return self.interval.inf
@property
def stop(self):
return self.interval.sup
@property
def length(self):
return oo
@property
def L(self):
return abs(self.period[1] - self.period[0]) / 2
def _eval_subs(self, old, new):
x = self.x
if old.has(x):
return self
[docs] def truncate(self, n=3):
"""
Return the first n nonzero terms of the series.
If n is None return an iterator.
Parameters
==========
n : int or None
Amount of non-zero terms in approximation or None.
Returns
=======
Expr or iterator
Approximation of function expanded into Fourier series.
Examples
========
>>> from sympy import fourier_series, pi
>>> from sympy.abc import x
>>> s = fourier_series(x, (x, -pi, pi))
>>> s.truncate(4)
2*sin(x) - sin(2*x) + 2*sin(3*x)/3 - sin(4*x)/2
See Also
========
sympy.series.fourier.FourierSeries.sigma_approximation
"""
if n is None:
return iter(self)
terms = []
for t in self:
if len(terms) == n:
break
if t is not S.Zero:
terms.append(t)
return Add(*terms)
[docs] def sigma_approximation(self, n=3):
r"""
Return :math:`\sigma`-approximation of Fourier series with respect
to order n.
Sigma approximation adjusts a Fourier summation to eliminate the Gibbs
phenomenon which would otherwise occur at discontinuities.
A sigma-approximated summation for a Fourier series of a T-periodical
function can be written as
.. math::
s(\theta) = \frac{1}{2} a_0 + \sum _{k=1}^{m-1}
\operatorname{sinc} \Bigl( \frac{k}{m} \Bigr) \cdot
\left[ a_k \cos \Bigl( \frac{2\pi k}{T} \theta \Bigr)
+ b_k \sin \Bigl( \frac{2\pi k}{T} \theta \Bigr) \right],
where :math:`a_0, a_k, b_k, k=1,\ldots,{m-1}` are standard Fourier
series coefficients and
:math:`\operatorname{sinc} \Bigl( \frac{k}{m} \Bigr)` is a Lanczos
:math:`\sigma` factor (expressed in terms of normalized
:math:`\operatorname{sinc}` function).
Parameters
==========
n : int
Highest order of the terms taken into account in approximation.
Returns
=======
Expr
Sigma approximation of function expanded into Fourier series.
Examples
========
>>> from sympy import fourier_series, pi
>>> from sympy.abc import x
>>> s = fourier_series(x, (x, -pi, pi))
>>> s.sigma_approximation(4)
2*sin(x)*sinc(pi/4) - 2*sin(2*x)/pi + 2*sin(3*x)*sinc(3*pi/4)/3
See Also
========
sympy.series.fourier.FourierSeries.truncate
Notes
=====
The behaviour of
:meth:`~sympy.series.fourier.FourierSeries.sigma_approximation`
is different from :meth:`~sympy.series.fourier.FourierSeries.truncate`
- it takes all nonzero terms of degree smaller than n, rather than
first n nonzero ones.
References
==========
.. [1] https://en.wikipedia.org/wiki/Gibbs_phenomenon
.. [2] https://en.wikipedia.org/wiki/Sigma_approximation
"""
terms = [sinc(pi * i / n) * t for i, t in enumerate(self[:n])
if t is not S.Zero]
return Add(*terms)
[docs] def shift(self, s):
"""Shift the function by a term independent of x.
f(x) -> f(x) + s
This is fast, if Fourier series of f(x) is already
computed.
Examples
========
>>> from sympy import fourier_series, pi
>>> from sympy.abc import x
>>> s = fourier_series(x**2, (x, -pi, pi))
>>> s.shift(1).truncate()
-4*cos(x) + cos(2*x) + 1 + pi**2/3
"""
s, x = sympify(s), self.x
if x in s.free_symbols:
raise ValueError("'%s' should be independent of %s" % (s, x))
a0 = self.a0 + s
sfunc = self.function + s
return self.func(sfunc, self.args[1], (a0, self.an, self.bn))
[docs] def shiftx(self, s):
"""Shift x by a term independent of x.
f(x) -> f(x + s)
This is fast, if Fourier series of f(x) is already
computed.
Examples
========
>>> from sympy import fourier_series, pi
>>> from sympy.abc import x
>>> s = fourier_series(x**2, (x, -pi, pi))
>>> s.shiftx(1).truncate()
-4*cos(x + 1) + cos(2*x + 2) + pi**2/3
"""
s, x = sympify(s), self.x
if x in s.free_symbols:
raise ValueError("'%s' should be independent of %s" % (s, x))
an = self.an.subs(x, x + s)
bn = self.bn.subs(x, x + s)
sfunc = self.function.subs(x, x + s)
return self.func(sfunc, self.args[1], (self.a0, an, bn))
[docs] def scale(self, s):
"""Scale the function by a term independent of x.
f(x) -> s * f(x)
This is fast, if Fourier series of f(x) is already
computed.
Examples
========
>>> from sympy import fourier_series, pi
>>> from sympy.abc import x
>>> s = fourier_series(x**2, (x, -pi, pi))
>>> s.scale(2).truncate()
-8*cos(x) + 2*cos(2*x) + 2*pi**2/3
"""
s, x = sympify(s), self.x
if x in s.free_symbols:
raise ValueError("'%s' should be independent of %s" % (s, x))
an = self.an.coeff_mul(s)
bn = self.bn.coeff_mul(s)
a0 = self.a0 * s
sfunc = self.args[0] * s
return self.func(sfunc, self.args[1], (a0, an, bn))
[docs] def scalex(self, s):
"""Scale x by a term independent of x.
f(x) -> f(s*x)
This is fast, if Fourier series of f(x) is already
computed.
Examples
========
>>> from sympy import fourier_series, pi
>>> from sympy.abc import x
>>> s = fourier_series(x**2, (x, -pi, pi))
>>> s.scalex(2).truncate()
-4*cos(2*x) + cos(4*x) + pi**2/3
"""
s, x = sympify(s), self.x
if x in s.free_symbols:
raise ValueError("'%s' should be independent of %s" % (s, x))
an = self.an.subs(x, x * s)
bn = self.bn.subs(x, x * s)
sfunc = self.function.subs(x, x * s)
return self.func(sfunc, self.args[1], (self.a0, an, bn))
def _eval_as_leading_term(self, x):
for t in self:
if t is not S.Zero:
return t
def _eval_term(self, pt):
if pt == 0:
return self.a0
return self.an.coeff(pt) + self.bn.coeff(pt)
def __neg__(self):
return self.scale(-1)
def __add__(self, other):
if isinstance(other, FourierSeries):
if self.period != other.period:
raise ValueError("Both the series should have same periods")
x, y = self.x, other.x
function = self.function + other.function.subs(y, x)
if self.x not in function.free_symbols:
return function
an = self.an + other.an
bn = self.bn + other.bn
a0 = self.a0 + other.a0
return self.func(function, self.args[1], (a0, an, bn))
return Add(self, other)
def __sub__(self, other):
return self.__add__(-other)
class FiniteFourierSeries(FourierSeries):
r"""Represents Finite Fourier sine/cosine series.
For how to compute Fourier series, see the :func:`fourier_series`
docstring.
Parameters
==========
f : Expr
Expression for finding fourier_series
limits : ( x, start, stop)
x is the independent variable for the expression f
(start, stop) is the period of the fourier series
exprs: (a0, an, bn) or Expr
a0 is the constant term a0 of the fourier series
an is a dictionary of coefficients of cos terms
an[k] = coefficient of cos(pi*(k/L)*x)
bn is a dictionary of coefficients of sin terms
bn[k] = coefficient of sin(pi*(k/L)*x)
or exprs can be an expression to be converted to fourier form
Methods
=======
This class is an extension of FourierSeries class.
Please refer to sympy.series.fourier.FourierSeries for
further information.
See Also
========
sympy.series.fourier.FourierSeries
sympy.series.fourier.fourier_series
"""
def __new__(cls, f, limits, exprs):
if not (type(exprs) == tuple and len(exprs) == 3): # exprs is not of form (a0, an, bn)
# Converts the expression to fourier form
c, e = exprs.as_coeff_add()
rexpr = c + Add(*[TR10(i) for i in e])
a0, exp_ls = rexpr.expand(trig=False, power_base=False, power_exp=False, log=False).as_coeff_add()
x = limits[0]
L = abs(limits[2] - limits[1]) / 2
a = Wild('a', properties=[lambda k: k.is_Integer, lambda k: k is not S.Zero, ])
b = Wild('b', properties=[lambda k: x not in k.free_symbols, ])
an = dict()
bn = dict()
# separates the coefficients of sin and cos terms in dictionaries an, and bn
for p in exp_ls:
t = p.match(b * cos(a * (pi / L) * x))
q = p.match(b * sin(a * (pi / L) * x))
if t:
an[t[a]] = t[b] + an.get(t[a], S.Zero)
elif q:
bn[q[a]] = q[b] + bn.get(q[a], S.Zero)
else:
a0 += p
exprs = (a0, an, bn)
args = map(sympify, (f, limits, exprs))
return Expr.__new__(cls, *args)
@property
def interval(self):
_length = 1 if self.a0 else 0
_length += max(set(self.an.keys()).union(set(self.bn.keys()))) + 1
return Interval(0, _length)
@property
def length(self):
return self.stop - self.start
def shiftx(self, s):
s, x = sympify(s), self.x
if x in s.free_symbols:
raise ValueError("'%s' should be independent of %s" % (s, x))
_expr = self.truncate().subs(x, x + s)
sfunc = self.function.subs(x, x + s)
return self.func(sfunc, self.args[1], _expr)
def scale(self, s):
s, x = sympify(s), self.x
if x in s.free_symbols:
raise ValueError("'%s' should be independent of %s" % (s, x))
_expr = self.truncate() * s
sfunc = self.function * s
return self.func(sfunc, self.args[1], _expr)
def scalex(self, s):
s, x = sympify(s), self.x
if x in s.free_symbols:
raise ValueError("'%s' should be independent of %s" % (s, x))
_expr = self.truncate().subs(x, x * s)
sfunc = self.function.subs(x, x * s)
return self.func(sfunc, self.args[1], _expr)
def _eval_term(self, pt):
if pt == 0:
return self.a0
_term = self.an.get(pt, S.Zero) * cos(pt * (pi / self.L) * self.x) \
+ self.bn.get(pt, S.Zero) * sin(pt * (pi / self.L) * self.x)
return _term
def __add__(self, other):
if isinstance(other, FourierSeries):
return other.__add__(fourier_series(self.function, self.args[1],\
finite=False))
elif isinstance(other, FiniteFourierSeries):
if self.period != other.period:
raise ValueError("Both the series should have same periods")
x, y = self.x, other.x
function = self.function + other.function.subs(y, x)
if self.x not in function.free_symbols:
return function
return fourier_series(function, limits=self.args[1])
[docs]def fourier_series(f, limits=None, finite=True):
"""Computes Fourier sine/cosine series expansion.
Returns a :class:`FourierSeries` object.
Examples
========
>>> from sympy import fourier_series, pi, cos
>>> from sympy.abc import x
>>> s = fourier_series(x**2, (x, -pi, pi))
>>> s.truncate(n=3)
-4*cos(x) + cos(2*x) + pi**2/3
Shifting
>>> s.shift(1).truncate()
-4*cos(x) + cos(2*x) + 1 + pi**2/3
>>> s.shiftx(1).truncate()
-4*cos(x + 1) + cos(2*x + 2) + pi**2/3
Scaling
>>> s.scale(2).truncate()
-8*cos(x) + 2*cos(2*x) + 2*pi**2/3
>>> s.scalex(2).truncate()
-4*cos(2*x) + cos(4*x) + pi**2/3
Notes
=====
Computing Fourier series can be slow
due to the integration required in computing
an, bn.
It is faster to compute Fourier series of a function
by using shifting and scaling on an already
computed Fourier series rather than computing
again.
e.g. If the Fourier series of ``x**2`` is known
the Fourier series of ``x**2 - 1`` can be found by shifting by ``-1``.
See Also
========
sympy.series.fourier.FourierSeries
References
==========
.. [1] mathworld.wolfram.com/FourierSeries.html
"""
f = sympify(f)
limits = _process_limits(f, limits)
x = limits[0]
if x not in f.free_symbols:
return f
if finite:
L = abs(limits[2] - limits[1]) / 2
is_finite, res_f = finite_check(f, x, L)
if is_finite:
return FiniteFourierSeries(f, limits, res_f)
n = Dummy('n')
neg_f = f.subs(x, -x)
if f == neg_f:
a0, an = fourier_cos_seq(f, limits, n)
bn = SeqFormula(0, (1, oo))
elif f == -neg_f:
a0 = S.Zero
an = SeqFormula(0, (1, oo))
bn = fourier_sin_seq(f, limits, n)
else:
a0, an = fourier_cos_seq(f, limits, n)
bn = fourier_sin_seq(f, limits, n)
return FourierSeries(f, limits, (a0, an, bn))
