計測  ·  31日 8月 2024

    BLEを下から見ていく(5)

     

    だいぶ一般的になっているBLEについて。すごく概念的なことはあちこちに書いてある。が、下から上まで(むしろ下だけ)知りたいのに、資料少なすぎ。なので自力で調べる無茶な挑戦。第５回＝CH37に飛んでるアドバタイズをじゃんじゃん見てみる

     

    Adalm-plutoを操作するPythonコードを変更してたっくさんのデータをとれるように
    サンプリングレートを8MS/sにして、取得データ長を10Mptsにする。
    1. import numpy as np
    2. import adi
    3. import matplotlib.pyplot as plt
    4. sample_rate = 8e6 # Hz
    5. center_freq = 2402e6 # Hz
    6. num_samps = 10000000 # number of samples per call to rx()
    7. sdr = adi.Pluto("ip:192.168.2.1")
    8. sdr.sample_rate = int(sample_rate)
    9. # Config Rx
    10. sdr.rx_lo = int(center_freq)
    11. #sdr.rx_rf_bandwidth = int(sample_rate)
    12. sdr.rx_rf_bandwidth = int(2e6)
    13. sdr.rx_buffer_size = num_samps
    14. sdr.gain_control_mode_chan0 = 'manual'
    15. sdr.rx_hardwaregain_chan0 = 6.0 # dB, increase to increase the receive gain, but be careful not to saturate the ADC
    16. # Clear buffer just to be safe
    17. for i in range (0, 10):
    18.     raw_data = sdr.rx()
    19. # Receive samples
    20. rx_samples = sdr.rx()
    21. print(rx_samples)
    22. t=np.linspace(0,1/sample_rate*num_samps,num_samps)
    23. np.save('o_t',t)
    24. np.save('o_rx_t',rx_samples)
    25. # Calculate power spectral density (frequency domain version of signal)
    26. psd = np.abs(np.fft.fftshift(np.fft.fft(rx_samples)))**2
    27. psd_dB = 10*np.log10(psd)
    28. f = np.linspace(sample_rate/-2, sample_rate/2, len(psd))
    29. # Plot time domain
    30. plt.figure(0)
    31. plt.plot(t,np.real(rx_samples))
    32. plt.plot(t,np.imag(rx_samples))
    33. plt.xlabel("Time")
    34. # Plot freq domain
    35. plt.figure(1)
    36. plt.plot(f/1e6, psd_dB)
    37. plt.xlabel("Frequency [MHz]")
    38. plt.ylabel("PSD")
    39. plt.xlim(-2,2)
    40. plt.show()
    まぁこんなかんじ。
    で、こうなる。
    nRF52840モジュールをAdalm-plutoにピタピタに置いているんで一番振幅が大きいのがきっとそれで、それ以外にも色々いるよね。大きいのも小さいのもできるだけ見ていくことにする。 で、全部復調復号してみると、エラーが出まくるので、データがありそうなところだけ復調復号していくことにする。で、最初は人力で範囲を指定して、、、って思ったけど、まぁまぁﾒﾝﾄﾞｸｻｲ。振幅とパタンで自動判定しよう！いやーこの辺の判断が難しいよねー。自動化に要する時間と人力で要する時間と一体どっちが早いのか。ある程度人力でやってみないと自動化の必要性がわからないってのは間違いない。今回はもう何回かやってたので、こりゃ人力じゃしんどいってわかったので、自動化や！
    1. import numpy as np
    2. import matplotlib.pyplot as plt
    3. from scipy import signal
    4.  
    5. def search_index_by_pattern(data,pattern):
    6.     ds=data.size
    7.     ps=pattern.size
    8.     cyc=ds-ps+1
    9.     ret=-1
    10.     if cyc>0:
    11.         for i in range(cyc):
    12.             data_part=data[i:i+ps]
    13.             if np.all(data_part==pattern):
    14.                 ret=i
    15.                 break
    16.     return ret
    17.  
    18. f_sym=1e6 #
    19. rx_all=np.load('o_rx_t.npy')
    20. t_all=np.load('o_t.npy')
    21. t_all=t_all-t_all[0]
    22. t_all_0=t_all[:-1:]
    23. t_all_1=t_all[1::]
    24. t_samp=np.mean(t_all_1-t_all_0)
    25. f_samp=1/t_samp
    26. print(f_samp/1e6) # should be 8MHz
    27.  
    28. # フィルター係数生成
    29. freq_cutoff=f_sym*2
    30. w_cutoff=freq_cutoff/(f_samp/2)
    31. b,a=signal.butter(1,w_cutoff,'lowpass')
    32.  
    33. # フィルター適用
    34. rx_all=signal.lfilter(b,a,rx_all)
    35.  
    36. rx_all_amp=np.abs(rx_all)
    37. rx_all_amp_dig=np.where(rx_all_amp>5,1,0)
    38. rx_all_amp_dig_0=rx_all_amp_dig[:-1:]
    39. rx_all_amp_dig_1=rx_all_amp_dig[1::]
    40. rx_all_amp_dig_10=np.append(rx_all_amp_dig_1-rx_all_amp_dig_0,0)
    41. rx_all_amp_dig_rising_edge=np.where(rx_all_amp_dig_10>0)[0]
    42. rx_all_amp_dig_falling_edge=np.where(rx_all_amp_dig_10<0)[0]
    43. print(rx_all_amp_dig_rising_edge.size)
    44. print(rx_all_amp_dig_falling_edge.size)
    45.  
    46. pa_and_aa=np.array([0,1,0,1,0,1,0,1,0,1,1,0,1,0,1,1,0,1,1,1,1,1,0,1,1,0,0,1,0,0,0,1,0,1,1,1,0,0,0,1])
    47.  
    48. for i in range(rx_all_amp_dig_rising_edge.size):
    49.     if rx_all_amp_dig_falling_edge[i]-rx_all_amp_dig_rising_edge[i]>40:
    50.         rx=rx_all[rx_all_amp_dig_rising_edge[i]:rx_all_amp_dig_falling_edge[i]]
    51.         t=t_all[rx_all_amp_dig_rising_edge[i]:rx_all_amp_dig_falling_edge[i]]
    52.         # demod
    53.         rx_t0=rx[:-1:]
    54.         rx_t1=rx[1::]
    55.         rx_dem=-1*np.angle(rx_t1/rx_t0)
    56.         rx_dem=np.append(0,rx_dem)
    57.         # normalization
    58.         rx_dem=rx_dem/(f_sym/f_samp*np.pi*1)
    59.         # digitization
    60.         rx_dig=np.where(rx_dem>0,0,1)
    61.         # counter top value
    62.         counter_top=8-1
    63.         counter_half=int(counter_top/2)
    64.         cnt_list=[]
    65.         cnt_list.append(int(0))
    66.         for i in range(rx_dig.size-1):
    67.             dig_pres=rx_dig[i+1]
    68.             dig_prev=rx_dig[i]
    69.             if (dig_pres!=dig_prev):
    70.                 cnt_list.append(int(0))
    71.             else:
    72.                 if (cnt_list[i]<counter_top):
    73.                     cnt_list.append(cnt_list[i]+int(1))
    74.                 else:
    75.                     cnt_list.append(int(0))
    76.         cnt=np.array(cnt_list)
    77.         smp_pos=np.where(cnt==counter_half)[0]
    78.         rd_dig_smp=rx_dig[smp_pos]
    79.         pp_and_aa_pos=search_index_by_pattern(rd_dig_smp,pa_and_aa)
    80.         if(pp_and_aa_pos>0):
    81.             print("t_start=",end="")
    82.             print(t[0])
    83.             print("t_end=",end="")
    84.             print(t[-1])
    85.             print(pp_and_aa_pos)
    86.             rd_dig_smp=rd_dig_smp[pp_and_aa_pos:]
    87.             pos=0
    88.             while pos<rd_dig_smp.size:
    89.                 for n in range(8):
    90.                     print(rd_dig_smp[pos],end='')
    91.                     pos=pos+1
    92.                     if pos==rd_dig_smp.size:
    93.                         break
    94.                 print("")
    95.                 
    こんなんで。
    振幅が5以上の範囲で、復調復号して、プリアンブルとAccess Addressを既知のパタンとして検索して存在すれば結果を表示する。ってしている。ndarrayからパタンマッチするインデックスを探すってのが機能としてすでにありそうなのに見つからなかったので自力で作成している。まぁ絶妙にローテク。それと、ndarrayを直接printすると、"["がついてたり、勝手に改行したり後で面倒なので、これまたローテクで、ついでに8シンボルごとに改行してprintしている。ちなみに今度は行数が多すぎて、VSCodeのコンソールの行バッファが足りなくなるので、そこは都度、よしなに。ちなみに個人的にはリダイレクトでファイルに出力するって手段を好んで使っている。
    で、こんな感じで出力される。
    で、数えると24フレームとれたようだ。実際の受信波形を見ると、振幅大と振幅中を数えると24個ある。ちなみに振幅小部分は拡大すると、
    って感じで振幅が動いているので、きっと隣のチャンネルだね。で、隣はアドバタイズじゃないからAccess Addressが見つからなかったってことでしょ(そもそもFSK復調できたのかどうかもわからんが)。
    で、解析してみると、こんなんで、
    結局、
    6024DE30CD58ABDF0201061AFF590002150112233445566778899AABBCCDDEEFF0DEADBEEFCA54545E
    ってのと、
    4225E6BDBB8333081EFF060001092022A0DA83015D56DBE5C18DB0E031F8F60DC69DA1F4E008C09DAD78
    ってのが繰り返し来てるってのがわかった。
    で、1つ目はSeeed XIAO BLE nRF52840のbeaconで2つ目は何かわからん(いや、いつかわかるようになりたい)。ところで、これをどう解釈するのかってので、Android用にnRF Connectっていうツールがある(アイホン用ももちろんあるだろうけど林檎法人は商品戦略が好きじゃないので敢えて知らんふうを装っている)ので利用してみると、
    RAWを押すと、
    つまり
    Header : 60
    PDU Length : 24 (36bytes)
    AdvA : DE30CD58ABDF
    AdvData : 
    02 01 06
    1A FF 590002150112233445566778899AABBCCDDEEFF0DEADBEEFCA
    CRC : 54545E
    で、AdvDataについて、Assigned Numbers(https://www.bluetooth.com/ja-jp/specifications/assigned-numbers/)を見ると、
    ということらしい。で、電波からデータを取り出すってのは一部できたんだけど、Complete local Nameって飛んできてないんだけど、どうなってるん？って思ってる。
    カテゴリ： 現実