· 

numpyでCPFSK再び

  1. # 変数クリア
  2. from IPython import get_ipython
  3. get_ipython().magic('reset -sf')
  4.  
  5. import matplotlib.pyplot as plt
  6. import numpy as np
  7. import scipy.signal as signal
  8.  
  9. # シンボル列を引き延ばす
  10. def extend_symbol_to_simulation_sampling_rate(symbols,t_symbol,t_simulation_sampling):
  11.     # シンボル数
  12.     size_of_symbols=symbols.size
  13.     # 最後のシンボルが終わる時間
  14.     end_time=size_of_symbols*t_symbol
  15.     # 計算タイミング
  16.     t=np.arange(0,end_time,t_simulation_sampling)
  17.     # 各シンボルの開始(終了)時間
  18.     t_periods_of_symbols=np.append(0,np.arange(1,size_of_symbols,1)*t_symbol)
  19.     # 結果データ領域確保
  20.     ex_symbols=np.zeros_like(t)
  21.     # 先頭データは入れておく
  22.     ex_symbols[0]=symbols[0]
  23.     # 各シンボルについて、その期間のインデックスを取得し、データを入れ込む
  24.     for i in range(0,size_of_symbols):
  25.         ex_symbols[np.where(t>t_periods_of_symbols[i])]=symbols[i]
  26.     return ex_symbols,end_time
  27.          
  28. # よく使う変数
  29. pi=np.pi
  30. deg2rad=pi/180.0
  31. twopi=2*pi
  32.  
  33. # 設定
  34. f_carrier=433.92e6            # 搬送波周波数 in Hz
  35. f_deviation=1e6         # 周波数偏差 in Hz
  36. f_sym=1e6                 # シンボルレート in Hz
  37. f_samp_simulation=f_carrier*10     # 計算のサンプリングレート in Hz
  38. phase_initial=45*deg2rad  # 搬送波の初期位相(何でもいい) in rad
  39. symbols_in_base=signal.max_len_seq(7)[0] # シンボル
  40. #symbols_in_base[0::2]=1
  41. #symbols_in_base[1::2]=0
  42.  
  43. # 後で使う変数
  44. t_sym=1/f_sym
  45. t_samp_simulation=1/f_samp_simulation
  46. omega_carrier=twopi*f_carrier
  47. omega_deviation=twopi*f_deviation
  48.  
  49. # シンボルを計算用に拡張
  50. [symbols,t_end]=extend_symbol_to_simulation_sampling_rate(symbols_in_base,t_sym,t_samp_simulation)
  51. symbols=symbols*2-1 # -1 or +1 にする
  52.  
  53. # 計算するタイミング
  54. t=np.arange(0,t_end,t_samp_simulation)
  55.  
  56. # 全計算タイミングにおける瞬間的な角速度
  57. omega_in_moment=omega_carrier+omega_deviation*symbols
  58.  
  59. # 全計算タイミングにおける瞬間的な位相変化量
  60. phase_change_in_moment=omega_in_moment*t_samp_simulation
  61.  
  62. # 全計算タイミングでの位相
  63. phase_t=np.cumsum(phase_change_in_moment)+phase_initial
  64.  
  65. # 波形生成
  66. amp_t=np.exp(1j*phase_t)
  67.  
  68. np.save('output/fsk_o_t',t)
  69. np.save('output/fsk_o_amp_t',amp_t)
  70.  
  71. # プロット
  72. fig=plt.figure()
  73. fig.add_subplot(3,1,1)
  74. plt.plot(t,symbols)
  75. fig.add_subplot(3,1,2)
  76. plt.plot(t,amp_t)
  77. fig.add_subplot(3,1,3)
  78. plt.plot(t,amp_t)
  79. plt.xlim(t_sym*7.95,t_sym*8.05) # ちょうどシンボルの変わり目を表示

 

 

まぁ、画像じゃ周波数がどこで変わったのかはわかりません。