質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

90.84%

  • Python 3.x

    4821questions

    Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

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

解決済

回答 2

投稿 編集

  • 評価
  • クリップ 0
  • VIEW 273

gymgym

score 59

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をちゃんと求めるようにして実行できるように組んでみたのですが、やはりうまくいきません。

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

  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

回答 2

+2

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

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

checkベストアンサー

+1

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

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

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

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

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

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

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2018/04/27 16:12

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

    キャンセル

  • 2018/04/27 16:25

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

    キャンセル

  • 2018/05/01 12:41

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

    キャンセル

  • 2018/05/01 12:55

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

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 90.84%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

同じタグがついた質問を見る

  • Python 3.x

    4821questions

    Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。