要旨
Phased Array System Toolboxでmatchedfilterを使うとき、デフォルトではなく必要な設定をきちんとしよう。
環境
- MATLAB R2019b (Update1)
- Phased Array System Toolbox
想定読者
- チャープ、マッチドフィルタについて理解がある方
- MATLABユーザ
はじめに
仕事でPhased Array System Toolboxを使うことがあったので、勉強がてらお家芸であるSARのシミュレータでも作ってみようと思ったのですが、まぁつまずくつまずく。とりあえずマッチドフィルタとるところに一癖あることがわかったので、他の方がつまずかないように共有します。
ポイントはデフォルト設定ではなく、きちんと設定しよう。
やりたいこと
この記事でやることは、下記の2点です。
- チャープ信号を作る
- それをマッチドフィルタを使ってパルス圧縮する。
はい、とても簡単です。
普通に書いてみましょう
まずチャープはLinearFMWaveformというオブジェクトを使って生成することができます。必要なプロパティとしてはPRFとかバンド幅とか、スイープ時間とか、サンプリングレートとなります。
ちょっとだけpalsar-2に似せて、下記設定で作っていきます。
fs = 50e6; % サンプリングレート [Hz]
nof = 1000;
prf = fs/nof % PRF [Hz]
prf = 50000
pw = 10e-6; % スイープ時間 [s]
BW = 20e6; % 帯域 [Hz]
waveform = phased.LinearFMWaveform('SampleRate',fs, ...
'PulseWidth',pw,'SweepBandwidth', BW, ...
'PRF',prf,'NumPulses',1);
確認はplot関数でできます。
plot(waveform)
ここで出てくるのはチャープ部分のみなので、PRF分見たい場合は、いったん出力してから普通にプロットしてください。
t = [0:(nof-1)]/fs; % 時間軸 [s]
plot(t, real(waveform()));
xlabel('時間[s]');
ylabel('振幅 [-]')
よし、いいかんじですね。
次はマッチドフィルタを作っていきます。もうめちゃくちゃ簡単で感動さえしたんですが、こちらもオブジェクトを使っていきます。
coef = getMatchedFilter(waveform); % 参照信号から自動で係数を取得！
mf = phased.MatchedFilter("Coefficients", coef, "SpectrumWindow", "Hamming");
以上！簡単！
作ったフィルタを使って参照信号を圧縮します。
x = waveform(); % 参照信号
compressed = mf(x); % フィルタリング
か、感動的に簡単・・・！！
結果を見てみます。
subplot(2,1,1),plot(real(x))
xlabel('サンプル [-]')
ylabel('振幅')
title('参照信号')
subplot(2,1,2),plot(real(compressed))
xlabel('サンプル [-]')
ylabel('振幅')
title('圧縮後');
できた！！！
sinc波形汚くない？
できた・・・けどなんか変じゃないでしょうか。sinc関数きったな！本当のsincはこんな感じです。
figure, plot(sinc(-10:0.01:10));
はぁ~、きれい。こうなってほしいんですよね。
絶対変なことになってるわ。ということで、時間波形見ていてもわからないのでスペクトルを確認してみましょう。
まずは参照信号のスペクトル
freq = [0:(nof-1)]*fs/nof;
figure, plot(freq, abs(fft(x)));
xlabel('周波数 [Hz]')
うん？・・・あ、いいのか。なんか0Hzから始まるチャープを見慣れてないので、変な感じしましたが、大丈夫ですね。周波数軸の表示違うじゃねぇか、っていう指摘はその通りです。確認できればいいのでこのままでいきます。
それでは、マッチドフィルタ後の信号のスペクトルを見てみます。当然これにハミングがかかった感じじゃないといけないのですが…。
figure, plot(freq, abs(fft(compressed)));
xlabel('周波数 [Hz]');
な・・・なんだこれは…。どうなってるんだ。ほとんどのスペクトル成分が消えてしまっていますね。なぜか5MHzくらいまででスパッと切れていて、そこにハミング窓がかかってしまっている感じです。
MatchedFilterのデフォルトプロパティを見直してみよう
ヘルプドキュメントを見てみますと次のようになっているのがわかりました。
SpectrumRangeのとこ、
Defalut: [0 1e5]
こ・・・これか。coefficientのところでwaveformを丸々入れて作った値を入れていたので、他のプロパティも勝手にオブジェクトが読み取ってよしなにしているのかと自己判断していました。そんなうまい話はないか。でもこれwaveformをいれたら一括でやってくれたら嬉しいなぁ。
関連してその下のSampleRateもこれ入れなきゃいけないんじゃね？ということで、次のように書き直してみます。
mf = phased.MatchedFilter("Coefficients", coef, "SpectrumWindow", "Hamming", ...
"SpectrumRange", [0 BW], ...
"SampleRate", fs);
SpectrumRangeは、[0 fs]と入れかけて、慌てて修正。窓関数を適用する範囲を入れろ、ってなってますね。確かに[0 fs]だと、samplerateで判別できるだろってことか。これでもう一度やってみますと、
compressed = mf(x); % フィルタリング
subplot(2,1,1),plot(real(x))
xlabel('サンプル [-]')
ylabel('振幅')
title('参照信号')
subplot(2,1,2),plot(real(compressed))
xlabel('サンプル [-]')
ylabel('振幅')
title('圧縮後');
お！パルスがずっと細いぞ！ちょっとアップにします。
figure, plot(real(compressed))
axis([400 600 -100 150])
う~ん、まぁこんなもんかな。スペクトルを見てみます。
figure, plot(freq, abs(fft(compressed)));
xlabel('周波数 [Hz]');
ばっちりですね。元のスペクトル幅そのままにハミングがかかった感じになっています。きれいです。
ついでに正負対象のチャープにしたいな…
元々の目的はSARのシミュレーションです。SAR chirpで検索するとこんな感じなんですよ。真ん中あたりに0Hzが来ている。つまりup chirpなら負の周波数から始まって0を通り越して正の周波数で終わる。実運用でどうなってるかは正直知らないのですが、この形式にしたい。そこでLinearFMWaveformのプロパティを確認してみます。
関係ありそうなところはSweepDirectionかなぁ、と確認してみるもUp と Downしかありません。これは負から正にスイープするか、その逆かを決めるものと思われるので関係ないです。マジかよ、設定できないのかなと思いつついろいろ見てたところ、その下のSweepIntervalっていうの見てくださいよ。「'Symmetric'に設定すると-B/2からB/2までスイープします」って書いている！ こ・・・これや。でもSweepIntervalって名称わかりにくくないですか？！意味捉えた後なら解釈もできますが・・・！これがネイティブとの差か・・・。
ということで、これを変えてもう一度やってみます。
waveform = phased.LinearFMWaveform('SampleRate',fs, ...
'PulseWidth',pw,'SweepBandwidth', BW, ...
'SweepInterval', 'Symmetric', ...
'PRF',prf,'NumPulses',1);
plot(waveform)
ごちそうさまです！連動して、matchedfilterの方も変えなくてはいけないですね。coefficientとSpectrumRangeが変わります。
coef = getMatchedFilter(waveform);
mf = phased.MatchedFilter("Coefficients", coef, "SpectrumWindow", "Hamming", ...
"SpectrumRange", [-BW/2 BW/2], ...
"SampleRate", fs);
それでは結果を見てみましょう。
x = waveform();
compressed = mf(x); % フィルタリング
subplot(2,1,1),plot(real(x))
xlabel('サンプル [-]')
ylabel('振幅')
title('参照信号')
subplot(2,1,2),plot(real(compressed))
xlabel('サンプル [-]')
ylabel('振幅')
title('圧縮後');
波形違うな…。またアップにしてみます。
figure, plot(real(compressed))
axis([400 600 -100 150])
スペクトルを見ます。
figure, plot(freq, abs(fft(x)));
xlabel('周波数 [Hz]')
plot(freq, abs(fft(compressed)));
xlabel('周波数 [Hz]');
スペクトルはOKですね。きちんと正負対象になっています。波形の違いは周波数成分が違うので当然ですね、よく考えれば。包絡線取ったら同じになるはずです。
おわりに
ということで、マッチドフィルタのところまではできるようになりました。SARのシミュレーションのためにはまだまだステップがありますが、引続き取り組んで、できた段階でまた共有させていただきます。
謝辞
本記事もeigsさんのlivescript2markdownを使わせていただいてます。