合成データおよび取得データを使用した P25 スペクトルセンシング

スクリプトを開く

この例では、周期定常性特徴検出を使用して、P25 信号 [1] を含む、異なる変調スキームの信号を識別する方法を示します。定義する事例の信号は次の 4 つです。定義する事例の信号は次の 4 つです。ノイズのみの信号、C4FM、CQPSK および任意の種類の信号。この例では、異なる SNR 値をもつ信号および取得した実際の P25 信号に検出アルゴリズムを適用し、4 つの種類うちの 1 つとして信号を分類します。すべての事例で検出アルゴリズムが成功したという結果が視覚的に表示されます。

Project 25 (P25)

Project 25 (P25 または APCO-25) は、北米における連邦、州、郡および地域の公共保安機関が使用するデジタル無線通信の一連の標準規格です。緊急事態が発生した場合、このプロトコルスイートにより行政機関と共済対応チームとの通信が可能になります。その意味では、P25 は欧州の TErrestrial Trunked RAdio (TETRA) [2] プロトコルと同じ役割を果たします。ただし、両者間の標準規格には相互運用性はありません。北米では、P25 は公共保安、防衛、公共サービスおよび商用での利用において幅広く使用されています [1]。

プロジェクト 25 は 2 段階で展開されます。段階 1 では、P25 は C4FM を使用します。C4FM は 4 値周波数変調 (compatible 4 level frequency modulation) の頭文字です。最も単純な形式において、シンボルの表現に 4 つの異なる周波数を使用する 4FSK 変調の特別な種類です。段階 1 ではこの変調スキームを使用して 12.5 kHz チャネル上でデジタル情報を送信します。

段階 2 では、互換性のある直交位相偏移変調 (CQPSK) 変調形式を使用して 6.25 kHz チャネル上でデジタル情報を送信します。CQPSK 変調は、本質的には符号化が対称となる pi/4 差動直交位相偏移変調 (pi/4 DQPSK) であり、次の図に示すように -135 度、-45 度、+45 度および +135 度の位相変化値を使用します。

この図では、必ず、赤い点の次の状態は緑の点、緑の点の次の状態は赤い点になります。シンボルごとのデータ転送速度とデータビットは同一ですが、2 つの変調スキームの主な違いは C4FM がシンボルを表すために周波数シフトを使用する点です。これにより固定された振幅信号が提供されます。それとは対照的に、CQPSK は位相シフトを使用してシンボルを表します。これが信号に振幅成分を与えます。

周期定常性特徴検出

変調認識と信号分類は、20 年以上にわたって多くの研究が行われてきたテーマです。分類スキームは一般的に、次の 2 つの大きなカテゴリのいずれかに分けられます。尤度ベース (LB) アプローチおよび特徴ベース (FB) アプローチ [3]。周期定常性特徴検出は FB 手法の 1 つであり、通信信号は正確には定常とされず、周期定常としてより適切にモデル化されるという事実に基づくものです [4]。

周期定常プロセスは、時間とともに循環的に変化する統計的な特性をもつ信号です [5]。このような周期性は、サンプリング、スキャン、変調、多重化および符号化などの処理のため、明確な形で信号に現れます。この信号の周期性特質を利用して未知の信号の変調スキームを判別することができます [4]。

周期定常性特徴検出は、信頼性の高いスペクトルセンシング手法です。ノイズとは異なり、変調された情報は周期定常プロセスだからです。結果として、周期検出器は低 SNR 環境下でも正常に動作します。

ノイズのみの事例

ノイズのみの事例では、電力が 1 dBW のホワイトガウスノイズの (4*N) 行 1 列のベクトルを生成します。1/(4*N) は、commP25ssca.m のスペクトル自己相関関数 (SAF) の算出に使用される循環分解能です。

N = 4096;
input = wgn(4*N,1,1);

時間領域のスペクトル自己相関関数を使用して信号 x(t) の周期定常性特徴を解析します。入力信号に対してスペクトル自己相関関数 commP25ssca.m を実行します。この関数は、ストリップスペクトル相関アルゴリズム (SSCA) [3] の一時的な平滑化法を使用して、理想的なスペクトル自己相関関数を推定します。これは FFT ベースの時間平滑化アルゴリズムです。このアルゴリズムの実装の詳細は、[6] を参照してください。

プロット関数 commP25plot.m を実行します。この手順では、3 次元図であるスペクトル自己相関関数について説明します。x 軸は巡回周波数 (alpha) -1 ～ 1 を表します。y 軸はスペクトル周波数 (f) -0.5 ～ 0.5 を、z 軸 (Sx) は各 (alpha、f) ペアのスペクトル自己相関関数に対応する振幅をそれぞれ表します。巡回分解能 dalpha は 1/T です。ここで T はデータの観測時間です。スペクトル分解能 df は 1/Tw です。ここで Tw は複素復調を計算するための時間枠です [7]。T は Tw より大きいので、dalpha は df より小さくなります。alpha がゼロに等しくない場合、SAF 値はほぼゼロとなることに注意してください。

% 64 represents the window time Tw, 4*N represents the observation time T
[Sx,alphao,fo] = commP25ssca(input,1,1/64,1/(4*N));
fig1 = figure('Position',figposition([5 40 40 40]));
commP25plot(Sx,alphao,fo);

commP25decision_noise.m は入力信号がノイズのみを含んでいるかどうかを判別し、commP25decision_c4fm.m は入力信号が C4FM 信号かどうかを判別します。commP25decision_cqpsk.m は、入力信号が CQPSK 信号かどうかを判別します。これらの判定は、SAF のピークの位置に基づいています。この例では、コードは P25 信号が存在しないことを正しく結論付けています。

[c,d] = size(Sx);
[Ades,Index] = sort(Sx(:),'descend');   % sort Sx by its element and store in Ades
[Ridx,Cidx]  = ind2sub(size(Sx),Index); % corresponding row index and column index
leng = length(Ades);

noise_decision = commP25decision_noise(Ades,Ridx,Cidx,leng,c,d);
if noise_decision == 0
    c4fm_decision = commP25decision_c4fm(Ades,Ridx,leng,c);
    if c4fm_decision == 0
        commP25decision_cqpsk(Ades,Ridx,Cidx,leng,c,d);
    end
end

There is no P25 signal.

合成データを使用した C4FM の事例

[8] に従い、次の変調構造体は C4FM 出力信号を生成します。

ナイキストパルス整形の条件を満たす通常のレイズドコサインフィルターは、符号間干渉を最小限にします。レイズドコサインフィルターのパラメーターは、[8] のフィルターの仕様に基づいて選択されています。具体的には、このレイズドコサインフィルターのアップサンプリング係数は 4 で、ロールオフ係数は 0.2 です。C4FM 方式も逆 sinc フィルターをレイズドコサインフィルターの後に呼び出し、P25 受信機の積分およびダンプフィルターの sinc 応答を補正します。FM 変調器は 600 Hz の偏差があります。

設計判断へのノイズの影響を観測するため、-3 dB、3 dB および無限大 dB の SNR 値で検出を実行します。

% The length of input bits is N. The length of the output bits must also be
% N
x = randi([0,3],N,1);
sym = 2*x-3;                    % integer input

% Raised Cosine Filter
sampsPerSym = 4;                % Upsampling factor
% Design raised cosine filter with given order in symbols. Apply gain to
% the unit energy filter to obtain max amplitude of 1.
rctFilt = comm.RaisedCosineTransmitFilter(...
    'Shape', 'Normal', ...
    'RolloffFactor', 0.2, ...
    'OutputSamplesPerSymbol', sampsPerSym, ...
    'FilterSpanInSymbols', 60, ...
    'Gain', 1.9493);
c4fm_init = rctFilt(sym);

shape2 = 'Inverse-sinc Lowpass';
d2 = fdesign.interpolator(2, shape2);
intrpltr = design(d2, 'SystemObject', true);
c4fm_init = intrpltr(c4fm_init);

% Baseband Frequency Modulator
Fs = 4800;
freqdev = 600;
int_x = cumsum(c4fm_init)/Fs;
c4fm_output = exp(1i*2*pi*freqdev*int_x);
y = c4fm_output(1:N); % Ideal case, SNR = infinity
y1 = awgn(y,3);       % SNR = 3 dB
y2 = awgn(y,-3);      % SNR = -3 dB

対応するスペクトル自己相関関数が計算されプロットされます。SNR 値が減るにつれ、SAF ピークがより不明確になることに注目してください。

[Sx0,alphao0,fo0] = commP25ssca(y,1,1/64,1/(4*N));
[Sx1,alphao1,fo1] = commP25ssca(y1,1,1/64,1/(4*N));
[Sx2,alphao2,fo2] = commP25ssca(y2,1,1/64,1/(4*N));
fig2 = figure('Position',figposition([5 40 80 40]));
subplot(131);
commP25plot(Sx0,alphao0,fo0);
title('Ideal case');
subplot(132);
commP25plot(Sx1,alphao1,fo1);
title('SNR = 3 dB');
subplot(133);
commP25plot(Sx2,alphao2,fo2);
title('SNR = -3 dB');

このセクションは、前で行った手順と同じ手順で SNR 値ごとに分類結果を取得します。関数 commP25decision.m は、考えられるすべての入力信号タイプについてスペクトルセンシングによる分類を実行します。

commP25decision(Sx0); % Ideal case

There is signal present. Checking for presence of C4FM.
This is C4FM.

commP25decision(Sx1); % SNR = 3 dB

There is signal present. Checking for presence of C4FM.
This is C4FM.

commP25decision(Sx2); % SNR = -3 dB

There is signal present. Checking for presence of C4FM.
This is C4FM.

合成データを使用した CQPSK の事例

[8] に従い、次の変調構造体は CQPSK 出力信号を生成します。

CQPSK 変調器は同相および直交 (I および Q) 部分で構成されます。入力ビットは、ルックアップテーブル [8] によって処理され、5 値 I/Q 信号を発生させます。ルックアップテーブルの仕様は pi/4 DQPSK と同等であるため、この例では DQPSK 変調器 System object™ を使用してこのルックアップテーブルを実装します。I/Q 信号は次に、前の事例で説明したレイズドコサインフィルターでフィルター処理されます。

% The size of input bits is 2*N, the size of output is 4*N
x = randi([0,1],2*N,1);

% Create a DQPSK modulator System object with bits as inputs, phase
% rotation of pi/4 and Gray-coded constellation
dqpskMod = comm.DQPSKModulator(pi/4,'BitInput',true);
% Modulate and filter
modout = dqpskMod(x);
release(rctFilt);
cqpsk_output = rctFilt(modout);
y = cqpsk_output;   % Ideal case, SNR = infinity
y1 = awgn(y,3);     % SNR = 3 dB
y2 = awgn(y,-3);    % SNR = -3 dB

対応するスペクトル自己相関関数を計算してプロットします。

[Sx0,alphao0,fo0] = commP25ssca(y,1,1/64,1/(4*N));
[Sx1,alphao1,fo1] = commP25ssca(y1,1,1/64,1/(4*N));
[Sx2,alphao2,fo2] = commP25ssca(y2,1,1/64,1/(4*N));
fig3 = figure('Position',figposition([5 40 80 40]));
subplot(131);
commP25plot(Sx0,alphao0,fo0);
title('Ideal case');
subplot(132);
commP25plot(Sx1,alphao1,fo1);
title('SNR = 3 dB');
subplot(133);
commP25plot(Sx2,alphao2,fo2);
title('SNR = -3 dB');

以下のコード出力は、3 つの異なる SNR 値の CQPSK 検出の結果を示しています。

commP25decision(Sx0); % Ideal case

There is signal present. Checking for presence of C4FM.
This is NOT C4FM. Checking for presence of CQPSK.
This is CQPSK.

commP25decision(Sx1); % SNR = 3 dB

There is signal present. Checking for presence of C4FM.
This is NOT C4FM. Checking for presence of CQPSK.
This is CQPSK.

commP25decision(Sx2); % SNR = -3 dB

There is signal present. Checking for presence of C4FM.
This is NOT C4FM. Checking for presence of CQPSK.
This is CQPSK.

合成データを使用した P25 信号以外の事例

この事例では、任意の信号タイプを定義し、P25 周期定常性検出器で処理して P25 信号かどうかを判定します。

FIR 等リップルローパスフィルターを設計して、乱数入力に適用します。この場合、信号にノイズを追加しないでください。他の信号タイプを使用してみて、周期定常特性検出器によりそれらを分類します。

bcoeffs = firpm(200,[0 0.2 0.22 1],[1 1 0 0]); % Set N to achieve 40 dB rejection
input = randn(N,1);
y = filter(bcoeffs,1,input);

その後、スペクトル自己相関関数を計算し、プロットします。各信号の異なる変調特性によって SAF が大幅に異なることに注目してください。

[Sx,alphao,fo] = commP25ssca(y,1,1/64,1/(4*N));
fig4 = figure('Position',figposition([5 40 40 40]));
commP25plot(Sx,alphao,fo);

同じ手順に従って、分類結果を取得します。

commP25decision(Sx);

There is signal present. Checking for presence of C4FM.
This is NOT C4FM. Checking for presence of CQPSK.
This is NOT CQPSK either, so it is not a P25 signal.

取得データを使用した C4FM の事例

この事例では、取得した実際の C4FM 信号に検出アルゴリズムを適用します。この信号は、446 MHz の P25 無線で送信され、USRP® 無線で受信され、その後 MATLAB® により capturedc4fm.mat に保存されます。同じ手順に従って、分類結果を取得します。

load capturedc4fm.mat;
y = y(1:4*N);
agc = comm.AGC;
y = 0.1*agc(y);
[Sx,alphao,fo] = commP25ssca(y,1,1/64,1/(4*N));
fig5 = figure('Position',figposition([5 40 40 40]));
commP25plot(Sx,alphao,fo);
commP25decision(Sx);

There is signal present. Checking for presence of C4FM.
This is C4FM.

まとめ

この例では、周期定常性検出を使用して異なる変調スキームの信号の判別方法を示しました。ここでのアルゴリズムは、スペクトル自己相関関数内のピーク位置に基づいて信号を分類します。周期定常性検出器はそのノイズに対する復元性から、エネルギー検出器などの検出器よりも利点があります。

付録

この例では以下のスクリプトと補助関数が使用されています。

commP25ssca.m

commP25plot.m

commP25decision.m

commP25decision_noise.m

commP25decision_c4fm.m

commP25decision_cqpsk.m

参考文献

P25 Technology Interest Group: https://www.project25.org/
TETRA <https://en.wikipedia.org/wiki/Terrestrial_Trunked_Radio https://en.wikipedia.org/wiki/Terrestrial_Trunked_Radio>
E. C. Like, "Non-Cooperative Modulation Recognition Via Exploitation of Cyclic Statistics".MS Thesis.2007
E. C. Like, V. D. Chakravarthy, P. Ratazzi, and Z. Wu, "Signal Classification in Fading Channels Using Cyclic Spectral Analysis", EURASIP Journal on Wireless Communications and Networking, Volume 2009, 2009.
W. A. Gardner, A. Napolitano and L. Paura, "Cyclostationarity: Half a century of research", Signal Processing, Vol. 86, No. 4, pp. 639-697, 2006.
E. L. Da Costa, "Detection and Identification of Cyclostationary Signals".MS Thesis.1996.
Antonio F. Lima, Jr., "Analysis of Low Probability of Intercept (LPI) Radio Signals using Cyclostationary Processing". MS Thesis.2002.
TIA Standard Project 25: <https://tiaonline.org/what-we-do/standards/ https://tiaonline.org/what-we-do/standards/>

合成データおよび取得データを使用した P25 スペクトル センシング