Airspy HF+からのIQ信号をプロットしてみる

私の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オフセットがありますが、良さそうに見えます。

Airspy HF+のホストソフトウェア

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からのUDPオーディオストリーム(WSJT-Xへ)

サーバー側:

Airspy HF+は、Linuxマシンに接続され、Gqrxもそこで走っています。


クライエント側:

WSJT-Xは、macOSマシン上で走ります。

$ python3 gqrxUdp5.py


私の関連する過去の記事を、検索して下さい。

今の所、4つあります。

Gqrxのリモート制御の設定

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()

PSK Reporterで2つの受信局を比較する

あなたは拡大鏡アイコンをクリックして、地図にズームインすることが出来ます。

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()

PSK Reporterから自前で地図を描く(2)

もし、あなたが日本に住んでいるのならこの方が見やすいかもしれません。

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)

PSK Reporterから自前で地図を描く

丸の色と大きさは、受信された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に対応します。

PSK ReporterとCartopyパッケージ

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()