ブログ

割とコンピュータよりの情報をお届けします。

2019年1月23日

Python numpyでFFTを実行する.周波数軸はどうしたら

numpyでFFTを実行する場合,周波数軸はnumpy.linspaceで設定するのが便利である.
よーく見るとずれてしまって困った話.

まず周波数解析の例を示す.周波数軸を作る部分は
np.linspace(0, サンプリング周波数, データ長さ)
で大体あっているのだが,微妙にずれる.
このずれの分が気になったのでよーく考えると以下の例の場合には
freq = np.linspace(0, fsmp - 1/(t_len + t[1] - t[0]), num=len(t));
ということになる.

ピったしの周波数軸が作れる.

import numpy as np
import matplotlib.pyplot as plt

# 例の信号を作成
t = np.linspace(0.001, 4, 4000);
z = 0.1 + 0.2 * np.sin(t * 10 * 2 * np.pi) + 0.2 * np.sin(t * 33 * 2 * np.pi);

# サンプリング周波数
fsmp = 1 / (t[1] - t[0]);
# 解析時間
t_len = t.max() - t.min();

Fz = np.fft.fft(z) / z.shape[0] * 2; # 折り返すのでパワーが2分の1になっている.
Fz[0] = Fz[0] / 2; # 平均成分は折り返さない.
Fz_abs = np.abs(Fz);

freq = np.linspace(0, fsmp, num=len(t));

fname = 'text.csv';
np.savetxt(fname, np.array([freq, np.real(Fz), np.imag(Fz)]).transpose(), fmt='%.18g', delimiter=',');

plt.figure(1)
plt.subplot(311)
plt.plot(t, z)

plt.subplot(312)
plt.plot(freq, Fz_abs)
plt.xlim([0, fsmp/2]);

#136/4
plt.subplot(313)
plt.plot(range(len(t)), Fz_abs)
#plt.xlim([0, 500]); # 後半(高周波側)は前半の折り返し

plt.show();

≫ 続きを読む

2019/01/23 コンピュータ   TakeMe
タグ:Python

Pythonでは複素数の演算が行えることを知った

最近の言語は複素数が使えることが多いが,Pythonも複素数演算が扱える.

例えば以下のような演算ができる.
以下の例では虚部と実部を分けて指定する方法と極座標表示から得る方法の例.

import cmath

# 虚部と実部を分けて指定
a1 = complex(1, 0);
a2 = complex(2, 0);

# 偏角を使った指定
mag1 = 1.0;
arg1 = 0.0;
a1 = complex(mag1 * cmath.cos(arg1), mag1 * cmath.sin(arg1))

mag2 = 2.0;
arg2 = 0.0;
a2 = complex(mag2 * cmath.cos(arg2), mag2 * cmath.sin(arg2))


trial = complex(100, 0);

a21 = a2 - a1;
b = a21 / trial

print(b)

いったん複素数になってしまえば掛け算,割り算は可能.
(例の偏角はラジアンなので注意)

≫ 続きを読む

2019/01/23 コンピュータ   TakeMe
タグ:Python