インバーター作るぞー！とか言いつつ全く進んでいませんが，気が向いたので先に波形表示するプログラムを書いてしまいました。
こんな感じ
このプログラムを使って制御するわけではないですが，単純な正弦波を見ているよりは何倍も面白いので作って正解だったと思います。
それと，任意の波形を表示するだけだとつまらないので，リアルタイム表示するようにしました。計算量が多いのでヌルヌルとは動きませんけど。
参考にしたサイト
今回書いたプログラム
- Calculation_wave.py
- plot_wave.py
１で波形表示に必要な計算をして，２で表示しています。
１．Calculation_wave.py
import numpy as np
import matplotlib.pyplot as plt
class Inverter_wave:
def __init__(self, E, h=2**-7, N=40):
self.E = E # インバーターへの入力
self.h = h # 計算の刻み幅
self.N = N # 三角波の定数
self.pi = np.pi
self.x_num = int(2*self.pi / h) # x方向のの分割数
# 三角変調波（正弦波）を計算
def sin_wave(self, F, Es=200, P_d=0):
"""
F : 変調波信号の周波数[Hz]
Es : 変調波の振幅[V]
P_d : 位相差(Phase difference)[rad]
"""
self.F = F
self.Es = Es
x_list = np.linspace(0, 2*self.pi, self.x_num)
y_list = Es * np.sin((x_list * F) + P_d)
return (x_list, y_list)
# 三角搬送波(triangular wave)を計算
def tri_wave(self, P, Ec=200, P_d=0):
"""
P : キャリア周波数[Hz]
Ec: 搬送波（キャリア）信号の振幅[V]
P_d : 位相差
"""
self.P = P
self.Ec = Ec
def f(x, P, N=self.N):
y = 0
for n in range(1, N+1):
y += ((np.sin(n*self.pi/2) * np.sin(n*(x + P_d)*P)) / n**2) * (8/self.pi) / 3
return y
x_list = np.linspace(0, 2*self.pi, self.x_num)
y_list = []
for s,t in enumerate(x_list):
y_list.append(f(t, P) * Ec)
return (x_list, y_list)
# 矩形波(rectangle wave)の座標を計算
def rectangle_wave(self, sin_x, sin_y, tri_x, tri_y, V=200):
# sin_y > tri_y のとき V/2
# sin_y <= tri_y のとき -V/2
x_list = sin_x
y_list = []
for i in range(self.x_num):
if sin_y[i] > tri_y[i]:
y_list.append(V/2)
else:
y_list.append(-V/2)
return (x_list, y_list)
def plot_info(self):
# 変調率(modulation factor)
M = self.Es / self.Ec
# 相電圧
Vph = M * self.E / 2
# 線間電圧
V1rms = np.sqrt(3) * Vph
print("インバーターへの入力電圧：%d[V]" % self.E)
if (1 < M):
print("変調率 ：{:.2f} (過変調PWM)".format(M))
else:
print("変調率 ：{:.2f}".format(M))
print("相電圧 ：{:.2f}[V]".format(Vph))
print("線間電圧：{:.2f}[V]".format(V1rms))
print('変調波周波数：{:.2f}[Hz]'.format(self.F))
print('キャリア周波数：{:.2f}[Hz]'.format(self.P))
# example of use
if __name__ == '__main__':
plt.title('PWM')
plt.ylabel('Voltage')
plt.xlabel('Range : 0-2π[rad]')
plt.xlim(0, 2*np.pi)
wave = Inverter_wave(E=1)
sin_x, sin_y = wave.sin_wave(F=2, Es=0.5)
tri_x, tri_y = wave.tri_wave(P=12, Ec=1)
rec_x, rec_y = wave.rectangle_wave(sin_x, sin_y, tri_x, tri_y)
plt.plot(sin_x, sin_y, color='green')
plt.plot(tri_x, tri_y, color='red')
plt.plot(rec_x, rec_y, color='blue')
plt.show()
# 電圧などの情報を表示
wave.plot_info()
それぞれの関数は，0〜2πの範囲を計算して，要素数を揃えたNumpy配列の形で返しています。
式中に位相を含めてあるので，波形が進んでいく様子を見ることもできます。
※追記(2018/6/3)：numpyを活用して上記のプログラムを改善しました．このページの一番下に改善したものを載せておきます．
２．plot_wave.py
import numpy as np
import matplotlib.pyplot as plt
from Calculation_wave import *
pi = np.pi
def wave_plot(E=200, F=1, Es=160, P=18, Ec=200, dif=0):
"""
E ： インバーターへの入力直流電圧
F : 三角変調波の周波数
Es: 三角変調波の振幅
P : 三角搬送波の周波数
Ec: 三角搬送波の振幅
dif : 三角変調波の位相
"""
wave = Inverter_wave(E)
# 変調率(Modulation factor)
M = Es / Ec
# 相電圧
Vph = M * E / 2
# 線間電圧
V1rms = np.sqrt(3) * Vph
# 三角変調波
# U相
sin1_x, sin1_y = wave.sin_wave(F, Es, P_d= 0 + dif)
# V相
sin2_x, sin2_y = wave.sin_wave(F, Es, P_d= 2/3*pi + dif)
# W相
sin3_x, sin3_y = wave.sin_wave(F, Es, P_d= 4/3*pi + dif)
# 三角搬送波
tri_x, tri_y = wave.tri_wave(P, Ec, P_d=dif)
# 矩形波
rec1_x, rec1_y = wave.rectangle_wave(sin1_x, sin1_y, tri_x, tri_y, Vph)
rec2_x, rec2_y = wave.rectangle_wave(sin2_x, sin2_y, tri_x, tri_y, Vph)
rec3_x, rec3_y = wave.rectangle_wave(sin3_x ,sin3_y, tri_x, tri_y, Vph)
# 軸線を表示
x_axis = np.linspace(0, 2*pi, 2)
y_axis = np.zeros(2)
# U-V線間電圧
rec_U_V = (np.array(rec1_y) - np.array(rec2_y)) / Vph * V1rms
# V-W線間電圧
rec_V_W = (np.array(rec2_y) - np.array(rec3_y)) / Vph * V1rms
# W-U線間電圧
rec_W_U = (np.array(rec3_y) - np.array(rec1_y)) / Vph * V1rms
# ---------- 各相のスイッチング波形を表示 ---------- #
ax1 = plt.subplot(421) # 合計4行2列,１番目
ax1.set_title('PWM control')
ax1.plot(sin1_x, sin1_y, color='blue')
ax1.plot(sin2_x, sin2_y, color='green')
ax1.plot(sin3_x, sin3_y, color='red')
ax1.plot(tri_x, tri_y, color='orange')
plt.ylabel('{:.1f}Hz'.format(F), rotation=0)
plt.xlim(0, 2*pi)
plt.grid(True)
ax2 = plt.subplot(423)
ax2.plot(rec1_x, rec1_y, color='blue')
plt.ylabel('U', rotation=0)
plt.xlim(0, 2*pi)
plt.grid(True)
ax3 = plt.subplot(425)
ax3.plot(rec2_x, rec2_y, color='green')
plt.ylabel('V', rotation=0)
plt.xlim(0, 2*pi)
plt.grid(True)
ax4 = plt.subplot(427)
ax4.plot(rec3_x, rec3_y, color='red')
plt.ylabel('W', rotation=0)
plt.xlim(0, 2*pi)
plt.grid(True)
# ---------- 線間電圧を表示 ---------- #
ax5 = plt.subplot(422)
ax5.set_title('Line-to-line voltage')
ax5.plot(x_axis, y_axis, color='blue')
ax5.plot(rec1_x, rec_U_V, color='blue')
plt.ylabel('U-V', rotation=0)
plt.xlim(0, 2*pi)
ax6 = plt.subplot(424)
ax6.plot(x_axis, y_axis, color='green')
ax6.plot(rec2_x, rec_V_W, color='green')
plt.ylabel('V-W', rotation=0)
plt.xlim(0, 2*pi)
ax7 = plt.subplot(426)
ax7.plot(x_axis, y_axis, color='red')
ax7.plot(rec3_x, rec_W_U, color='red')
plt.ylabel('W-U', rotation=0)
plt.xlim(0, 2*pi)
plt.pause(0.01)
# グラフを消去
ax1.cla()
ax2.cla()
ax3.cla()
ax4.cla()
ax5.cla()
ax6.cla()
ax7.cla()
# 電圧などの情報を表示
wave.plot_info()
if __name__ == '__main__':
E = 200 # インバーターの入力電圧
F = 0 # 出力周波数
d = 0 # 位相差
for i in range(50):
wave_plot(E, F, Es=0.7, P=20, Ec=1, dif=d)
F += 0.1
d += 1
wave_plot関数を繰り返し実行すればアニメーションになります。
ポイントなど
特別難しいことをしているわけではないですが，三角波をどうやって計算しようかというのは結構悩みました。最初は，範囲を区切って数値解法で計算していましたが，この方法だと周波数を増やしたときにy軸方向のズレが大きくなってしまうことがわかりました。計算の刻み幅を小さくすれば大分マシになりましたが，それでも精度には限界があったり，そもそも式がスマートじゃないということで別の方法を考えることにしました。
三角波の計算式はいくつか調べましたが，最終的には以下のサイトでフーリエ級数展開された式を使うことにしました。
三角波のフーリエ級数展開（参考にしたサイトのリンクと同じ）
この式のおかげでプログラムが大分スッキリしました。ありがとうございます。
矩形波を作るところは，三角波と正弦波のy座標を一つずつ比べて波形を生成しています。もう少し効率良くできそうな気がしますが，今回は処理内容を理解して画面表示することが目標だったので（今からやるのが面倒というのもありますが…）取り敢えずはこのままで。
他にも改善・修正点はいくつかあるので，気が向いたときに直していこうと思います。
改善版（Calculation_wave_new.py）
import numpy as np
import matplotlib.pyplot as plt
class Inverter_wave:
def __init__(self, E, h=2**-7, N=40):
self.E = E # インバーターへの入力
self.h = h # 計算の刻み幅
self.N = N # 三角波の定数
self.pi = np.pi
self.x_num = int(2*self.pi / h) # x方向のの分割数
# 三角変調波（正弦波）を計算
def sin_wave(self, F, Es=200, P_d=0):
"""
F : 変調波信号の周波数[Hz]
Es : 変調波の振幅[V]
P_d : 位相差(Phase difference)[rad]
"""
self.F = F
self.Es = Es
x_list = np.linspace(0, 2*self.pi, self.x_num)
y_list = Es * np.sin((x_list * F) + P_d)
return (x_list, y_list)
# 三角搬送波(triangular wave)を計算
def tri_wave(self, P, Ec=200, P_d=0):
"""
P : キャリア周波数[Hz]
Ec: 搬送波（キャリア）信号の振幅[V]
P_d : 位相差
"""
self.P = P
self.Ec = Ec
def f(x, P, N=self.N):
y = 0
for n in range(1, N+1):
y += ((np.sin(n*self.pi/2) * np.sin(n*(x + P_d)*P)) / n**2) * (8/self.pi) / 3
return y
x_list = np.linspace(0, 2*self.pi, self.x_num)
y_list = f(x_list, P) * Ec
return (x_list, y_list)
# 矩形波(rectangle wave)の座標を計算
def rectangle_wave(self, sin_x, sin_y, tri_x, tri_y, V=200):
# sin_y > tri_y のとき V/2
# sin_y <= tri_y のとき -V/2
x_list = sin_x
y_list = np.where(sin_y > tri_y, V/2, -V/2)
return (x_list, y_list)
def plot_info(self):
# 変調率(modulation factor)
M = self.Es / self.Ec
# 相電圧
Vph = M * self.E / 2
# 線間電圧
V1rms = np.sqrt(3) * Vph
print("インバーターへの入力電圧：%d[V]" % self.E)
if (1 < M):
print("変調率 ：{:.2f} (過変調PWM)".format(M))
else:
print("変調率 ：{:.2f}".format(M))
print("相電圧 ：{:.2f}[V]".format(Vph))
print("線間電圧：{:.2f}[V]".format(V1rms))
print('変調波周波数：{:.2f}[Hz]'.format(self.F))
print('キャリア周波数：{:.2f}[Hz]'.format(self.P))