ドキュメンテーション

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

tfridge

説明

fridge = tfridge(tfm,f) は、時間-周波数行列 tfm と周波数ベクトル f から最大エネルギーの時間-周波数リッジを抽出し、時間依存周波数 fridge を出力します。

[fridge,iridge] = tfridge(tfm,f) は最大エネルギー リッジに対応する行インデックス ベクトルも返します。

[fridge,iridge,lridge] = tfridge(tfm,f) は線形インデックス lridge も返します。この場合、tfm(lridge) は最大エネルギー リッジに沿った tfm の値になります。

[___] = tfridge(tfm,f,penalty) は、周波数ビン間の平方距離を penalty によってスケーリングすることにより、周波数の変更にペナルティを課します。

[___] = tfridge(___,'NumRidges',nr) は、最大エネルギーの nr 時間-周波数リッジを抽出します。この構文は、前の構文の入力引数の任意の組み合わせを受け入れます。

[___] = tfridge(___,'NumRidges',nr,'NumFrequencyBins',nbins) は、複数のリッジを抽出するとき、tfm から除去されるリッジの周囲の周波数ビンの数を指定します。

すべて折りたたむ

鋭いリッジをもつ時間-周波数行列に類似する行列を作成します。行列を 3 次元で可視化します。

t = 0:0.05:10;
f = 0:0.2:8;
rv = 1;

[F,T] = ndgrid(f,t);

S = zeros(size(T));
S(abs((F-4)-cos((T-6).^2))<0.1) = rv;

mesh(T,F,S)
view(-30,60)

行列にノイズを追加してプロットを再表示します。

S = S+rand(size(S))/10;

mesh(T,F,S)
view(-30,60)
xlabel('Time')
ylabel('Frequency')

リッジを抽出して結果をプロットします。

[fridge,~,lridge] = tfridge(S,f);

rvals = S(lridge);

hold on
plot3(t,fridge,rvals,'k','linewidth',4)
hold off

3 kHz で 1 秒間サンプリングされる信号を生成します。信号は 2 つのトーンと 2 次チャープで構成されます。

  • 最初のトーンは周波数 1000 Hz、単位振幅です。

  • 2 番目のトーンは周波数 1200 Hz、単位振幅です。

  • チャープの初期周波数は 500 Hz で、サンプリングの最後には 750 Hz に到達します。振幅は 6 です。

fs = 3000;
t = 0:1/fs:1-1/fs;

x1 = 6*chirp(t,fs/6,t(end),fs/4,'quadratic');
x2 = sin(2*pi*fs/3*t);
x3 = sin(2*pi*fs/2.5*t);

x = x1+x2+x3;

信号のフーリエ シンクロスクイーズド変換を計算し、表示します。

[sst,f] = fsst(x,fs);
mx = max(abs(sst(:)))*ones(size(t));

mesh(t,f,abs(sst))
view(2)

最も高いエネルギーの信号成分を 2 つ抽出し、プロットします。周波数の変更に対するペナルティを設定しません。

penval = 0;

fridge = tfridge(sst,f,penval,'NumRidges',2);

hold on
plot3(t,fridge,mx,'w','linewidth',5)
hold off

2 つのトーンは同じ振幅をもち、アルゴリズムはそれらの間をジャンプします。周波数の変更に対してペナルティを 1 に設定します。

penval = 1;

fridge = tfridge(sst,f,penval,'NumRidges',2);

mesh(t,f,abs(sst))
view(2)
xlabel('Time (s)')
ylabel('Frequency (Hz)')

hold on
plot3(t,fridge,mx,'w','linewidth',5)
hold off

比較のために、ペナルティを高い値に設定します。チャープの周波数が一定でないため、ペナルティが課されます。

penval = 1000;

fridge = tfridge(sst,f,penval,'NumRidges',2);

mesh(t,f,abs(sst))
view(2)
xlabel('Time (s)')
ylabel('Frequency (Hz)')

hold on
plot3(t,fridge,mx,'w','linewidth',5)
hold off

2 つの 2 次チャープで構成された信号を生成します。信号は 1 kHz で 3 秒間サンプリングされます。チャープでは、瞬時周波数がサンプリング間隔の中間点について対称です。一方のチャープは凹であり、他方のチャープは凸です。凹チャープは、凸チャープの 2 倍の振幅があります。

fs = 1e3;
t = 0:1/fs:3;

x =   chirp(t-1.5,100,1.1,200,'quadratic',[],'convex');
y = 2*chirp(t-1.5,300,1.1,400,'quadratic',[],'concave');

% To hear, type soundsc(x+y,fs)

信号のフーリエ シンクロスクイーズド変換を計算し、表示します。

sig = x+y;

[sst,f,t] = fsst(sig,fs);

fsst(sig,fs,'yaxis')

最大のエネルギーをもつ 2 つの時間-周波数リッジを抽出します。周波数の変更に対してペナルティ 1 を指定します。2 番目のリッジを抽出する前に、最大エネルギーをもつリッジ周囲にある 1 周波数ビンを削除します。リッジをプロットします。

nfb = 1;
[fr,ir] = tfridge(sst,f,1,'NumRidges',2,'NumFrequencyBins',nfb);

plot(t,fr)

1 ビンでは十分ではありません。関数が、最初のリッジの勾配に、2 番目のリッジの一部が存在することを検出しました。削除するビン数を 50 に増やして、計算を繰り返します。

nfb = 50;
[fr,ir] = tfridge(sst,f,1,'NumRidges',2,'NumFrequencyBins',nfb);

plot(t,fr)

削除するビン数が多すぎると、低いエネルギーのリッジが歪みます。数を 15 に減らして、計算を繰り返します。

nfb = 15;
[fr,ir] = tfridge(sst,f,1,'NumRidges',2,'NumFrequencyBins',nfb);

plot(t,fr)

2 つのリッジに対応する変換の逆変換を行います。リッジを追加して、信号を再構成します。再構成した信号とチャープの差をプロットします。

itr = ifsst(sst,[],ir,'NumFrequencyBins',nfb);
xrec = sum(itr');

plot(t,xrec-(x+y))
ylim([-.1 .1])

% To hear, type soundsc(xrec,fs)

一致度は、ほとんどの時間で良好ですが、周波数が最も急激に変化する両端では劣化しています。

入力引数

すべて折りたたむ

時間-周波数行列。行列として指定します。

例: fsst(cos(pi/4*(0:159))) は、正弦波のシンクロスクイーズド変換を指定します。

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

サンプリング周波数。ベクトルとして指定します。f の長さは、tfm の行数と等しくなければなりません。

データ型: single | double

周波数の変更に対するペナルティ。非負の実数スカラーとして指定します。

データ型: single | double

抽出する時間-周波数リッジの数。'NumRidges' と正の整数スカラーとで構成されるコンマ区切りのペアとして指定します。この名前と値のペアは、入力引数リストで任意の場所に tfm に続けて指定できます。

nr が 1 より大きい場合、tfridge は以下を行います。

  1. 最大エネルギーの時間-周波数リッジを抽出

  2. tfm から、抽出されたリッジおよびそのリッジのいずれかの側の隣接する周波数のビン nbins 個に含まれるエネルギーを削除

  3. 変更された tfm 内の最大エネルギー リッジを抽出

  4. nr 個のリッジが抽出されるまで繰り返す

データ型: single | double

複数のリッジを抽出するときに削除するビンの数。'NumFrequencyBins' および正の整数スカラーで構成されるコンマ区切りのペアとして指定します。nbins は、サンプリング周波数の 1/4 より小さくなければなりません。片側のビンの数が nbins より少ない周波数エッジの近傍のインデックスは、より少ないビンの数を使用して再構成されます。

データ型: single | double

出力引数

すべて折りたたむ

時間-周波数リッジ。nr 列の行列として返されます。fridge の行数は、tfm の列数と同じになります。最初の列は、最大エネルギーのリッジに対応する周波数を含みます。後続の列は、エネルギーの高い順に他のリッジの周波数を含みます。

リッジの行インデックス。nr 列の行列として返されます。iridge の行数は、tfm の列数と同じになります。最初の列は、最大エネルギーのリッジに対応するインデックスを含みます。後続の列は、エネルギーの高い順に他のリッジのインデックスを含みます。

リッジの線形インデックス。nr 列の行列として返されます。lridge は、tfm(lridge) がリッジに沿った tfm の振幅になるように定義されます。lridge の行数は、tfm の列数と同じになります。最初の列は、最大エネルギーのリッジに対応するインデックスを含みます。後続の列は、エネルギーの高い順に他のリッジのインデックスを含みます。

例: 関数 lridge は、関数 sub2ind(size(tfm),iridge,repmat((1:size(tfm,2))',1,nr)) と等価です。

アルゴリズム

関数は、ペナルティが課された順方向/逆方向のグリーディ アルゴリズムを使用して、最大エネルギー リッジを時間-周波数行列から抽出します。アルゴリズムは、各時間点で –ln A を最小化することによって最大時間-周波数リッジを検索します。ここで、A は行列の絶対値です。–ln A を最小化することは A の値を最大化することと等価です。アルゴリズムは、オプションで、周波数ビン間の距離に比例するペナルティをもつ周波数のジャンプを制約します。

次の例は、周波数ビン間の距離の 2 倍のペナルティを使用した時間-周波数リッジのアルゴリズムを示しています。具体的には、要素 (j,k)(m,n) 間の距離は、(j-m)2 として定義されます。時間-周波数の行列には、3 つの周波数ビンと 3 つのタイム ステップがあります。行列の列はタイム ステップに対応し、行列の行は周波数ビンに対応します。2 つ目の行の値は正弦波を表します。

  1. 行列は次のとおりとします。

    1   4   4
    2   2   2
    5   5   4

  2. (1,2) 要素の値の更新は次のとおりです。

    1. 最初の時間点の値はそのまま変更しません。行列の (1,2) 要素をもつアルゴリズムを開始します。これは、2 つ目の時間点の最初の周波数ビンを表します。ビンの値は 4 です。(1,2) 要素からそれらの距離を基にして最初の列の値にペナルティを課します。最初の列にペナルティを適用して次を作成します。

      original value + penalty × distance
      
      1 + 2 × 0 =  1
      2 + 2 × 1 =  4
      5 + 2 × 4 = 13
      
       1   4
       4   2
      13   5
      最初の列の最小値は 1 です。これはビン 1 にあります。

    2. 列 1 の最小値を現在のビンの値 4 に追加します。(1,2) の更新された値は 5 になります。これはビン 1 から派生しています。

  3. 列 2 の残りの要素に対して値を次のように更新します。

    ステップ 2a と同じプロセスを使用してペナルティ係数をもつ、元の列 1 の値を再計算します。ステップ 2b と同じプロセスを使用して残りの 2 つ目の列の値を取得します。たとえば、(2,2) 要素を更新する場合 (ビンの値 2 があります)、列にペナルティを適用すると以下が得られます。

    original value + penalty × distance
    
    1 + 2 × 1 =  3
    2 + 2 × 0 =  2
    5 + 2 × 1 =  7
    
    最小値 2 を現在のビン値に追加します。(2,2) の更新された値が 4 になります。(3,2) 要素を更新した後、行列は以下のようになります。
    1   5(1)  4
    2   4(2)  2
    5   9(2)  4
    2 つ目の列のみ更新されています。添字は値の派生元の列にあるビンのインデックスを示します。

  4. 3 つ目の列に対してステップ 2 を繰り返します。ただし、ペナルティは更新された 2 つ目の列に適用されます。たとえば、(1,3) 要素を更新すると、ペナルティは次のようになります。

    5 + 2 × 0 =  5
    4 + 2 × 1 =  6
    9 + 2 × 4 = 17
    
    最小値 5 (最初のビンにあります) が (1,3) ビン値に追加されます。3 つ目の列の値すべてを更新した後、最後の行列は次のようになります。
    1   5(1)   9(1)
    2   4(2)   6(2)
    5   9(2)  10(2)

  5. 行列の最後の列から始めて、最小値を見つけます。現在のビンから前の時間点のそのビンの原点に移行することによって行列を通じて時間順に戻ります。ビンのインデックスを継続して追跡します。これはリッジを構成するパスを形成します。アルゴリズムは最小値をもつビンの代わりに元のビンを使用して遷移を平滑化します。この例の場合、リッジのインデックスは、222 です。これは、ステップ 1 で示す行列の行 2 で正弦波のエネルギー パスと一致します。

複数のリッジを抽出している場合、アルゴリズムは時間-周波数の行列から最初のリッジを削除して、プロセスを繰り返します。

R2016b で導入