フィルタの構成方法と適用方法について記録を残す。で、FIRフィルタを試してみる。
ところで、フィルタってホントいろいろあって、学問としてどういう風に学ぶのが正解なのかわからん。必要な時に必要なことを学ぶとして、理論はあいかわらずそっちのけですすめる。
FIRフィルタについては、あっちこっちに説明があるので、ここでは一切触れない。まじで、一切。こちらのscipy.signalのリファレンスをのんびり眺めると、firwinってのがきっとそうだろう。
で、これで思った通りの波形になるのかどうか、、、
まずはベースバンドにフィルタしてどうなるのか、
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=1e7
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=zeta*f_carrier
print(fcut)
wcut=fcut/(f_sampling_simulation/2)
fircoeff=signal.firwin(1001,wcut)
#fircoeff = signal.firwin(numtaps=21, cutoff=fcut, fs=f_sampling_simulation)
print(fircoeff)
w,h=signal.freqz(fircoeff,1,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(fircoeff,1,baseband_t)
#変調フィルタ復調した波形をファイルからロード
npz = np.load('o_t_baseband_filtered_t.npz')
t_from_file=npz['arr_0']
baseband_filtered_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,baseband_filtered_t_from_file)
fig.tight_layout()
plt.show()
なんか違う雰囲気です。(燈色が期待する波形)
では変調波でもやってみる。
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
def generate_ask(f_carrier,f_sym,f_samp_generation,phase_carrier_initial,symbols_in_base):
t_sym=1/f_sym
t_samp_generation=1/f_samp_generation
omega_carrier=twopi*f_carrier
# シンボルを計算用に拡張
[symbols,t_end]=extend_symbol_to_simulation_sampling_rate(symbols_in_base,t_sym,t_samp_generation)
# 計算するタイミング
t=np.arange(0,t_end,t_samp_generation)
# 波形生成
amp_t=np.exp(1j*omega_carrier*t+phase_carrier_initial)*symbols
return t,symbols,amp_t
f_sampling_simulation=1e7
f_carrier=500e3
#generate ask waveform
f_symbol=10e3 # シンボルレート in Hz
t_sampling_simulation=1/f_sampling_simulation
phase_initial=0*deg2rad # 搬送波の初期位相(何でもいい) in rad
symbols_in_base=np.array([1,0,1,0,1,0,1,0,1,0,1,0])
t,baseband_t,amp_t=generate_ask(f_carrier,f_symbol,f_sampling_simulation,phase_initial,symbols_in_base)
#amp_t=np.real(amp_t) #後で復調するので複素数のままにする
number_of_samples=t.size
#フィルタ生成
q_value=10
zeta=1/(2*q_value)
fcut=np.array([f_carrier-zeta*f_carrier,f_carrier+zeta*f_carrier])
print(fcut)
wcut=fcut/(f_sampling_simulation/2)
fircoeff=signal.firwin(1001,wcut,pass_zero=False)
#fircoeff = signal.firwin(numtaps=21, cutoff=fcut, fs=f_sampling_simulation)
print(fircoeff)
w,h=signal.freqz(fircoeff,1,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
#フィルタ適用
amp_filtered_t=signal.lfilter(fircoeff,1,amp_t)
#復調
amp_demod_t=np.abs(amp_filtered_t)
np.savez('o_t_amp_demod_t_fir',t,amp_demod_t)
#プロット
fig = plt.figure()
ax1=plt.subplot(4,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_2.set_xlabel('frequency [Hz]')
ax1_2.set_ylabel('Phase [deg.]')
ax2=plt.subplot(4,1,2)
ax2.plot(t,baseband_t)
ax3=plt.subplot(4,1,3)
ax3.plot(t,np.real(amp_t))
ax3.plot(t,np.imag(amp_t))
ax4=plt.subplot(4,1,4)
ax4.plot(t,np.real(amp_filtered_t))
ax4.plot(t,np.imag(amp_filtered_t))
ax4.plot(t,amp_demod_t)
fig.tight_layout()
plt.show()
やっぱちゃうなー。
矩形波に対してFIRフィルタを適用している例をwebで見つけられなかったので、これが正しいかどうかはわからないけど、欲しい波形とは違う。理論には触れない。とりあえず期待と違う。
地面テラス星6アーマーガアはてんねんのヘイラッシャで楽勝(あまごいは必要)。なぜ今まで気づかなかったのか、、、
コメントをお書きください