前回、ソフトウエア無線受信用USBドングルRTL-SDRのデータをPythonで読み取ったものをScilabへ「輸出」。これでScilabでもRTL-SDRのデータを処理できる（リアルタイムじゃないけど）とほくそ笑みました。今回は「輸出」に使ったPythonスクリプトに手を入れデシメーション機能など追加。益々繁盛？
※「手習ひデジタル信号処理」投稿順 Indexはこちら
Pythonからの複素数データの出力方法をチョイとScilabに合わせるだけでRTL-SDRのIQ信号をScilabで読み取れるようになりました。ただ「生データ」を素のままで送るのも芸がない上、サイズもデカいです。今回は、前回作成スクリプトのオプションとしてデシメーション処理した後のデータを出力する機能を追加してみました。これでデータ量は何分の１だかに削減可能。ついでに取得したデータのPSDプロットを表示するオプションも追加してみました。
Scipy使ったデシメーション処理
デシメーションそのものは一定の間隔でデータを間引く処理なのでPythonでは一撃ですが、メンドイのがデシメーション・フィルタです。デシメーションによりサンプリング周波数は 1/n (nは整数）に低下してしまうので、信号帯域を新しいサンプリング周波数の1/2（ナイキスト周波数）以下に制限するためにLPFかける必要があるからです。
ところがPythonでの定番Scipyモジュールの信号処理機能の中にデシメーション・フィルタも「ついでにかけてくれる」機能を発見。これ使えばお楽じゃね。
Signal processing (scipy.signal)
デシメーション機能については以下です。
それにしてもデシメーション、恐ろしい言葉なんだが。。。
機能追加したPythonスクリプト
前回スクリプトに対して追加した機能は以下です。
-
- “–D f ” オプション、デシメーション・ファクタfを指定します。このオプションがついていたらデシメーション済のScilab可読なファイルも生成します。ファイル名は、”–FNAME” で指定したファイル名が”test.csv”だとすると”test_decimated.csv” のようになります。
- “–G opt” オプション、opt=psdのとき元のIQ信号のPSD（パワースペクトル密度）プロット、opt=decivatedのときデシメーション済のIQ信号のPSDプロットを表示します。
なお、コマケー話ですが以下変更点です。
-
- “–N n” サンプルデータ点数の指定。前回は素のnどおりの点数でしたが、今回からは n * 1024 と解釈するようにいたしました。
#!/usr/bin/python
# coding: utf-8
### @file pyrtlsdr2sci.py
### @brief save RTL-SDR IQ signal for scilab
###
### @author: Joseph Halfmoon
### @date: June 7, 2023
"""
Save RTL-SDR IQ signal for scilab v0.2
input
=====
--FS s sample rate (default=2.048[MHz])
--FC c center_freq (defalut=80[MHz]
--N n number of samples (defalut=256, means 256*1024)
--FNAME x File name to save (defalut=toScilab.csv)
--G opt 'psd' 'decimated'
--D f decimation factor
-V VERSION
"""
import argparse
import csv
import matplotlib.pyplot as plt
import numpy as np
import os.path
import scipy.signal
import sys
from rtlsdr import RtlSdr
versionSTR = "pyrtlsdr2sci.py v0.2"
def parse2float(st):
"""parse string to float
Returns at Fail: None
"""
try:
work = float(st)
except:
return None
return work
def chkVal(arg, dval):
vl = parse2float(arg)
if vl is None:
nLen = len(arg)
if nLen > 1:
bdy = arg[0:nLen-1]
lastChar = arg[nLen-1]
vl = parse2float(bdy)
if vl is None:
vl = dval
else:
lastChar = "0"
return (vl, lastChar)
def freqIN(arg, dval):
prefixes = 'mkMG'
vl, lastChar = chkVal(arg, dval)
pfx = prefixes.find(lastChar)
if pfx == 0:
ex = 1e-3
elif pfx == 1:
ex = 1e3
elif pfx == 2:
ex = 1e6
elif pfx == 3:
ex = 1e9
else:
ex = 1
return vl * ex
def parse2int(st):
"""parse string to int
Returns at Fail: None
"""
try:
work = int(st)
except:
return None
return work
class rtlsdrBuffer(object):
"""rtlsdr buffer class
read & save RTL-SDR IQ signal
Fs : sampling Freq[Hz]
Fc : center Freq[Hz]
"""
def __init__(self, Fs, Fc, N):
self.sdr = RtlSdr()
self.sdr.sample_rate = Fs
self.sdr.center_freq = Fc
self.sdr.freq_correction = 60
self.sdr.gain = 'auto'
self.N = N
self.buffer = np.zeros(self.N, dtype=np.complex128)
def read(self):
self.buffer = self.sdr.read_samples(self.N)
self.sdr.close()
def writeSCI(self, fnam, x):
self.fnam = fnam
np.savetxt(self.fnam, x, fmt='%.18e %+.18ej')
def plotPSD(self, x):
plt.psd(x, NFFT=1024, Fs=self.sdr.sample_rate/1e6, Fc=self.sdr.center_freq/1e6)
plt.xlabel('Frequency [MHz]')
plt.ylabel('Relative power [dB]')
plt.show()
def decimation(self, N):
self.decimated = scipy.signal.decimate(self.buffer, N)
def main():
"""main program.
process options
"""
parser = argparse.ArgumentParser(description=versionSTR)
parser.add_argument('--FS', nargs=1, help='sample rate (default=2.048[MHz])')
parser.add_argument('--FC', nargs=1, help='center_freq (defalut=80[MHz])')
parser.add_argument('--N', nargs=1, help='number of samples (defalut=256, means 256*1024)')
parser.add_argument('--FNAME', nargs=1, help='File name to save (defalut=toScilab.csv)')
parser.add_argument('--G', nargs=1, help='select options: psd, decimated')
parser.add_argument('--D', nargs=1, help='decimation factor (defalut=8)')
parser.add_argument('-V', dest='VERSION', help='Show Version, then exit', action='store_true', default=False)
args = parser.parse_args()
if args.VERSION:
print( versionSTR )
sys.exit(0)
Fs = 2.048e6
if args.FS is not None:
Fs = freqIN(args.FS[0], Fs)
Fc = 80e6
if args.FC is not None:
Fc = freqIN(args.FC[0], Fc)
N = 256*1024
if args.N is not None:
tmpN = parse2int(args.N[0])
if tmpN is not None:
N = tmpN * 1024
FNAME = "toScilab.csv"
if args.FNAME is not None:
FNAME = args.FNAME[0]
f1, fext = os.path.splitext(FNAME)
f1 = f1 + "_decimated"
FNAME_D = f1 + fext
sdr = rtlsdrBuffer(Fs, Fc, N)
sdr.read()
sdr.writeSCI(FNAME, sdr.buffer)
if args.D is not None:
D = parse2int(args.D[0])
if D is not None:
sdr.decimation(D)
sdr.writeSCI(FNAME_D, sdr.decimated)
if args.G is not None:
if args.G[0].startswith('psd'):
sdr.plotPSD(sdr.buffer)
if args.G[0].startswith('decimated'):
sdr.plotPSD(sdr.decimated)
sys.exit(0)
if __name__ == "__main__":
main()
実行結果
以下のコマンドラインから上記スクリプトを実行することでScilab 可読のファイル２つを生成できます。
また、デシメーション済のデータのPSDプロットも以下のように表示されます。
元のIQファイルのサンプリング周波数2.048MHzに対して、８：１のデシメーションかけているのでデシメーション後のIQファイルのサンプリング周波数は256kHzです。