pythonでのフーリエ変換の結果について

Question

前提・実現したいこと

FFTを使わずにノイズの乗った三角関数のフーリエ変換を行い結果をプロットしています。
該当のコードのように変換を行いfr、fiをそれぞれグラフにしているのですが、
①乱数の生成にsin波の周波数は関係ないはずなのに
低周波数：フーリエ変換後の実部にノイズの影響が出やすい
高周波数：フーリエ変換後の虚部にノイズの影響が出やすい
となるのはなぜか。
②sin波をcos波に変えると①の結果の逆になるのはなぜか
上の２つについて教えていただきたいです。よろしくお願いいたします。

該当のソースコード

f = 50 #or450　周波数
N=1000

x = np.linspace(0, 1, N)
noise = np.random.rand(N)
y = np.cos( 2 * math.pi * f * x )
y1=y+noise-0.5
fr=np.zeros(N)
fi=np.zeros(N)

for i in range(N):
  for j in range(N):
    fr[i]+=y1[j]*np.cos(-2*math.pi*i*j/N) #y1のフーリエ変換後の実数部分
    fi[i]+=y1[j]*np.sin(-2*math.pi*i*j/N) #y1のフーリエ変換後の虚数部分

Answer 1

実はノイズの量は同じで、グラフの縦軸の表示範囲が違うからノイズの量が違って「見えてる」だけ、というのは、jeanbiegoさんの回答と同じです
グラフの縦軸の表示範囲を全部揃えてグラフを描けば、そのことはすぐに分かります

さらに、コードを下記のように変更してから、実行・グラフ表示させてみてください

x = np.linspace(0, 1, N)

↓ 変更

x = np.arange(N) / N

そうすれば、元信号yがsinの場合は周波数によらずフーリエ変換後は虚部のみに信号の振幅が現れ、元信号yがcosの場合は周波数によらずフーリエ変換後は実部のみに信号の振幅が現れることが分かります

また、上記コード変更をした場合に、あえてグラフの縦軸の表示範囲を指定せずにグラフを描いたら、どうなるのかもやってみてください

Comments

[ryu-sei]

ご回答いただきありがとうございました。
提案してくださったコードも試してみます！

[jbpb0]

詳細はjeanbiegoさんの回答のコメントに書きましたけど、私の回答のようにコードを書いていたら、元信号の周波数を変えてもグラフは(だいたい)同じになるので、グラフの見た目に惑わされずに、実はノイズの量は同じということに(ここで質問するまでもなく)気づいてたかもしれません

実際にコードを変更して実行して表示されるグラフを見たら、私が言ってることが分かると思うので、やってみてください

Answer 2

①乱数の生成にsin波の周波数は関係ないはずなのに

低周波数：フーリエ変換後の実部にノイズの影響が出やすい
高周波数：フーリエ変換後の虚部にノイズの影響が出やすい
となるのはなぜか。

縦軸がバラバラなせいでノイズの周波数成分に大小あるようにみえていますが、実際はどれも同じamplitudeです。つまり、ノイズは実は一切関係ありません。

②sin波をcos波に変えると①の結果の逆になるのはなぜか

そしてなぜ縦軸がバラバラなのかというと、正弦波の変換結果のamplitudeに大きく差が出ているからなので、つまりはΣcos(2πfj/N)×cos(-2πfij/N)とΣcos(2πfj/N)×sin(2πfij/N)の差ということなのですが…
プログラミングの質問というより数式処理の話なので、ご自身で手を動かしてみたほうが早いかもしれません。（すみませんが、私も数式処理は上手くないので他の方にお任せします）

Comments

[ryu-sei]

ご回答いただきありがとうございました。
計算したらわかってきました。

[jbpb0]

jeanbiegoさん

質問に掲載のコードでの、
・「x = np.linspace(0, 1, N)」の「x」
・「for j in range(N):」の「j」を「N」で割ったもの
は、実際に値を見たら分かりますが違うので、

fr[i]+=y1[j]*np.cos(-2*math.pi*i*j/N) #y1のフーリエ変換後の実数部分
fi[i]+=y1[j]*np.sin(-2*math.pi*i*j/N) #y1のフーリエ変換後の虚数部分

は

> Σcos(2πfj/N)×cos(-2πfij/N)とΣcos(2πfj/N)×sin(2πfij/N)

とはならないと思います

たとえば、

N=1000
f = 50
x = np.linspace(0, 1, N)
y = np.cos( 2 * math.pi * f * x )

の場合に、

i = f
j = 800
print(y[j])
print(np.cos(-2*math.pi*i*j/N))

を実行して、二つの「print」の結果を比べたら違いますよね

私の回答に書いてるように、

x = np.linspace(0, 1, N)
↓ 変更
x = np.arange(N) / N

とすれば、二つの「print」の結果が(f, jによらず)一致するので、その場合はjeanbiegoさんが回答に書かれてる数式の通りになると思いますが

[jeanbiego]

>jbpb0さん
ご指摘ありがとうございます。
そういえば、range(1000)/1000は
0 to 0.999で、
x = np.linspace(0, 1, 1000)は
0 to 1だから、
途中の要素も同じ値なわけはありませんね。

[jbpb0]

はい
そこがずれてるため、本来は0.0になるはずの「Sin x Cos の積分」が値を持っていて、下記に元信号の振幅の一部が現れてます
・元信号が「sin」の場合の実部
・元信号が「cos」の場合の虚部

それがどれくらいになるのかが周波数で違います
ずれの影響は、低周波は小さく高周波は大きいためです

ずれないようにコードを書けば、元信号の周波数によらず結果は変わりません
質問に掲載されてるグラフのように、元信号が低周波と高周波の場合で様子が変わるのは、ずれてしまうようにコードを書いたためなので、

> プログラミングの質問というより数式処理の話

とは言い切れないと思います

[jeanbiego]

> そこがずれてるため、本来は積分したら0.0になるはずの「Sin x Cos」が値を持っていて
おお、なるほど！　
cos(θ)sin(iθ)、θ=-2π j/Nとして、jを0からNまでとるというのは、θを0から2πで積分したときの値に近いんですね。iが整数だから、基本は0になるはずと。
わかりやすいです。

>> プログラミングの質問というより数式処理の話
>とは言い切れないと思います
確かにおっしゃる通りで、私の早とちりでした。

[jbpb0]

ああ、そうですね
「f」が整数では無い場合も考えたら、-∞〜∞の範囲で積分しないと「Sin x Cos の積分」は0.0にはなりませんね

「f」が整数の場合に限れば、データが存在する範囲内だけでの積分(実際はN個の総和)が0.0になります
波長が N/f なので、N個のデータがSinとCosの途中で終わることが無いためです