このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。
tfridge
時間-周波数リッジ
構文
説明
例
複素チャープの瞬時周波数
この例では、フーリエ シンクロスクイーズド変換を使用して、信号の瞬時周波数を計算する方法を示します。
正弦関数的に変化する周波数成分をもつチャープを生成します。信号はホワイト ガウス ノイズに組み込まれ、3 kHz で 1 秒間サンプリングされます。
fs = 3000; t = 0:1/fs:1-1/fs; x = exp(2j*pi*100*cos(2*pi*2*t)) + randn(size(t))/100;
信号のフーリエ シンクロスクイーズド変換を計算し、プロットします。時間が x 軸に、周波数が y 軸に表示されます。
fsst(x,fs,'yaxis')
信号の瞬時周波数を求めるために、フーリエ シンクロスクイーズド変換の最大エネルギーの時間-周波数リッジを抽出します。
[sst,f,tfs] = fsst(x,fs); fridge = tfridge(sst,f);
リッジを変換プロットに重ね合わせます。時間をミリ秒に、周波数を kHz に変換します。
hold on plot(t*1000,fridge/1000,'r') hold off
実信号の場合、関数 instfreq を使用して瞬時周波数をより簡単に検出できます。たとえば、解析信号を計算し、その位相を微分することで複素チャープの実数部の瞬時周波数を表示します。
ax = real(x); instfreq(ax,fs,'Method','hilbert')
ノイズを含む信号のリッジを検出
鋭いリッジをもつ時間-周波数行列に類似する行列を作成します。行列を 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)一致度は、ほとんどの時間で良好ですが、周波数が最も急激に変化する両端では劣化しています。
入力引数
出力引数
アルゴリズム
関数は、ペナルティが課された順方向/逆方向の貪欲法を使用して、最大エネルギー リッジを時間-周波数行列から抽出します。アルゴリズムは、各時間点で –ln A を最小化することによって最大時間-周波数リッジを検索します。ここで、A は行列の絶対値です。–ln A を最小化することは A の値を最大化することと等価です。アルゴリズムは、オプションで、周波数ビン間の距離に比例するペナルティをもつ周波数のジャンプを制約します。
次の例は、周波数ビン間の距離の 2 倍のペナルティを使用した時間-周波数リッジのアルゴリズムを示しています。具体的には、要素 (j,k) と (m,n) 間の距離は、(j-m)2 として定義されます。時間-周波数の行列には、3 つの周波数ビンと 3 つのタイム ステップがあります。行列の列はタイム ステップに対応し、行列の行は周波数ビンに対応します。2 つ目の行の値は正弦波を表します。
行列は次のとおりとします。
1 4 4 2 2 2 5 5 4
(1,2) 要素の値の更新は次のとおりです。
最初の時間点の値はそのまま変更しません。行列の (1,2) 要素をもつアルゴリズムを開始します。これは、2 つ目の時間点の最初の周波数ビンを表します。ビンの値は 4 です。(1,2) 要素からそれらの距離を基にして最初の列の値にペナルティを課します。最初の列にペナルティを適用して次を作成します。
original value + penalty × distance 1 + 2 × 0 = 1 2 + 2 × 1 = 4 5 + 2 × 4 = 13
最初の列の最小値は 1 です。これはビン 1 にあります。1 4 4 2 13 5
列 1 の最小値を現在のビンの値 4 に追加します。(1,2) の更新された値は 5 になります。これはビン 1 から派生しています。
列 2 の残りの要素に対して値を次のように更新します。
ステップ 2a と同じプロセスを使用してペナルティ係数をもつ、元の列 1 の値を再計算します。ステップ 2b と同じプロセスを使用して残りの 2 つ目の列の値を取得します。たとえば、(2,2) 要素を更新する場合 (ビンの値 2 があります)、列にペナルティを適用すると以下が得られます。
最小値 2 を現在のビン値に追加します。(2,2) の更新された値が 4 になります。(3,2) 要素を更新した後、行列は次のようになります。original value + penalty × distance 1 + 2 × 1 = 3 2 + 2 × 0 = 2 5 + 2 × 1 = 7
2 つ目の列のみ更新されています。添字は値の派生元の列にあるビンのインデックスを示します。1 5(1) 4 2 4(2) 2 5 9(2) 4
3 つ目の列に対して手順 2 を繰り返します。ただし、ペナルティは更新された 2 つ目の列に適用されます。たとえば、(1,3) 要素を更新すると、ペナルティは次のようになります。
最小値 5 (最初のビンにあります) が (1,3) ビン値に追加されます。3 つ目の列の値すべてを更新した後、最後の行列は次のようになります。5 + 2 × 0 = 5 4 + 2 × 1 = 6 9 + 2 × 4 = 17
1 5(1) 9(1) 2 4(2) 6(2) 5 9(2) 10(2)
行列の最後の列から始めて、最小値を見つけます。現在のビンから前の時間点のそのビンの原点に移行することによって行列を通じて時間順に戻ります。ビンのインデックスを継続して追跡します。これはリッジを構成するパスを形成します。アルゴリズムは最小値をもつビンの代わりに元のビンを使用して遷移を平滑化します。この例の場合、リッジのインデックスは、
2、2、2です。これは、ステップ 1 で示す行列の行 2 で正弦波のエネルギー パスと一致します。
複数のリッジを抽出している場合、アルゴリズムは時間-周波数の行列から最初のリッジを削除して、プロセスを繰り返します。
拡張機能
バージョン履歴
R2016b で導入
You can also select a web site from the following list:
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)