最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

高分解能スペクトル解析

この例では、チャネライザーと呼ばれることもある効率のよいフィルター バンクを使用して高分解能スペクトル解析を実行する方法を説明します。比較のために、従来の平均修正ピリオドグラム (ウェルチ) 手法も一緒に説明します。

スペクトル解析における分解能

ここでの分解能とは、近傍に位置する 2 つのスペクトル成分を識別できる能力を表します。分解能はスペクトルの計算に使用する時間領域セグメントの長さによって異なります。修正ピリオドグラムに見られるように、ウィンドウ処理を時間領域セグメントで使用すると、使用するウィンドウのタイプも分解能に影響を与えます。

異なるウィンドウによる典型的なトレードオフには、分解能とサイドローブ減衰の関係が挙げられます。箱型ウィンドウは最も高い分解能を提供しますが、サイドローブ減衰が非常に小さくなります (~14 dB)。サイドローブ減衰が小さいとスペクトル成分がウィンドウ処理によって埋もれてしまうので、望ましくありません。ハン ウィンドウは周波数分解能を低く抑えることによって適切なサイドローブ減衰を提供します。カイザーなどのパラメーター化可能なウィンドウでは、ウィンドウ パラメーターを変更してトレードオフを制御できます。

平均修正ピリオドグラム (ウェルチ法) を使用する代わりにアナログ スペクトル アナライザーの動作をエミュレートするフィルター バンク アプローチを使って、より高い分解能の推定を達成することができます。主な考え方としては、フィルター バンクを使って信号を異なる周波数ビンに分割し、各サブバンド信号の平均強度を計算することです。

フィルター バンク ベースのスペクトル推定

この例では、箱型ウィンドウと同じ分解能を取得するために 512 の異なるバンドパス フィルターを使用する必要があります。512 のバンドパス フィルターを効率よく実装するため、ポリフェーズ解析フィルター バンク (チャネライザーとも呼ばれる) を使用します。これは、Fs/N (N は目的の周波数分解能、この例では 512) の帯域幅をもつプロトタイプのローパス フィルターを使用し、FIR 間引きの実装と同じようにフィルターをポリフェーズ形式で実装することで動作します。間引きの場合のようにすべての分岐の結果を追加するのではなく、各分岐は N 点の FFT への入力として使用されます。FFT の各出力が変調されたローパス フィルターに対応し、したがってバンドパス フィルターを実装していることがわかります。フィルター バンク アプローチの主な短所は、ポリフェーズ フィルターによって計算が増加し、そのフィルターの状態のために変化する信号への適応が遅くなることです。詳細については、書籍、『Multirate Signal Processing for Communications Systems』(fredric j. harris 著、Prentice Hall PTR、2004) を参照してください。

この例では、全体を通してスペクトル推定の 100 の平均値が使用されています。サンプリング周波数は 1 MHz に設定されています。フレームのサンプル数は 64 と想定されており、スペクトル推定を実行するにはバッファーする必要があります。

NAvg = 100;
Fs = 1e6;
FrameSize = 64;
NumFreqBins = 512;
filterBankRBW = Fs/NumFreqBins;

dsp.SpectrumAnalyzer は、Method が適切に設定されていれば、フィルター バンク ベースのスペクトル推定器を実装します。内部では、ポリフェーズ フィルター処理と FFT を実装する (マルチキャリア通信などスペクトル解析以外の他の応用事例にも使用できる) dsp.Channelizer が使用されます。

filterBankSA = dsp.SpectrumAnalyzer(...
    'Method','Filter bank',...
    'NumTapsPerBand',24,...
    'SampleRate',Fs,...
    'RBWSource','Property',...
    'RBW',filterBankRBW,...
    'SpectralAverages',NAvg,...
    'PlotAsTwoSidedSpectrum',false,...
    'YLimits',[-150 50],...
    'YLabel','Power',...
    'Title','Filter bank Power Spectrum Estimate',...
    'Position',[50 375 800 450]);

テスト信号

この例では、サンプル数が 64 のフレームでテスト信号を取得しています。スペクトル解析では、フレームが大きくなるほど分解能も向上します。

テスト信号は 2 つの正弦波とホワイト ガウス ノイズから構成されています。周波数ビンの数、振幅、周波数およびノイズ パワーの各値を変更することが推奨されます。

sinegen = dsp.SineWave('SampleRate',Fs,'SamplesPerFrame',FrameSize);

最初のテスト ケース

まず、振幅が 1 と 2、周波数が 200 kHz と 250 kHz の正弦波のフィルター バンク スペクトル推定をそれぞれ計算します。ホワイト ガウス ノイズの平均強度 (分散) は 1e-12 です。スペクトル推定で片側のノイズ フロアが -114 dBm と正確に示されていることに注意してください。

release(sinegen)
sinegen.Amplitude = [1 2];
sinegen.Frequency = [200000 250000];

noiseVar = 1e-12;
noiseFloor = 10*log10((noiseVar/(NumFreqBins/2))/1e-3); % -114 dBm onesided
fprintf('Noise Floor\n');
fprintf('Filter bank noise floor = %.2f dBm\n\n',noiseFloor);

timesteps = 10 * ceil(NumFreqBins / FrameSize);
for t = 1:timesteps
    x = sum(sinegen(),2) + sqrt(noiseVar)*randn(FrameSize,1);
    filterBankSA(x);
end

release(filterBankSA)
Noise Floor
Filter bank noise floor = -114.08 dBm

スペクトル推定器を使用した数値計算

dsp.SpectrumEstimator は、フィルター バンク スペクトル推定を計算するために使用できます。

より長いフレームをスペクトル推定器に供給するために、スペクトル推定を計算する前に 512 のサンプルがバッファーで収集されます。この例では使用していませんが、バッファーでオーバーラッピングを使用すれば、特定のデータ集合から取得される平均の数を増やすことができます。

filterBankEstimator = dsp.SpectrumEstimator(...
    'Method','Filter bank',...
    'NumTapsPerBand',24,...
    'SampleRate',Fs,...
    'SpectralAverages',NAvg,...
    'FrequencyRange','onesided',...
    'PowerUnits','dBm');

buff    = dsp.AsyncBuffer;

release(sinegen)

timesteps = 10 * ceil(NumFreqBins / FrameSize);
for t = 1:timesteps
    x     = sum(sinegen(),2) + sqrt(noiseVar)*randn(FrameSize,1);
    write(buff,x);      % Buffer data
    if buff.NumUnreadSamples >= NumFreqBins
        xbuff = read(buff,NumFreqBins);
        Pfbse = filterBankEstimator(xbuff);
    end
end

異なる方法を使用したスペクトル推定の比較

振幅が 1 と 2、周波数が 200 kHz と 250 kHz の正弦波のウェルチ スペクトル推定とフィルター バンク スペクトル推定をそれぞれ計算します。ホワイト ガウス ノイズの平均強度 (分散) は 1e-12 です。

release(sinegen)
sinegen.Amplitude = [1 2];
sinegen.Frequency = [200000 250000];

filterBankSA.RBWSource = 'Auto';
filterBankSA.Position = [50 375 400 450];

welchSA = dsp.SpectrumAnalyzer(...
    'Method','Welch',...
    'SampleRate',Fs,...
    'SpectralAverages',NAvg,...
    'PlotAsTwoSidedSpectrum',false,...
    'YLimits',[-150 50],...
    'YLabel','Power',...
    'Title','Welch Power Spectrum Estimate',...
    'Position',[450 375 400 450]);

noiseVar = 1e-12;

timesteps = 500 * ceil(NumFreqBins / FrameSize);
for t = 1:timesteps
    x = sum(sinegen(),2) + sqrt(noiseVar)*randn(FrameSize,1);
    filterBankSA(x);
    welchSA(x);
end

release(filterBankSA)

RBW = 488.28;
hannNENBW = 1.5;

welchNSamplesPerUpdate = Fs*hannNENBW/RBW;
filterBankNSamplesPerUpdate = Fs/RBW;

fprintf('Samples/Update\n');
fprintf('Welch       Samples/Update = %.3f Samples\n',welchNSamplesPerUpdate);
fprintf('Filter bank Samples/Update = %.3f Samples\n\n',filterBankNSamplesPerUpdate);

welchNoiseFloor = 10*log10((noiseVar/(welchNSamplesPerUpdate/2))/1e-3);
filterBankNoiseFloor = 10*log10((noiseVar/(filterBankNSamplesPerUpdate/2))/1e-3);

fprintf('Noise Floor\n');
fprintf('Welch noise floor       = %.2f dBm\n',welchNoiseFloor);
fprintf('Filter bank noise floor = %.2f dBm\n\n',filterBankNoiseFloor);
Samples/Update
Welch       Samples/Update = 3072.008 Samples
Filter bank Samples/Update = 2048.005 Samples

Noise Floor
Welch noise floor       = -121.86 dBm
Filter bank noise floor = -120.10 dBm

ウェルチ スペクトル推定とフィルター バンク ベースのスペクトル推定の両方で、200 kHz と 250 kHz の 2 つのトーンが検出されました。フィルター バンク ベースのスペクトル推定の方がトーンの分離に優れています。同じ分解能帯域幅 (RBW) について、平均修正ピリオドグラム (ウェルチ) 手法はスペクトルの計算に 3073 のサンプルを必要とするのに比べ、フィルター バンク ベースの推定では 2048 のサンプルを必要とします。フィルター バンク スペクトル推定で片側のノイズ フロアが -120 dBm と正確に示されていることに注意してください。

異なるウィンドウを使用した修正ピリオドグラムの比較

使用するウィンドウ (箱型またはハン) のみが異なる 2 つのスペクトル アナライザーを考えてみます。

rectRBW = Fs/NumFreqBins;
hannNENBW = 1.5;
hannRBW = Fs*hannNENBW/NumFreqBins;

rectangularSA = dsp.SpectrumAnalyzer(...
    'SampleRate',Fs,...
    'Window','Rectangular',...
    'RBWSource','Property',...
    'RBW',rectRBW,...
    'SpectralAverages',NAvg,...
    'PlotAsTwoSidedSpectrum',false,...
    'YLimits',[-50 50],...
    'YLabel','Power',...
    'Title','Welch Power Spectrum Estimate using Rectangular window',...
    'Position',[50 375 400 450]);

hannSA = dsp.SpectrumAnalyzer(...
    'SampleRate',Fs,...
    'Window','Hann',...
    'RBWSource','Property',...
    'RBW',hannRBW,...
    'SpectralAverages',NAvg,...
    'PlotAsTwoSidedSpectrum',false,...
    'YLimits',[-150 50],...
    'YLabel','Power',...
    'Title','Welch Power Spectrum Estimate using Hann window',...
    'Position',[450 375 400 450]);

release(sinegen)
sinegen.Amplitude = [1 2]; % Try [0 2] as well
sinegen.Frequency = [200000 250000];

noiseVar = 1e-12;
timesteps = 10 * ceil(NumFreqBins / FrameSize);
for t = 1:timesteps
    x = sum(sinegen(),2) + sqrt(noiseVar)*randn(FrameSize,1);
    rectangularSA(x);
    hannSA(x);
end
release(rectangularSA)
release(hannSA)

箱型ウィンドウはサイドローブ減衰が小さくなりますが、狭いメインローブを提供します。対照的に、ハン ウィンドウではメインローブはより広くなりますが、これと引き換えにはるかに大きなサイドローブ減衰を提供します。メインローブが広くなるのは 250 kHz 地点で特に顕著になります。両方のウィンドウとも正弦波がある周波数の辺りで大きなロールオフを示しています。これにより、ノイズ フロアより上のパワーの低い対象信号がマスクされます。この問題は、フィルター バンクの場合にはほぼ見られません。

振幅を [1 2] ではなく [0 2] に変更することは、実際には単一の 250 kHz の正弦波がノイズと共に存在することを表します。箱型ウィンドウは 200 kHz の正弦波が干渉しない場合に特に動作がよくなるため、このケースは興味深いものになっています。これは、250 kHz が 1 MHz を均等に分割する 512 の周波数の 1 つであることが理由です。その場合、FFT 固有の周波数サンプリングによって生成される時間領域の複製により、パワー スペクトル計算に使用される時間制限のあるデータ セグメントが完全な周期的拡張をもつようになります。一般に、任意の周波数をもつ正弦波においては、こうなりません。信号干渉の影響の受けやすさと連動する、正弦波の周波数の依存性は、修正ピリオドグラム アプローチのもう 1 つの欠点です。

分解能帯域幅 (RBW)

各アナライザーの分解能帯域幅は、入力の長さがわかった時点で計算できます。RBW は電力成分が計算される帯域幅を示します。つまり、パワー スペクトル推定の各要素は、推定要素に対応する周波数を中心とする幅 RBW の帯域幅全体にわたり、パワーをワット、dBW または dBm の単位で表します。パワー スペクトル推定の各要素のパワー値は、パワー密度を RBW 値の範囲内にある周波数帯域で積分することによって求められます。パワーは細かいグリッド (小さい帯域幅) 上で計算されるため、RBW が小さいほど分解能が高くなります。箱型ウィンドウがすべてのウィンドウの中で最も高い分解能をもちます。カイザー ウィンドウの場合、RBW は使用されるサイドローブ減衰によって異なります。

fprintf('RBW\n')
fprintf('Welch-Rectangular  RBW = %.3f Hz\n',rectRBW);
fprintf('Welch-Hann         RBW = %.3f Hz\n',hannRBW);
fprintf('Filter bank        RBW = %.3f Hz\n\n',filterBankRBW);
RBW
Welch-Rectangular  RBW = 1953.125 Hz
Welch-Hann         RBW = 2929.688 Hz
Filter bank        RBW = 1953.125 Hz

振幅を [0 2] に設定する、つまり 250 kHz の単一の正弦波がある場合、RBW とノイズ フロアの関係は興味深いものになります。予想されるノイズ フロアは 10*log10((noiseVar/(NumFreqBins/2))/1e-3) またはおよそ -114 dBm です。箱型ウィンドウに対応するスペクトル推定では予想されるノイズ フロアとなりますが、ハン ウィンドウを使用したスペクトル推定のノイズ フロアは予想値よりも約 2 dBm 高くなります。これは、スペクトル推定は 512 周波数点で計算されますが、パワー スペクトルは特定のウィンドウの RBW で積分されることが理由です。箱型ウィンドウの場合、RBW はちょうど 1 MHz/512 となるため、スペクトル推定には各周波数ビンについて独立したパワーの推定値が含まれます。ハン ウィンドウの場合、RBW が大きいため、スペクトル推定には 1 つの周波数ビンから次のビンにオーバーラップしたパワーが含まれます。このオーバーラップしたパワーにより各ビンのパワー値が増加し、ノイズ フロアが上昇します。この量は、次のように解析的に計算できます。

hannNoiseFloor = 10*log10((noiseVar/(NumFreqBins/2)*hannRBW/rectRBW)/1e-3);
fprintf('Noise Floor\n');
fprintf('Hann noise floor = %.2f dBm\n\n',hannNoiseFloor);
Noise Floor
Hann noise floor = -112.32 dBm

相互に隣接する正弦波

分解能の問題を説明するため、次のケースを見てみましょう。正弦波周波数は 200 kHz と 205 kHz に変更されています。フィルター バンク推定の正確性は維持されます。ウィンドウベースの推定器に注目すると、ハン ウィンドウではメインローブがより広いため、箱型ウィンドウの推定と比較すると 2 つの正弦波の区別が難しくなっています。実際には、2 つの推定のいずれも特に正確ではありません。205 kHz は基本的に 200 kHz と区別できる限度に過ぎないことに注意してください。より近い周波数では、3 つの推定器はいずれも 2 つのスペクトル成分を分離することができません。より近い成分を分離する唯一の方法は大きいサイズのフレームをもつことであり、フィルター バンク推定器の場合は NumFrequencyBands の数値を上げることです。

release(sinegen)
sinegen.Amplitude = [1 2];
sinegen.Frequency = [200000 205000];

filterBankSA.RBWSource = 'Property';
filterBankSA.RBW       =  filterBankRBW;
filterBankSA.Position = [850 375 400 450];

noiseVar = 1e-10;
noiseFloor = 10*log10((noiseVar/(NumFreqBins/2))/1e-3); % -94 dBm onesided
fprintf('Noise Floor\n');
fprintf('Noise floor = %.2f dBm\n\n',noiseFloor);

timesteps = 500 * ceil(NumFreqBins / FrameSize);
for t = 1:timesteps
    x  = sum(sinegen(),2) + sqrt(noiseVar)*randn(FrameSize,1);
    filterBankSA(x);
    rectangularSA(x);
    hannSA(x);
end

release(filterBankSA)
release(rectangularSA)
release(hannSA)
Noise Floor
Noise floor = -94.08 dBm

パワーの低い正弦波成分の検出

次は、非常に小さい振幅をもつ 170 kHz の 3 番目の正弦波を追加して前述のシナリオを再実行します。この 3 番目の正弦波は箱型ウィンドウの推定とハン ウィンドウの推定では完全に無視されます。フィルター バンク推定では分解能およびトーンの分離が向上するため、これら 3 つの正弦波がはっきりと区別できます。

release(sinegen)
sinegen.Amplitude = [1e-5 1 2];
sinegen.Frequency = [170000 200000 205000];

noiseVar = 1e-11;
noiseFloor = 10*log10((noiseVar/(NumFreqBins/2))/1e-3); % -104 dBm onesided
fprintf('Noise Floor\n');
fprintf('Noise floor = %.2f dBm\n\n',noiseFloor);

timesteps = 500 * ceil(NumFreqBins / FrameSize);
for t = 1:timesteps
    x  = sum(sinegen(),2) + sqrt(noiseVar)*randn(FrameSize,1);
    filterBankSA(x);
    rectangularSA(x);
    hannSA(x);
end

release(filterBankSA)
release(rectangularSA)
release(hannSA)
Noise Floor
Noise floor = -104.08 dBm

Simulink バージョンのスペクトル アナライザー

上記と同様の高分解能スペクトル推定の設定は、Simulink では Spectrum Analyzer ブロックを使用してモデル化できます。SpectrumAnalyzerFilterBank モデルは、Simulink を使用したフィルター バンク ベースのスペクトル推定の高い分解能を示すものであり、ウェルチ法に比べてノイズ フロアが低いことがわかります。

次のケースを見てみましょう。正弦波は 170 kHz、200 kHz、205 kHz の 3 つで、振幅は [1e-5 1 2] です。最初の正弦波は箱型ウィンドウの推定では完全に無視されます。フィルター バンク推定では分解能および 3 つのトーンの分離が向上します。

モデル例

open_system('SpectrumAnalyzerFilterBank');
sim('SpectrumAnalyzerFilterBank');

モデルを閉じる

bdclose('SpectrumAnalyzerFilterBank');

Simulink バージョンのスペクトル推定器

上記の高分解能スペクトル推定の数値計算は、Simulink でも Spectrum Estimator ブロックを使用してモデル化できます。dspfilterbankspectrumestimation モデルは、Simulink を使用したフィルター バンク ベースの高分解能スペクトル推定を示したものであり、ウェルチ法に比べてノイズ フロアが低いことがわかります。結果の可視化に配列プロットが使用されています。配列プロットは、スペクトル推定をプロットするための便利な方法を提供します。値は dBm 単位で表示されますが、簡単にワット単位と dBW 単位に置き換えることができます。

モデル例

open_system('dspfilterbankspectrumestimation');
sim('dspfilterbankspectrumestimation');

モデルを閉じる

bdclose('dspfilterbankspectrumestimation');