pmusic

MUSIC アルゴリズムを使用した疑似スペクトル

構文

[S,wo] = pmusic(x,p)

[S,wo] = pmusic(x,p,wi)

[S,wo] = pmusic(___,nfft)

[S,wo] = pmusic(___,'corr')

[S,fo] = pmusic(x,p,nfft,fs)

[S,fo] = pmusic(x,p,fi,fs)

[S,fo] = pmusic(x,p,nfft,fs,nwin,noverlap)

[___] = pmusic(___,freqrange)

[___,v,e] = pmusic(___)

pmusic(___)

説明

[S,wo] = pmusic(x,p) では、多重信号分類 (MUSIC) アルゴリズムが実装され、入力信号 x の疑似スペクトル推定 S と、疑似スペクトルが評価される正規化周波数 (ラジアン/サンプル単位) のベクトル wo が返されます。入力引数 p を使用して、信号部分空間の次元を指定できます。

[S,wo] = pmusic(x,p,wi) では、ベクトル wi で指定された正規化周波数において計算された疑似スペクトルが返されます。ベクトル wi は 2 つ以上の要素を持たなければなりません。そうでない場合は、関数が nfft として解釈するためです。

[S,wo] = pmusic(___,nfft) では、疑似スペクトルの推定に使用される、整数値の FFT 長 nfft が指定されます。この構文には、前の構文の入力引数を任意に組み合わせて含めることができます。

[S,wo] = pmusic(___,'corr') では、入力引数 x が信号データの行列ではなく、相関行列として解釈されます。この構文では、x は正方行列でなければなりません。また、すべての固有値は非負でなければなりません。

[S,fo] = pmusic(x,p,nfft,fs) では、ベクトル fo で指定された周波数において計算された疑似スペクトルが返されます (Hz 単位)。サンプル レート fs を Hz 単位で指定します。

[S,fo] = pmusic(x,p,fi,fs) では、ベクトル fi で指定された周波数において計算された疑似スペクトルが返されます。ベクトル fi は 2 つ以上の要素を持たなければなりません。そうでない場合は、関数が nfft として解釈するためです。

[S,fo] = pmusic(x,p,nfft,fs,nwin,noverlap) では、ウィンドウ nwin およびオーバーラップ長 noverlap を使用して入力データ x をセグメント化することで、疑似スペクトル S が返されます。

[___] = pmusic(___,freqrange) では、fo または wo に含める周波数値の範囲が指定されます。

[___,v,e] = pmusic(___) では、ノイズの固有ベクトルの行列 v が、関連する固有値を含むベクトル e と共に返されます。

出力引数なしに pmusic(___) を使用すると、現在の Figure ウィンドウに疑似スペクトルがプロットされます。

例

サンプリングの指定がない pmusic

以下の例では、信号部分空間に 2 つの実数正弦波成分があると仮定して、ある信号ベクトル x を解析します。この場合、各実数正弦波は 2 つの複素指数の和で表されるので、信号部分空間の次元は 4 です。

n = 0:199;
x = cos(0.257*pi*n) + sin(0.2*pi*n) + 0.01*randn(size(n));
pmusic(x,4)      % Set p to 4 because there are two real inputs

サンプリング周波数と部分空間次元を指定

以下の例では、最小値を 10% 上回る固有値のカットオフで同じ信号ベクトル x を解析します。p(1) = Inf に設定することで、しきい値パラメーター p(2) に基づいて信号/ノイズ部分空間が決定されるようになります。引数 nwin を使用して長さ 7 の固有ベクトルを指定し、サンプリング周波数 fs を 8 kHz に設定します。

rng default
n = 0:199;
x = cos(0.257*pi*n) + sin(0.2*pi*n) + 0.01*randn(size(n));
[P,f] = pmusic(x,[Inf,1.1],[],8000,7); % Window length = 7
plot(f,20*log10(abs(P)))
xlabel 'Frequency (Hz)', ylabel 'Power (dB)'
title 'Pseudospectrum Estimate via MUSIC', grid on

相関行列の入力

スペクトル密度を推定するために、正定値相関行列 R を入力します。既定の 256 サンプルを使用します。

R = toeplitz(cos(0.1*pi*(0:6))) + 0.1*eye(7);
pmusic(R,4,'corr')

corrmtx で生成した信号データ行列の入力

関数 corrmtx を使用してデータから生成した信号データ行列、Xm を入力します。

n = 0:699;
x = cos(0.257*pi*(n)) + 0.1*randn(size(n));
Xm = corrmtx(x,7,'modified');
pmusic(Xm,2)

ウィンドウ処理を使用した信号データ行列の効果の作成

同じ信号を使用し、ウィンドウ処理の入力引数を使用して、pmusic で 100 行 7 列のデータ行列を作成します。さらに、FFT の長さを 512 に指定します。

n = 0:699;
x = cos(0.257*pi*(n)) + 0.1*randn(size(n));
[PP,ff] = pmusic(x,2,512,[],7,0);
pmusic(x,2,512,[],7,0)

入力引数

x — 入力信号

ベクトル | 行列

ベクトルまたは行列として指定される入力信号。x がベクトルの場合は、信号の 1 つの観測値として扱われます。x が行列の場合、x の各行は信号の個別の観測値を表します。たとえば、x'*x が相関行列の推定となるような、配列処理の場合のように各行がセンサー配列の 1 つの出力である場合などです。

メモ

corrmtx の出力を使用して、x を生成できます。

複素数のサポート: あり

p — 部分空間の次元

実数の正の整数 | 2 要素ベクトル

部分空間の次元。実数の正の整数または 2 要素ベクトルとして指定します。p が実数の正の整数である場合は、部分空間の次元として扱われます。p が 2 要素ベクトルの場合は、p の 2 番目の要素は、信号の相関行列での推定固有値の最小値、λ_minを乗算したしきい値を表しています。しきい値 λ_min*p(2) より小さい固有値は、ノイズ部分空間に割り当てられます。この場合、p(1) では、信号部分空間の最大次元が指定されます。p の 2 番目のエントリにある追加のしきい値パラメーターを使用すると、ノイズ部分空間と信号部分空間の割り当てを柔軟にコントロールできるようになります。

メモ

pmusic の入力が実数正弦波である場合は、p の値は正弦波の数の 2 倍に設定します。入力が複素正弦波である場合は、p は正弦波数と等しく設定します。

複素数のサポート: あり

wi — 入力正規化周波数

ベクトル

入力正規化周波数。ベクトルとして指定します。

データ型: double

nfft — DFT 点の数

256 (既定値) | 整数 | []

正の整数として指定する DFT 点の数。nfft が空として指定されている場合、既定の nfft が使用されます。

fs — サンプル レート

1 (既定値) | 正のスカラー | []

サンプル レート。Hz 単位の正のスカラーとして指定します。fs を空ベクトル [] として指定すると、サンプル レートの既定値である 1 Hz が使用されます。

fi — 入力周波数

ベクトル

入力周波数。ベクトルとして指定します。ベクトルで指定された周波数において疑似スペクトルが計算されます。

nwin — 箱型ウィンドウの長さ

2*p(1) (既定値) | 非負の整数

箱型ウィンドウの長さ。非負の整数として指定します。

noverlap — オーバーラップするサンプル数

nwin-1 (既定値) | 非負の整数

オーバーラップするサンプル数。ウィンドウの長さより小さな非負の整数として指定します。

メモ

引数 nwin と noverlap は、構文で 'corr' を使用した場合は無視されます。

freqrange — 疑似スペクトル推定の周波数範囲

'half' | 'whole' | 'centered'

疑似スペクトル推定の周波数範囲。'half'、whole、または 'centered' のいずれかとして指定します。

• 'half' — 実数入力信号 x について半分のスペクトルを返します。nfft が偶数の場合、S は長さ nfft/2 + 1、計算区間は [0,π] です。nfft が奇数の場合、S の長さは (nfft + 1)/2、周波数範囲は [0,π) です。fs を指定すると、その区間は nfft が偶数の場合 [0,fs/2)、奇数の場合 [0,fs/2] になります。

• 'whole' — 実数または複素数の入力 x についてスペクトル全体を返します。この場合、S の長さは nfft、計算区間は [0, 2π) です。fs を指定する場合、周波数範囲は [0,fs) です。

• 'centered' — 実数または複素数の入力 x について中央に揃えたスペクトル全体を返します。この場合、S の長さは nfft で、計算区間は nfft が偶数の場合 (–π,π]、nfft が奇数の場合 (–π,π) になります。fs を指定すると、周波数範囲は nfft が偶数の場合 (–fs/2,fs/2]、奇数の場合 (–fs/2,fs/2) になります。

メモ

引数 freqrange または 'corr' は、入力引数リストの p の後で任意の位置に配置できます。

出力引数

S — 疑似スペクトル推定

ベクトル

疑似スペクトル推定。ベクトルとして返されます。疑似スペクトルは、入力データ x に関連した相関行列の固有ベクトルの推定を使用して計算されます。

wo — 出力正規化周波数

ベクトル

出力正規化周波数。ベクトルとして指定します。S と wo は同じ長さです。一般的に、FFT の長さと入力値 x により、計算される S の長さと対応する正規化周波数の範囲が決まります。この表は、最初の構文での S (と wo) の長さと、対応する正規化周波数の範囲を示しています。

長さ 256 (既定の設定)の FFT に対する S 特性

入力データ型	S と w0 の長さ	対応する正規化周波数の範囲
実数	129	[0, π]
複素数	256	[0, 2π)

nfft が指定されている場合、以下の表は S と wo の長さ、および wo に対する周波数範囲を示しています。

S と周波数ベクトルの特性

入力データ型	nfft 偶数または奇数	S と w の長さ	w の範囲
実数	偶数	(nfft/2 )+ 1	[0, π]
実数	奇数	(nfft + 1)/2	[0, π)
複素数	偶数または奇数	nfft	[0, 2π)

fo — 出力周波数

ベクトル

出力周波数。ベクトルとして返されます。fo に対する周波数範囲は、nfft、fs、および入力 x の値に依存します。S (と fo) の長さは、上の表のS と周波数ベクトルの特性と同じです。以下の表は、nfft と fs が指定されている場合の fo に対する周波数範囲を示しています。

fs を設定した S と周波数ベクトルの特性

入力データ型	nfft 偶数/奇数	f の範囲
実数	偶数	[0, fs/2]
実数	奇数	[0, fs/2)
複素数	偶数または奇数	[0, fs)

さらに、nwin と noverlap も指定されている場合は、相関行列固有値を推定するために使用される行列が定式化される前に、入力データ x が分割され、ウィンドウが適用されます。データのセグメンテーションは、nwin、noverlap、および x の形式に依存します。次の表に、結果として得られるウィンドウを適用したセグメントに関するコメントを記載します。

x と nwin に依存するウィンドウを適用したデータ

x の形式	nwin の形式	ウィンドウが適用されたデータ
データ ベクトル	スカラー	長さが nwin。
データ ベクトル	係数ベクトル	長さが length(nwin)。
データ行列	スカラー	データには、ウィンドウが適用されません。
データ行列	係数ベクトル	length(nwin) は、x の列の長さと同じでなければなりません。また、noverlap は使用されません。

この構文の関連情報は、入力データと構文に依存した固有ベクトルの長さを参照してください。

v — ノイズの固有ベクトル

行列

ノイズの固有ベクトル。行列として返されます。v の列は、次元が size(v,2) のノイズ部分空間をにまたがっています。信号部分空間の次元は、size(v,1)-size(v,2) です。

e — 推定固有値

ベクトル

相関行列の推定固有値。ベクトルとして返されます。

ヒント

疑似スペクトルの推定の過程で、関数 pmusic は、信号の相関行列の推定された固有ベクトル v_j および固有値 λ_j からノイズと信号のそれぞれの部分空間を計算します。これらの固有値の中の最小のものが、しきい値パラメーター p(2) と共に使用され、場合によってはノイズ部分空間の次元に影響を与えます。

関数 pmusic で計算される固有ベクトルの長さ n は、信号部分空間の次元とノイズ部分空間の次元との和です。この固有ベクトルの長さは、入力 (信号データまたは相関行列) と使用する構文に依存します。

次の表は、入力引数への固有ベクトル長の依存度をまとめたものです。

入力データと構文に依存した固有ベクトルの長さ

入力データ x の形式	構文に関するコメント	固有ベクトルの長さ n
行ベクトルまたは列ベクトル	nwin は、スカラー整数として指定します。	nwin
行ベクトルまたは列ベクトル	nwin は、ベクトルとして指定します。	length(nwin)
行ベクトルまたは列ベクトル	nwin は指定しません。	2 × p(1)
l 行 m 列の行列	nwin がスカラーとして指定された場合は、使用されません。nwin がベクトルとして指定された場合、length(nwin) は m と等しくなければなりません。	m
m 行 m 列の非負定値行列	'corr' が指定され、nwin は使用されません。	m

p(2) > 1 を有効にするには、nwin > p(1) または length(nwin) > p(1) のいずれかを設定しなければなりません。

アルゴリズム

多重信号分類 (MUSIC) アルゴリズムでは、シュミットの固有空間解析法を使用して、信号または相関行列から疑似スペクトルが推定されます [1]。このアルゴリズムでは、信号の周波数成分を推定するために、信号の相関行列の固有空間解析が行われます。このアルゴリズムは、正弦波に加法性ホワイト ガウス ノイズを加算した信号に対して、特に適しています。相関行列を指定しない場合は、信号の相関行列の固有値と固有ベクトルが推定されます。

MUSIC 疑似スペクトル推定の計算に使用されるアルゴリズムは、信号部分空間とノイズ部分空間の相対サイズを考慮します。

• 信号部分空間の次元が、ノイズ部分空間の次元より小さい場合、MUSIC 疑似スペクトル推定は次のように与えられます。

  $$P_{\text{MUSIC}}(f)=\frac{1}{{\displaystyle \sum _{k=1}^p\left(1-{\left|v_k^{H}e\left(f\right)\right|}^2\right)}},$$

  ここで、p は信号部分空間の次元、v_k は相関行列の k 番目の固有ベクトルです。和において使用される固有ベクトルは、最大の固有値に対応し、信号部分空間を範囲とします。

• ノイズ部分空間の次元が、信号部分空間の次元より小さい場合、MUSIC 疑似スペクトル推定は次のように与えられます。

  $$P_{\text{MUSIC}}(f)=\frac{1}{e^H(f)\left({\displaystyle \sum _{k=p+1}^N{v}_k{v}_k^H}\right)e(f)}=\frac{1}{{\displaystyle \sum _{k=p+1}^N{\left|v_k^{H}e(f)\right|}^2}},$$

  ここで、N は信号部分空間の次元にノイズ部分空間の次元を加えた値、v_k は相関行列の k 番目の固有ベクトルです。和において使用される固有ベクトルは、最小の固有値に対応し、ノイズ部分空間を範囲とします。

各場合において、ベクトル e(f) は、複素指数で構成されるため、内積 v_k^He(f) は、フーリエ変換になります。これは、PSD 推定の計算で使用されます。FFT は、各 v_k に対して計算され、振幅の 2 乗の和が計算されます。

参照

[1] Marple, S. Lawrence. Digital Spectral Analysis. Englewood Cliffs, NJ: Prentice-Hall, 1987, pp. 373–378.

[2] Schmidt, R. O. “Multiple Emitter Location and Signal Parameter Estimation.” IEEE^® Transactions on Antennas and Propagation. Vol. AP-34, March, 1986, pp. 276–280.

[3] Stoica, Petre, and Randolph L. Moses. Spectral Analysis of Signals. Upper Saddle River, NJ: Prentice Hall, 2005.

拡張機能

C/C++ コード生成

MATLAB® Coder™ を使用して C および C++ コードを生成します。

使用上の注意および制限:

コード生成時に nfft または nwin が可変サイズである場合、実行時にスカラーまたは空の配列に変換してはいけません。

バージョン履歴

R2006a より前に導入

参考

corrmtx | pburg | peig | periodogram | pmtm | prony | pwelch | rooteig | rootmusic