Pythonのフーリエ級数展開について

Question

Pythonでフーリエ級数展開をしたいと考えています

# -*- coding: utf-8 -*-
from pylab import *
from scipy.integrate import quad

x_min = -pi - 0.5
x_max = pi*2 + 0.5
y_min = -0.25
y_max = 1.25
axis([x_min, x_max, y_min, y_max])
xs = linspace(x_min, x_max, 256)

n_max = 5

def fun(x):
  x = (x + pi) % (pi * 2) - pi
  f = abs(x)

  return f
  """
  if x >= 0:
    return 1
  else:
    return 0
  """

def fourier(fun, n_max):
  a = []
  b = []
  res, err = quad(lambda x:fun(x), -pi, pi)
  a += [res/pi]
  b += [0.0]
  for n in range(1, 1+n_max):
    res, err = quad(lambda x:fun(x)*cos(n*x), -pi, pi)
    a+=[res/pi]
  for n in range(1, 1+n_max):
    res, err = quad(lambda x:fun(x)*sin(n*x), pi, pi)
    b+=[res/pi]

  return a, b

a, b = fourier(fun, n_max)

def fourier_function(x,a,b):
   f = 0.5 * a[0]
   for n in range(1, len(a)):
     f += a[n]*cos(n*x) + b[n]*sin(n*x)
   return f

f_fn = amap(lambda x:fourier_function(x,a,b),xs)


#plot(xs, amap(f, xs), 'b:', lw=1)
plot(xs, f_fn, 'r-', lw=1)
show()

a, bをちゃんと求めるようにして実行できるように組んでみたのですが、やはりうまくいきません。

何度も申し訳ないのですがアドバイスいただけたらありがたいです。
よろしくお願いします。

Answer 1

矩形波ならフーリエ係数を数式で簡単に求められるので、まずはその振幅を正しく計算できていることを確認するのが先でしょうね。a, bは正しい値を持っていますか？

Answer 2

(tachikomaさんご指摘のとおり)fourierに間違いがありますね。

プログラミングの問題というよりはフーリエ級数展開の数式そのものの間違いなので、ヒントだけ申し上げますと間違いは2点あります。

おそらく矩形波のフーリエ級数展開の係数を手計算により求めるような解説をご覧になったのだろうと想像します。

手計算では「偶数項が０になるから奇数項だけ残せばよい」みたいな三角関数の計算上の最適化を工夫したりしますが、例えば次数0～N-1の奇数項だけ計算すればよいといった工夫をそのままプログラミングするならrange(N)ではなくrange(1, N, 2)のようにしないとまずいですよね？

プログラムでは数値計算的に係数を計算していますので、そういう工夫は必ずしも必要ありません。かえって混乱しやすくなると思います。手計算で行うような最適化は考えずに、理論式を単純に当てはめればそういう間違いをすることもなく、すんなりでてくると思います。

2点のうち1点は単純なケアレスミスなのでヒントをいっちゃうと答えになってしまいます。ご自分で見つけてみてください。

Comments

[gymgym]

せっかくアドバイスいただいてやってみたのですが、まだうまく実行できませんでした。
もう少しアドバイスいただけないでしょうか。
よろしくお願いします。

[KSwordOfHaste]

quadは定積分を計算する関数ということはおわかりですよね？

[gymgym]

間違いに気づきました。
お手数をおかけしてすみませんでした。

[KSwordOfHaste]

気づいてしまえば「あ・・・しまった」と思われたことでしょう。
quad(f, pi, pi)だと積分区間の幅が0なのでbの項が全部0になってしまいます。