· 

Pythonでフィルタ(2)

フィルタの構成方法と適用方法について記録を残す。ASK変調前のベースバンド波形に対して狙った特性のフィルタをかけるにはどうしたらいいのか。

ちなみに、今週は遠出して日常とは違う生活を送っていたので、頭の中はすっかりリセットされて、この取り組みの最終目的が何だったのか思い出せない。とりあえずだらだらやる。

で、ベースバンド波形に、、、ですが、この取り組み単体の成果として、AM変調機能のついているSGやらFGを使って、フィルタを通した後の波形を出力することができるようになります。(ただし、ベースバンド波形を出す任意波形発生器(通常、FGのArbitraryWaveform機能を使う)のメモリーの制約を受けて、波形の長さやサンプリング周波数が制限されます。)

一つの方法として、Pythonでフィルタ(1)で生成したフィルタのかかった変調波形を復調すればいいってことになります。復調するには変調波形をヒルベルト変換して複素信号の大きさを出す(AMなので)方法と、最初から複素数で変調して大きさを出す方法があります。今回は完全に計算だけの話なので最初から複素数で変調することにします。

Pythonでフィルタ(1)では実部だけを取り出していたので複素数にします。

実部虚部を同じチャート上に表示します。

3,4番目のチャートに実部虚部がプロットされているはずだけどつぶれてて見えない。

matplotlibの拡大機能で4番目のチャートのt=0.0002付近を見ると、

となっているので、うまくできているようです。では復調してみます。

ていうか、これだけ。

で復調波形もプロットします。

こうなる。そして拡大する。

すばらしい!

これでも全然問題ないんだけど、やっぱりまどろっこしいので、ベースバンドに直フィルタをかけます。ちなみに理論そっちのけの結果重視で取り組んでいきます。

まず現状の結果をファイルに書き出します。

これで、tとamp_demod_tがo_t_amp_demod_t.npzに保存されます。

で、ベースバンドにフィルタをかけるとなると、どうするか、

変調波形にバンドパスフィルタを適用する際、

\[ \zeta=\frac{1}{2Q} \] \[ f_{cut}=\left[\left(1-\zeta\right) f_{carrier},\left(1+\zeta\right) f_{carrier}\right] \]

としたわけなので、

\[ \zeta=\frac{1}{2Q} \] \[ f_{cut}=\zeta f_{carrier} \]

でlowpassフィルタをかければいけそうな気がします(徹底的に理論回避)。

で、さっそく、

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

# よく使う変数
pi=np.pi
deg2rad=pi/180.0
rad2deg=180.0/pi
twopi=2*pi
 
# シンボル列を引き延ばす
def extend_symbol_to_simulation_sampling_rate(symbols,t_symbol,t_simulation_sampling):
    # シンボル数
    size_of_symbols=symbols.size
    # 最後のシンボルが終わる時間
    end_time=size_of_symbols*t_symbol
    # 計算タイミング
    t=np.arange(0,end_time,t_simulation_sampling)
    # 各シンボルの開始(終了)時間
    t_periods_of_symbols=np.append(0,np.arange(1,size_of_symbols,1)*t_symbol)
    # 結果データ領域確保
    ex_symbols=np.zeros_like(t)
    # 先頭データは入れておく
    ex_symbols[0]=symbols[0]
    # 各シンボルについて、その期間のインデックスを取得し、データを入れ込む
    for i in range(0,size_of_symbols):
        ex_symbols[np.where(t>t_periods_of_symbols[i])]=symbols[i]
    return ex_symbols,end_time
        
f_sampling_simulation=1e9
f_carrier=500e3

#generate baseband waveform
f_symbol=10e3                 # シンボルレート in Hz
t_sampling_simulation=1/f_sampling_simulation
symbols_in_base=np.array([1,0,1,0,1,0,1,0,1,0,1,0])
baseband_t,t_end=extend_symbol_to_simulation_sampling_rate(symbols_in_base,1/f_symbol,t_sampling_simulation)
t=np.arange(0,t_end,t_sampling_simulation)
number_of_samples=t.size

#フィルタ生成
q_value=10
zeta=1/(2*q_value)
fcut=np.array([0+zeta*f_carrier])
print(fcut)
wcut=fcut/(f_sampling_simulation/2)
b,a=signal.butter(1,wcut,'lowpass')
w,h=signal.freqz(b,a,worN=2**20)
freqs=w*f_sampling_simulation/twopi
h_abs=np.abs(h)
amps_db=20*np.log10(h_abs)
angles=np.angle(h)*rad2deg

#フィルタ適用
baseband_filtered_t=signal.lfilter(b,a,baseband_t)

#変調フィルタ復調した波形をファイルからロード
npz = np.load('o_t_amp_demod_t.npz')
t_from_file=npz['arr_0']
amp_demod_t_from_file=npz['arr_1']

#プロット
fig = plt.figure()
ax1=plt.subplot(3,1,1)
ax1.grid(True)
ax1.plot(freqs,amps_db, '-')
ax1.set_xlim(400e3,600e3)
ax1.set_ylim(-50,0)
ax1.set_ylabel('Amplitude [dB]')
ax1_2=ax1.twinx()
ax1_2.grid(True)
ax1_2.plot(freqs,angles, 'r-')
ax1_2.set_xlim(0,1e6)
ax1_2.set_ylim(-180,180)
ax1_2.set_yticks(np.arange(-180,181,step=90))
ax1.set_xlabel('frequency [Hz]')
ax1_2.set_ylabel('Phase [deg.]')

ax2=plt.subplot(3,1,2)
ax2.plot(t,baseband_t)

ax3=plt.subplot(3,1,3)
ax3.plot(t,baseband_filtered_t)
ax3.plot(t_from_file,amp_demod_t_from_file)

fig.tight_layout()
plt.show()

3つめのチャートが

ベースバンド-->ローパスフィルタ

としたものと

ベースバンド-->変調-->バンドパスフィルタ-->復調

としたものを重ね書きしたもの、うまく逝っているみたい。

拡大してみる。

ほぼ重なっているので、うまくいっているのだろう。

ベースバンド-->ローパスフィルタのほうが滑らかな波形になっているのはなぜかはわからない。が、こっちのほうがよさそう。計算もまどろっこしくない。

やっているうちに思い出した、この取り組みの最終目的をもう忘れないようにメモっておくと、

FPGAにベースバンド信号を入力してフィルタをかけてDAC回路に出力する

ってのでした。これでFGの最大ポイント数に泣かされずに好きなフィルタをかけたベースバンド波形を(ほぼ)リアルタイムで生成できる。

しかし、道のりは遠い。、、、死ぬまでに終わるかわからん。

目標=長生き!からの孤独死!٩(。•̀ω•́。)و