python：正弦波の合成の仕方

Question

素人質問で恐縮です。
周波数を設定した正弦波を合成したいのですが、どうすればいいか分かりません。

a = 1 # 振幅
fs = 10000 # サンプリング周波数
f1 = 5 # 周波数
f2 = 15
f3 = 25
sec = 5 # 秒

swav = []
swav1 = []
swav2 = []
swav3 = []

for n in np.arange(fs * sec):
s = a * np.sin(2.0 * np.pi * f1 * n / fs)
swav1.append(s)

for n in np.arange(fs * sec):
s = a * np.sin(2.0 * np.pi * f2 * n / fs)
swav2.append(s)

for n in np.arange(fs * sec):
s = a * np.sin(2.0 * np.pi * f3 * n / fs)
swav3.append(s)

swav = swav1 + swav2 + swav3※

plt.plot(swav[0:10000])
plt.show()

上記のように周波数を設定し、足し合わせ方が分かりません。

宜しくお願い致します。

Answer 1

他の皆さんが「どう直せばいいか」をコメントしておられますが、「なぜうまくいかないか」について若干コメントします。（もうお分かりと思いますが・・・）

質問者さんのコードの問題は「List同士を+で演算する意味」を勘違いしている点にあったと思います。

[1, 2] + [3, 4]

の結果は

[4, 6]ではなく
[1, 2, 3, 4]なのですね。

listに対する+は加算ではなく連結のため、先頭の10000要素を取り出してプロットしても周波数f1の波形しかでてこないわけです。

ところで配列の要素同士の演算(本件の加算のように)をしているコードをよく見かけます。Pythonには型宣言がないため、自分などは何型を使っているのかピンときてませんでした。しかし今にしてようやくそれらがlistでなくnumpy.ndarrayだからできるのだということがわかってきました。

np.array([1, 2]) + np.array([3, 4])

上記のようにすると結果は

array([4, 6])

になります。質問者さんの期待はこのようになることですよね。

---

listとして扱うなら要素に対する演算、それをするための要素への分解、結果をlistへ再構成することを組みで考える必要があると思います。それがKojiDoiさん、can110さんの回答になるのですね。

一方numpy.ndarrayとして扱うなら配列同士に対して直接演算すれば自動的に要素同士の演算をして元の形で結果を返してくれるという大変便利なことになっていますね。それがmagichanさん、mkgreiさんの回答になるのだと思います。

[...]はPythonではlist型で、np.array([...])とやるとそれはnumpy.ndarray型になりますが、前者を文字通り「リスト」、後者を数学的な「配列」と捉えると、この両者の使い分けを考えやすい気がしました。

[...]がなぜarrayではなくlistと呼ばれるか、自分はここ何件かのQ&Aでようやく腑に落ちた気がします。

Answer 2

そもそも、ループは必要ないのでは？

import numpy as np
import matplotlib.pyplot as plt

a = 1  # 振幅
fs = 10000  # サンプリング周波数
f1 = 5  # 周波数
f2 = 15
f3 = 25
sec = 5  # 秒

t = np.arange(fs * sec)
swav1 = a * np.sin(2.0 * np.pi * f1 * t / fs)
swav2 = a * np.sin(2.0 * np.pi * f2 * t / fs)
swav3 = a * np.sin(2.0 * np.pi * f3 * t / fs)
swav = swav1 + swav2 + swav3

plt.plot(swav[:10000])
plt.show()

Answer 3

１つのforループで足すのが正攻法ですが、もしどうしてもというなら以下でもできます。

swav1,swav2,swav3 = [1,2,3],[4,5,6],[7,8,9]

# swav = swav1 + swav2 + swav3※
swav = [s1+s2+s3 for s1,s2,s3 in zip(swav1,swav2,swav3)]
print(swav) # [12, 15, 18]

Answer 4

forループは一つにまとめられますね。
そのループの中で足し算を行なえばいいだけだと思います。

Answer 5

関数で攻めることもできます。

import numpy as np

a = 1.
fs = 10000
f1 = 5.
f2 = 15.
f3 = 25.
sec = 5.

def f(x):
  return a * (np.sin(2.0 * np.pi * f1 * x / fs) +
              np.sin(2.0 * np.pi * f2 * x / fs) +
              np.sin(2.0 * np.pi * f3 * x / fs))

t = np.arange(fs * sec)
swav = f(t)

print(swav)

Answer 6

Sin関数の結果同士を数学的に加算するのであれば、ほかの回答者の方の言う通りです。しかしながら、対象が音波のような場合、周波数の異なる音波同士が強めあったり、打ち消しあったりするので計算式そのものが異なるものになるようです。
私自身、その分野の専門家ではないので具体的な理論および公式を紹介できなくて申し訳ないのですが、音波を対象にしているのであれば、専門書を確認することをお勧めします。

Comments

[KSwordOfHaste]

おそらく位相のことをおっしゃっている気がします。本件では全ての周波数の位相があっていますね(t=0で位相=0)。位相を加味するなら各周波数のSinの引数に位相項を追加することになりますが、計算式自体は大きく変わらないと思います。

[mkgrei]

https://ja.wikipedia.org/wiki/重ね合わせの原理
成り立つのは重ね合わせの原理であって、数学的に加算した結果、波全体が異なる周期で振動しているようにみえるはずです。観測者は波全体しか見れませんので、強めあったり弱め合ったりしているように見えます。ポイントは「みえる」ということであって、複雑な波形をしていても実は単純なSin関数の加算なのです。
https://ja.wikipedia.org/wiki/フーリエ級数
今のシステムの場合、周期的境界条件が入っておりますので、位相をずらしたとしても全体がある定数分時間軸がシフトしたような形になります。

なおSin波とCos波は合成することができます。これは高校数学の話だと記憶しています。
一般に解析関数として表現する場合には合成波は複雑になります。
http://www.riruraru.com/cfv21/phys/superposition.htm
それでも実は周期的境界条件のもとで許される高調波の合成は結構すっきりした式になるはずです。
http://www.geisya.or.jp/~mwm48961/kou2/d_triangle0.html

いずれにしても今は数値的に関数を表現しようとしているので、解析関数としての複雑性を無視して、単純に加算するだけで期待する結果が得られます。

[KSwordOfHaste]

＞　周期的境界条件が入っておりますので、位相をずらしたとしても全体がある定数分時間軸がシフトしたような形になります

おっしゃるとおりと思います。周波数の関係が倍数にならないような一般ケースの場合は位相によって微妙な変化がある・・・程度の考えでコメントしてました～。

[mkgrei]

なにかもやもやしていたら、固定端の条件も今のケースでは重要でした。
つまり、f(t=0)=f(t=5)=0の条件です。

単なる周期的境界条件だけだと、sin関数では自由端の関数系を表現できませんので。

[KSwordOfHaste]

あ、２つ上の自分のコメントはトンチンカンなこといってました。周波数が倍数関係にある合成波の場合が位相差がきいてくるのでした。初期位相の違いが時間経過により変化しないから。失礼しました。