私のPythonプログラムで、IQ信号を正常に受け取れるか試して見ましょう。
弱いCW信号を、Airspy HF+に与えます。
user1@Asrock ~ % airspyhf_rx -z -d -r stdout -f 7.025 -m on > iq.bin user1@Asrock ~ % python3 airspyhf.py
DCオフセットがありますが、良さそうに見えます。
Ham Radio Blog
Airspy HF+のユーザーモードドライバがここで提供されています:
https://github.com/airspy/airspyhf
あなたは、サイトに示された手順に従って、ソフトウェアをビルドすることが出来ます。そして、あなたは以下のものを得るでしょう:
user1@Asrock ~ % ls -l /usr/local/bin/airspy* -rwxr-xr-x 1 root root 17184 Jan 7 22:54 /usr/local/bin/airspyhf_gpio -rwxr-xr-x 1 root root 17504 Jan 7 22:54 /usr/local/bin/airspyhf_info -rwxr-xr-x 1 root root 16608 Jan 7 22:54 /usr/local/bin/airspyhf_lib_version -rwxr-xr-x 1 root root 27320 Jan 7 22:54 /usr/local/bin/airspyhf_rx
ツールの1つを試してみましょう。
user1@Asrock ~ % airspyhf_info AirSpy HF library version: 1.4.2 S/N: .......... Part ID: 0x00000001 Firmware Version: R1.00.00 Available sample rate: 768 kS/s
別のツール、airspyhf_rxは、32ビット浮動小数点数のIQ信号を出力します。
user1@Asrock ~ % airspyhf_rx -z -d -r stdout -f 0.810 -m on | csdr octave_complex_c 1024 120000 --2d | octave -i > /dev/null airspyhf_rx Frequency: -f 0.810000MHz (810000 Hz) 0.768000 MS/s IQ Device Serial Number: .......... Stop with Ctrl-C Streaming at 0.768 MS/s Streaming at 0.768 MS/s Streaming at 0.768 MS/s Streaming at 0.658 MS/s Streaming at 0.680 MS/s Streaming at 0.697 MS/s Streaming at 0.697 MS/s Streaming at 0.711 MS/s Streaming at 0.723 MS/s Streaming at 0.732 MS/s Streaming at 0.739 MS/s Streaming at 0.745 MS/s Streaming at 0.749 MS/s Streaming at 0.753 MS/s
最初の行のCSDRはSDR用のDSPタスクを実行するためのコマンドラインツールです。
GNU octaveは、科学プログラミング言語で、プロッティングとビジュアリゼーションのツールを含んでいます。
さて、私はこのAirspy HF+からのIQ信号で何が出来るでしょうか?
Gqrxは、TCP接続を用いて制御を行うことが可能です。
http://gqrx.dk/doc/remote-control
これは、ローカルホストからです。
そして、これはリモートホストからです。
どちらの場合でも、あなたは、リモート制御の設定として、IPv4-mapped IPv6アドレスを与える必要があります。
user1@Asrock ~/Downloads/gqrx-sdr-2.11.5 % grep -r "ffff:127" . ./src/applications/gqrx/remote_control.cpp:#define DEFAULT_RC_ALLOWED_HOSTS "::ffff:127.0.0.1" user1@Asrock ~/Downloads/gqrx-sdr-2.11.5 % grep -r "AF_INET" . user1@Asrock ~/Downloads/gqrx-sdr-2.11.5 %
これは、Airspy HF+をGqrxを介して制御するPythonプログラムの例です。
import numpy as numpy import matplotlib.pyplot as plt import telnetlib HOST = "192.168.0.90" EXPECTED =b"dummy" TOUT = 0.5 tn = telnetlib.Telnet(HOST, 7356) xvalue = [] yvalue = [] tn.read_until(EXPECTED, timeout=TOUT*5) tn.write(b"M CW 1000\n") print(tn.read_until(EXPECTED, timeout=TOUT).decode('ascii'), end="") for freq_hz in range(500000, 1700000, 1000): freq=b"F "+str(freq_hz).encode('ascii')+b"\n" tn.write(freq) tn.read_until(EXPECTED, timeout=TOUT) tn.write(b"f\n") f = int(tn.read_until(EXPECTED, timeout=TOUT).decode('ascii')) xvalue.append(f) print(f) tn.write(b"l\n") l = float(tn.read_until(EXPECTED, timeout=TOUT).decode('ascii')) yvalue.append(l) print(l) plt.figure(1, figsize=(8,8)) plt.subplot(111) plt.title('Airspy HF+') plt.xlabel('Freq [Hz]') plt.ylabel('Signal Strength [dBm]') plt.grid(True) plt.plot(xvalue, yvalue, color='blue', linewidth=1) plt.show()
私は、ゆっくりとGitHubからGitLabへ移行中です。
ここに、私がRigExpert AA-30.Zeroのために書いた(ほぼ)全てのプログラムがあります。
あなたは、あなたのプロジェクトをGitHubからGitLabへインポートすることが出来ます。
あなたは拡大鏡アイコンをクリックして、地図にズームインすることが出来ます。
import cartopy import cartopy.crs as ccrs import cartopy.feature as cfeature import shapely.geometry as sgeom import cartopy.io.shapereader as shpreader import matplotlib.pyplot as plt import numpy as np import re import datetime def maiden2lonlat(maiden: str): n = len(maiden) maiden = maiden.lower() aaa = ord('a') lon = -180.0 lat = -90.0 lon += (ord(maiden[0])-aaa)*20.0 lat += (ord(maiden[1])-aaa)*10.0 lon += int(maiden[2])*2.0 lat += int(maiden[3])*1.0 if n >= 6: lon += (ord(maiden[4])-aaa) * 5.0/60.0 lat += (ord(maiden[5])-aaa) * 2.5/60.0 if n >= 8: lon += int(maiden[6]) * 5.0/600.0 lat += int(maiden[7]) * 2.5/600.0 return lon, lat def parse_adif(fn): raw = re.split('<eor>|<eoh>',open(fn).read() ) raw.pop(0) raw.pop() logbook =[] for record in raw: qso = {} tags = re.findall('<(.*?):\d+.*?>([^<\t\n\r\f\v]+)',record) for tag in tags: qso[tag[0]] = tag[1] logbook.append(qso) return logbook def mydatetime(date, time): dt = datetime.datetime(int(date[0:4]), int(date[4:6]), int(date[6:8]), \ int(time[0:2]), int(time[2:4]), int(time[4:6])) return dt fnames = ['station1.adif', 'station2.adif'] dt0 = datetime.datetime(2001, 1, 1, 0, 0 , 0) dt9 = datetime.datetime(2099, 12, 31, 23, 59, 59) fig = plt.figure(figsize=(16, 9)) fcount = 0 for fn in fnames: x = [] y = [] r = [] c = [] log = parse_adif(fn) scount = 0 for qso in log: if ('GRIDSQUARE' in qso): dt = mydatetime(qso['QSO_DATE'], qso['TIME_ON']) mytime = qso['TIME_ON'] myhour = float(mytime[0:2]) if dt >= dt0 and dt <=dt9: grid = qso['GRIDSQUARE'] mylon, mylat = maiden2lonlat(grid) if ('APP_PSKREP_SNR' in qso): snr = float(qso['APP_PSKREP_SNR']) print(fcount, scount, grid, mylon, mylat, snr, mytime, myhour) x.append(mylon) y.append(mylat) r.append(50.0+2.0*snr) c.append(myhour/24.0) scount += 1 cent_lon = 152.5 ax = fig.add_subplot(2, 3, 1+3*fcount, projection=ccrs.PlateCarree(central_longitude=cent_lon)) ax.stock_img() ax.gridlines() ax.coastlines() ax.scatter(np.subtract(x,cent_lon), y, c=c, s=r, cmap='hsv', alpha=0.7) cent_lon = 139.7 cent_lat = 35.7 ax = fig.add_subplot(2, 3, 2+3*fcount, projection=ccrs.AzimuthalEquidistant(central_longitude=cent_lon, central_latitude=cent_lat)) ax.stock_img() ax.gridlines() ax.coastlines() ax.scatter(x, y, c=c, s=r, cmap='hsv', alpha=0.7, transform=ccrs.Geodetic()) ax = fig.add_subplot(2, 3, 3+3*fcount, projection=ccrs.PlateCarree()) ax.stock_img() ax.coastlines() ax.set_extent([-130.0, -66.5, 20.0, 50.0], ccrs.Geodetic()) shapename = 'admin_1_states_provinces_lakes_shp' states_shp = shpreader.natural_earth(resolution='110m', category='cultural', name=shapename) for state in shpreader.Reader(states_shp).geometries(): facecolor = [0.9375, 0.9375, 0.859375] edgecolor = 'black' ax.add_geometries([state], ccrs.PlateCarree(), facecolor=facecolor, edgecolor=edgecolor, alpha=0.1) ax.scatter(x, y, c=c, s=r, cmap='hsv', alpha=0.9) fcount += 1 plt.show()
もし、あなたが日本に住んでいるのならこの方が見やすいかもしれません。
cent_lon = 152.5 ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree(central_longitude=cent_lon)) x.append(mylon-cent_lon)
色は1日の中の時間を表しています。点の大きさはSNRです。
カラーマップは、”hsv”です。左端が00Zに、右端が23Zに対応しています。
cent_lon = 139.7 cent_lat = 35.7 ax = fig.add_subplot(1, 1, 1, projection=ccrs.AzimuthalEquidistant(central_longitude=cent_lon, central_latitude=cent_lat)) x.append(mylon) ax.scatter(x, y, c=c, s=r, cmap='hsv', alpha=0.7, transform=ccrs.Geodetic())
ax.set_extent([-135, -66.5, 20, 55], ccrs.Geodetic()) shapename = 'admin_1_states_provinces_lakes_shp' states_shp = shpreader.natural_earth(resolution='110m', category='cultural', name=shapename) for state in shpreader.Reader(states_shp).geometries(): facecolor = [0.9375, 0.9375, 0.859375] edgecolor = 'black' ax.add_geometries([state], ccrs.PlateCarree(), facecolor=facecolor, edgecolor=edgecolor, alpha=0.1)
丸の色と大きさは、受信されたFT8信号のSNRを表しています。
import cartopy.crs as ccrs import matplotlib.pyplot as plt import numpy as np import re import datetime def maiden2lonlat(maiden: str): n = len(maiden) maiden = maiden.lower() print(maiden, n) aaa = ord('a') lon = -180.0 lat = -90.0 lon += (ord(maiden[0])-aaa)*20.0 lat += (ord(maiden[1])-aaa)*10.0 lon += int(maiden[2])*2.0 lat += int(maiden[3])*1.0 if n >= 6: lon += (ord(maiden[4])-aaa) * 5.0/60.0 lat += (ord(maiden[5])-aaa) * 2.5/60.0 if n >= 8: lon += int(maiden[6]) * 5.0/600.0 lat += int(maiden[7]) * 2.5/600.0 return lon, lat ax = plt.axes(projection=ccrs.PlateCarree()) ax.stock_img() fnames = ['station1.adif'] def parse_adif(fn): raw = re.split('<eor>|<eoh>',open(fn).read() ) raw.pop(0) raw.pop() logbook =[] for record in raw: qso = {} tags = re.findall('<(.*?):\d+.*?>([^<\t\n\r\f\v]+)',record) for tag in tags: qso[tag[0]] = tag[1] logbook.append(qso) return logbook def mydatetime(date, time): dt = datetime.datetime(int(date[0:4]), int(date[4:6]), int(date[6:8]), \ int(time[0:2]), int(time[2:4]), int(time[4:6])) return dt dt0 = datetime.datetime(2001, 1, 1, 0, 0 , 0) dt9 = datetime.datetime(2099, 12, 31, 23, 59, 59) x = [] y = [] r = [] c = [] fcount = 0 for fn in fnames: log = parse_adif(fn) scount = 0 for qso in log: if ('GRIDSQUARE' in qso): dt = mydatetime(qso['QSO_DATE'], qso['TIME_ON']) if dt >= dt0 and dt <=dt9: grid = qso['GRIDSQUARE'] mylon, mylat = maiden2lonlat(grid) if ('APP_PSKREP_SNR' in qso): snr = float(qso['APP_PSKREP_SNR']) print(fcount, scount, grid, mylon, mylat, snr) x.append(mylon) y.append(mylat) r.append(50.0+2.0*snr) c.append(-snr*30.0) scount += 1 fcount += 1 plt.scatter(x, y, c=c, s=r, cmap='hsv', alpha=0.5) plt.show()
上の図では、大きさはSNRを、色は1日の時間を表しています。
カラーマップは、”jet”という名称のものを使いました。
左端が00Zに、右端が23Zに対応します。
Cartopyは、地理空間情報を処理して地図を生成するなどのデータ解析を行うためのPythonのパッケージです。
Basemapは、Cartopyプロジェクトに置き換えられます。新規のソフトウェア開発は、可能な限りCartopyを用いるべきです。
import cartopy.crs as ccrs import matplotlib.pyplot as plt import numpy as np ax = plt.axes(projection=ccrs.PlateCarree()) ax.stock_img() n=50 x = 360.0 * np.random.rand(n) - 180.0 y = 120.0 * np.random.rand(n) - 60.0 r = 50.0 * np.random.rand(n) + 20.0 plt.scatter(x, y, c=x, s=r, cmap='hsv', alpha=0.75) plt.show()