このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

xwvd

交差 Wigner-Ville 分布と交差平滑化疑似 Wigner-Ville 分布

構文

d = xwvd(x,y)

d = xwvd(x,y,fs)

d = xwvd(x,y,ts)

d = xwvd(___,"smoothedPseudo")

d = xwvd(___,"smoothedPseudo",twin,fwin)

d = xwvd(___,"smoothedPseudo",NumFrequencyPoints=nf)

d = xwvd(___,MinThreshold=thresh)

[d,f,t] = xwvd(___)

xwvd(___)

説明

d = xwvd(x,y) は、x と y の交差 Wigner-Ville 分布を返します。

d = xwvd(x,y,fs) は、x と y がレート fs でサンプリングされたときの交差 Wigner-Ville 分布を返します。

d = xwvd(x,y,ts) は、x と y がサンプル間の時間間隔 ts でサンプリングされたときの交差 Wigner-Ville 分布を返します。

d = xwvd(___,"smoothedPseudo") は、x と y の交差平滑化疑似 Wigner-Ville 分布を返します。この関数は、入力信号の長さを使用して、時間と周波数の平滑化に使用されるウィンドウの長さを選択します。この構文には、前の構文の入力引数を任意に組み合わせて含めることができます。

d = xwvd(___,"smoothedPseudo",twin,fwin) は、平滑化に使用される時間ウィンドウ twin と周波数ウィンドウ fwin を指定します。時間または周波数の平滑化に既定のウィンドウを使用するには、対応する引数を空 ([]) に指定します。

d = xwvd(___,"smoothedPseudo",NumFrequencyPoints=nf) は、nf 個の周波数点を使用して交差平滑化疑似 Wigner-Ville 分布を計算します。この構文では、twin と fwin を指定できます。また、省略することもできます。

d = xwvd(___,MinThreshold=thresh) は、振幅が thresh より小さい d の要素をゼロに設定します。この構文は、交差 Wigner-Ville 分布と交差平滑化疑似 Wigner-Ville 分布の両方に適用されます。

[d,f,t] = xwvd(___) は、周波数のベクトル f と、d が計算される時間のベクトル t も返します。

出力引数なしで xwvd(___) を使用すると、現在の Figure に交差 Wigner-Ville 分布または交差平滑化疑似 Wigner-Ville 分布の実数部がプロットされます。

例

信号の交差 Wigner-Ville 分布

1 kHz で 1 秒間サンプリングされ、ホワイト ノイズに組み込まれた 2 つの信号を生成します。片方の信号は、周波数 150 Hz の正弦波です。もう片方の信号は、200 Hz と 400 Hz の間で周波数が正弦関数的に変化するチャープです。ノイズの分散は $0.1^2$ です。

fs = 1000;
t = (0:1/fs:1)';

x = cos(2*pi*t*150) + 0.1*randn(size(t));
y = vco(cos(3*pi*t),[200 400],fs) + 0.1*randn(size(t));

信号の和の Wigner-Ville 分布を計算します。

wvd(x+y,fs)

信号の交差 Wigner-Ville 分布を計算してプロットします。交差分布は、Wigner-Ville 分布の交差項に対応します。

xwvd(x,y,fs)

チャープの交差 Wigner-Ville 分布

2 つのチャープで構成された 2 チャネル信号を生成します。信号は、3 kHz で 1 秒間サンプリングされます。最初のチャープの初期周波数は 400 Hz で、サンプリングの最後には 800 Hz に到達します。2 番目のチャープは 500 Hz から開始し、最後には 1000 Hz に到達します。2 番目のチャープの振幅は最初のチャープ 2 倍です。

fs = 3000;
t = (0:1/fs:1-1/fs)';

x1 = chirp(t,1400,t(end),800);
x2 = 2*chirp(t,200,t(end),1000);

信号を timetable として保存します。2 つのチャネルの交差 Wigner-Ville 分布を計算してプロットします。

xt = timetable(seconds(t),x1,x2);

xwvd(xt(:,1),xt(:,2))

交差 Wigner-Ville 分布を使用した瞬時周波数の推定

既知の基準信号と交差 Wigner-Ville 分布を使用して、信号の瞬時周波数を計算します。

1 kHz で 1 秒間サンプリングされたガウス原子で構成される基準信号を作成します。ガウス原子は、ガウスにより変調された正弦波です。50 Hz の正弦波周波数を指定します。ガウスの中心は 64 ミリ秒で、$0.{01}^2$ の分散があります。

fs = 1e3;
t = (0:1/fs:1-1/fs)';

mu = 0.064;
sigma = 0.01;
fsin = 50;

xr = exp(-(t-mu).^2/(2*sigma^2)).*sin(2*pi*fsin*t);

チャープで構成される、解析用の "未知の" 信号を作成します。この信号は 0.4 秒で突然開始し、0.5 秒後に突然終了します。この間に、チャープの周波数は 400 Hz から 100 Hz に線形的に減少します。

f0 = 400;
f1 = 100;

xa = zeros(size(t));
xa(t>0.4 & t<=0.9) = chirp((0:1/fs:0.5-1/fs)',f0,0.5,f1);

未知の信号と基準信号の和で構成される 2 成分信号を作成します。結果の平滑化疑似 Wigner-Ville 分布は、"理想的な" 時間-周波数表現を提供します。

平滑化疑似 Wigner-Ville 分布を計算して表示します。

w = wvd(xa+xr,fs,"smoothedPseudo");

wvd(xa+xr,fs,"smoothedPseudo")

未知の信号と基準信号の交差 Wigner-Ville 分布を計算します。分布の絶対値を取り、振幅が 10 より小さい要素をゼロに設定します。交差 Wigner-Ville 分布は、2 成分信号の交差項に等しくなります。

交差 Wigner-Ville 分布の実数部をプロットします。

[c,fc,tc] = xwvd(xa,xr,fs);
c = abs(c);
c(c<10) = 0;

xwvd(xa,xr,fs)

理想的な時間-周波数表現を交差 Wigner-Ville 分布に追加することで、Wigner-Ville 交差項を強調します。Wigner-Ville 分布の交差項は、基準信号と未知の信号の中間に発生します。

d = w + c;

d = abs(real(d));

imagesc(tc,fc,d)
axis xy
colorbar

交差項に対応する高エネルギー リッジを特定してプロットします。リッジを分離するには、交差分布に非ゼロのエネルギーがある時間値を検索します。

ff = tfridge(c,fc);

tv = sum(c)>0;

ff = ff(tv);
tc = tc(tv);

hold on
plot(tc,ff,"r--",linewidth=2)
hold off

リッジと参照関数を使用して、未知の信号の瞬時周波数を再構成します。瞬時周波数を時間の関数としてプロットします。

tEst = 2*tc - mu;
fEst = 2*ff - fsin;

plot(tEst,fEst)

入力引数

x, y — 入力信号

ベクトル | timetable

入力信号。ベクトルまたはそれぞれが単一のベクトル変数を含む MATLAB^® の timetable として指定します。x と y は、両方ともベクトルであるか両方とも timetable である必要があり、同じ長さを持たなければなりません。

• x と y が timetable の場合、増加する有限の行時間を含んでいなければなりません。

• timetable が欠損している場合や時間点が重複している場合、欠損または重複する時間および非等間隔の時間をもつ timetable の整理のヒントを使用して修正できます。

入力信号の長さが奇数の場合は、関数によりゼロが追加され、長さが偶数になります。

例: cos(pi/8*(0:159))'+randn(160,1)/10 は、ホワイト ノイズに含まれる正弦波を指定します。

例: timetable(seconds(0:5)',rand(6,1)) は 1 Hz で 4 秒間サンプリングされた確率変数を指定します。

データ型: single | double
複素数のサポート: あり

fs — サンプル レート

2*pi (既定値) | 正の数値スカラー

サンプル レート。正の数値スカラーとして指定します。

ts — サンプル時間

duration スカラー

サンプル時間。duration スカラーで指定します。

twin, fwin — 時間ウィンドウと周波数ウィンドウ

長さが奇数のベクトル

平滑化に使用される時間ウィンドウと周波数ウィンドウ。長さが奇数のベクトルとして指定します。既定では、xwvd は形状係数β = 20 のカイザー ウィンドウを使用します。

• twin の既定の長さは、round(length(x)/10) 以上の最も小さい奇数の整数です。

• fwin の既定の長さは、nf/4 以上の最も小さい奇数の整数です。

各ウィンドウの長さは、2*ceil(length(x)/2) 以下でなければなりません。

例: kaiser(65,0.5) は、0.5 の形状係数を持つ 65 サンプルのカイザー ウィンドウを指定します。

nf — 周波数点の数

2*ceil(length(x)/2) (既定値) | 整数

周波数点の数。整数として指定します。この引数は、周波数のオーバーサンプリングの程度を制御します。周波数点の数は、(length(fwin)+1)/2 以上でなければならず、既定よりも大きくすることはできません。

thresh — 最小非ゼロ値

-Inf (既定値) | 実数スカラー

最小非ゼロ値。実数スカラーとして指定します。関数は、振幅が thresh より小さい d の要素をゼロに設定します。

出力引数

d — 交差 Wigner-Ville 分布

行列

交差 Wigner-Ville 分布。行列として返されます。時間は d の列方向に、周波数は行方向に下がって増加します。行列のサイズは N_f × N_t です。ここで、N_f は f の長さ、N_t は t の長さです。

f — 周波数

ベクトル

ベクトルとして返される周波数。

• 入力に時間情報がある場合、f は Hz 単位で表される周波数を含みます。

• 入力に時間情報がない場合、f はラジアン/サンプル単位で表される正規化周波数を含みます。

t — 時点

ベクトル

時点。ベクトルとして返されます。

• 入力に時間情報がある場合、t は秒単位で表される時間値を含みます。

• 入力に時間情報がない場合、t はサンプル数を含みます。

時間点の数は、4*ceil(length(x)/2) に固定されます。

詳細

交差 Wigner-Ville 分布

連続信号 x(t) および y(t) では、"交差 Wigner-Ville 分布" は次のように定義されます。

$${\mathrm{XWVD}}_{x,y}(t,f)={\displaystyle {\int }_{-\infty }^{\infty }x\left(t+\frac{\tau }{2}\right)y^{*}\left(t-\frac{\tau }{2}\right)e^{-j2\pi f\tau }d\tau .}$$

N 個のサンプルをもつ離散信号では、分布は次のようになります。

$${\text{XWVD}}_{x,y}(n,k)={\displaystyle \sum _{m=-N}^{N}x(n+m/2)y^{*}(n-m/2)e^{-j2\pi km/N}}.$$

奇数値の m では、この定義は半整数サンプル値での信号の評価を必要とします。これにより内挿が必要となるため、エイリアシングを避けるために離散フーリエ変換をゼロ パディングする必要があります。

交差 Wigner-Ville 分布には、その解釈を複雑にしがちな干渉項が含まれています。分布をシャープにするために、ローパス ウィンドウで定義をフィルター処理できます。交差平滑化疑似 Wigner-Ville 分布は、時間と周波数の平滑化に個別のウィンドウを使用します。

$${\text{XSPWVD}}_{x,y}^{g,H}(t,f)={\displaystyle {\int }_{-\infty }^{\infty }g(t)H(f)x\left(t+\frac{\tau }{2}\right)y^{*}\left(t-\frac{\tau }{2}\right)e^{-j2\pi f\tau }d\tau .}$$

参照

[1] Cohen, Leon. Time-Frequency Analysis: Theory and Applications. Englewood Cliffs, NJ: Prentice-Hall, 1995.

[2] Mallat, Stéphane. A Wavelet Tour of Signal Processing. Second Edition. San Diego, CA: Academic Press, 1999.

[3] Malnar, Damir, Victor Sucic, and Boualem Boashash. "A cross-terms geometry based method for components instantaneous frequency estimation using the cross Wigner-Ville distribution." In 11th International Conference on Information Sciences, Signal Processing and their Applications (ISSPA), pp. 1217–1222. Montréal: IEEE^®, 2012.

拡張機能

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

バージョン履歴

R2018b で導入

R2023a: コード生成向けの timetable のサポート

関数 xwvd は、コード生成向けに timetable 入力をサポートします。

参考

関数

• cpsd | fsst | pspectrum | xspectrogram | wvd

トピック

• 時間-周波数ギャラリー