Main Content

このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。

時間-周波数ギャラリー

このギャラリーでは、Signal Processing Toolbox™ および Wavelet Toolbox™ で使用可能な時間-周波数解析関数の概要を示します。説明と使用例で、信号解析に使用できるさまざまな方法を紹介しています。

メソッド特徴可逆性

短時間フーリエ変換 (スペクトログラム)

  • 短時間フーリエ変換 (STFT) の時間-周波数分解能は固定されます。

  • スペクトログラムは、STFT の振幅の 2 乗です。

  • stft:あり

  • spectrogram:なし

例:クジラの歌

連続ウェーブレット変換 (スカログラム)

  • 連続ウェーブレット変換 (CWT) の時間-周波数分解能は可変です。

  • CWT では、時間シフトと時間スケーリングが維持されます。

あり

例:ECG 信号

Wigner-Ville 分布

  • Wigner-Ville 分布 (WVD) は常に実数です。

  • 時間周辺密度と周波数周辺密度はそれぞれ、瞬時電力とスペクトル エネルギー密度に対応しています。

  • WVD の時間分解能は、入力サンプルの数に等しくなります。

なし

例:耳音響放射

再割り当てとシンクロスクイージング

  • 再割り当てにより、スペクトルの位置推定が鮮明になります。

  • シンクロスクイージングは、瞬時周波数の曲線の周囲の時間-周波数マップを "凝縮" します。

  • どちらの方法も、時間-周波数リッジの追跡と抽出に特に適しています。

  • pspectrum:なし

  • fsst, wsst:あり

例:反響定位パルス

定 Q ガボール変換

  • Q ガボール変換 (CQT) は、可変サイズのウィンドウで時間-周波数平面をタイリングします。

  • ウィンドウには、適応性のある帯域幅およびサンプリング密度があります。

  • すべてのウィンドウの帯域幅に対する中心周波数の比率 (Q 係数) は一定です。

あり

例:ロック音楽

経験的モード分解とヒルベルト・ファン変換

  • 経験的モード分解 (EMD) は、信号を固有モード関数に分解します。

  • ヒルベルト・ファン変換 (HHT) は、各経験的モードの瞬時周波数を計算します。

なし

例:ベアリング振動

短時間フーリエ変換 (スペクトログラム)

説明

  • "短時間フーリエ変換" は、非定常多成分信号の解析に役立つ線形時間-周波数表現です。

  • 短時間フーリエ変換は可逆変換です。

  • スペクトログラムは、STFT の振幅の 2 乗です。

  • 2 つの信号のクロス スペクトログラムを計算して、時間-周波数空間での類似度を調べることができます。

  • 信号の "パーシステンス スペクトル" は、与えられた周波数が信号内に存在する時間の割合を示す時間-周波数領域の表示です。パーシステンス スペクトルはパワー周波数空間のヒストグラムです。信号が変化する中で特定の周波数が信号内に存在する時間が長ければ長いほど、その時間の割合は大きくなるため、表示内の色が明るく ("熱く") なります。

考えられる用途

この時間-周波数手法の用途には以下のものが含まれますが、これらに限定されません。

  • 音声信号処理: 基本波周波数推定、クロス合成、スペクトル包絡線抽出、時間-スケール変更、タイム ストレッチ、ピッチ シフト。(詳細は異なる合成ウィンドウと解析ウィンドウをもつフェーズ ボコーダーを参照)。

  • きれつ検出: 超音波ラム波の分散曲線を使用した、アルミニウム板のきれつの検出。

  • センサー配列処理: ソナー探査、物理探査、ビームフォーミング。

  • デジタル通信: 周波数ホッピング信号の検出。

使用方法

  • stft は短時間フーリエ変換を計算します。短時間フーリエ変換の逆変換を行うには、関数 istft を使用します。

  • pspectrum または spectrogram はスペクトログラムを計算します。

  • xspectrogram は、2 つの信号のクロス スペクトログラムを計算します。

  • 信号アナライザーのスペクトログラム表示を使用して、信号のスペクトログラムを表示することもできます。

  • pspectrum または信号アナライザーのパーシステンス スペクトル オプションを使用して、他の信号の中に隠れている信号を識別します。

例:パルスと振動

5 kHz で 4 秒間サンプリングされた信号を生成します。信号は、振動する振幅の領域および増加傾向で変動する周波数の領域によって分離された、持続時間が減少していく一連のパルスで構成されています。

fs = 5000;
t = 0:1/fs:4-1/fs;

x = 10*besselj(0,1000*(sin(2*pi*(t+2).^3/60).^5));

信号の短時間フーリエ変換を計算してプロットします。形状係数 β=30 で 200 サンプルのカイザー ウィンドウを使用して、信号にウィンドウを適用します。

stft(x,fs,'Window',kaiser(200,30))

例:減少するチャープをもつ音声信号

減少する 2 つのチャープと広帯域のスプラッター音を含むオーディオ信号を読み込みます。

load splat

オーバーラップの長さを 96 サンプルに設定します。短時間フーリエ変換をプロットします。

stft(y,Fs,'OverlapLength',96)

例:クジラの歌

4 kHz でサンプリングされた太平洋のシロナガスクジラのオーディオ データを含むファイルを読み込みます。ファイルは、コーネル大学の生物音響学研究プログラムが管理する動物発声ライブラリのものです。データの時間スケールは、音の高さを上げ鳴き声を聞き取りやすくするために係数 10 で圧縮されています。

whaleFile = fullfile(matlabroot,'examples','matlab','data','bluewhale.au');
[w,fs] = audioread(whaleFile);

80% のオーバーラップ率でクジラの歌のスペクトログラムを計算します。スペクトログラムの最小しきい値を -50 dB に設定します。

pspectrum(w,fs,'spectrogram','Leakage',0.2,'OverlapPercent',80,'MinThreshold',-50)

例:過渡信号のパーシステンス スペクトル

広帯域信号に組み込まれた狭帯域の干渉信号を読み込みます。

load TransientSig

信号のパーシステンス スペクトルを計算します。両方の信号成分が明瞭に表示されます。

pspectrum(x,fs,'persistence', ...
    'FrequencyLimits',[100 290],'TimeResolution',1)

連続ウェーブレット変換 (スカログラム)

説明

  • ウェーブレット変換は、時間シフトと時間スケーリングを維持する線形時間-周波数表現です。

  • "連続ウェーブレット変換" は、非定常信号の過渡状態の検出や、瞬時周波数が急速に増大する信号に適しています。

  • CWT は可逆変換です。

  • CWT は、可変サイズのウィンドウで時間-周波数平面をタイリングします。ウィンドウの時間の幅は、低周波数の現象に適するように、自動的に広がります。また、高周波数の現象の場合は狭くなります。

考えられる用途

この時間-周波数手法の用途には以下のものが含まれますが、これらに限定されません。

  • 心電図 (ECG): ECG 信号の臨床的に最も有用な情報は、連続した波形の時間間隔と特徴によって定義される振幅にあります。ウェーブレット変換は ECG 信号をスケールに分解し、さまざまな周波数範囲の ECG 信号の解析を容易にします。

  • 脳波図 (EEG): 生の EEG 信号は空間分解能と S/N 比が低く、アーティファクトが発生します。ノイズ信号の連続ウェーブレット分解は、ノイズのランダム分布を変更せずに、大きな絶対値を持ついくつかのウェーブレット係数に固有の信号情報を集中させます。そのため、ウェーブレット係数をしきい値処理することでノイズ除去ができます。

  • 信号復調: 適応ウェーブレット構成手法を使用して、"拡張 2 位相シフト キーイング (EBPSK)" を復調します。

  • 深層学習: CWT を使用して、畳み込みニューラル ネットワークの学習に使用できる時間-周波数表現を作成できます。ウェーブレット解析と深層学習を使用した時系列の分類 (Wavelet Toolbox)は、スカログラムと転移学習を使用して ECG 信号を分類する方法を示しています。

使用方法

  • cwt (Wavelet Toolbox) は連続ウェーブレット変換を計算し、スカログラムを表示します。または、cwtfilterbank (Wavelet Toolbox) を使用して CWT フィルター バンクを作成し、関数 wt (Wavelet Toolbox) を適用します。この手法は、並列化アプリケーションで実行する場合や、ループ内にある複数の関数の変換を計算する場合に使用します。

  • icwt (Wavelet Toolbox) は連続ウェーブレット変換の逆変換を行います。

  • 信号アナライザーには、時系列の CWT を可視化するスカログラム表示があります。

例:ECG 信号

360 Hz でサンプリングされたノイズの多い ECG 波形を読み込みます。

load ecg
Fs = 360;

連続ウェーブレット変換を計算します。

cwt(ecg,Fs)

ECG データは MIT-BIH Arrhythmia Database から取得されたものです [2]。

Wigner-Ville 分布

説明

  • "Wigner-Ville 分布" (WVD) は、信号をその信号の時間と周波数が変換された複素共役バージョンと相関させることで計算される二次エネルギー密度です。

  • Wigner-Ville 分布は、信号が複素数の場合でも常に実数です。

  • 時間周辺密度と周波数周辺密度はそれぞれ、瞬時電力とスペクトル エネルギー密度に対応しています。

  • 瞬時周波数と群遅延は、Wigner 分布の局所的 1 次モーメントを使用して評価できます。

  • WVD の時間分解能は、入力サンプルの数に等しくなります。

  • Wigner 分布は、局所的に負の値を想定できます。

考えられる用途

この時間-周波数手法の用途には以下のものが含まれますが、これらに限定されません。

  • 耳音響放射 (OAE): OAE は蝸牛 (内耳) から発せられる狭帯域振動信号で、これがあることは正常な聴力を示します。

  • 量子力学: 古典的な統計力学に対する量子補正、電子伝達のモデル化、多体量子システムの静的特性と動的特性の計算。

使用方法

例:耳音響放射

20 kHz でサンプリングされた耳音響放射データを含むデータ ファイルを読み込みます。この放射は 25 ミリ秒から開始して 175 ミリ秒で終了するスティミュラスによって生成されました。

load dpoae
Fs = 20e3;

耳音響データの平滑化疑似 Wigner-Ville 分布を計算します。簡易プロットは、放射の周波数をほぼ期待値の 1.2 kHz で分離します。

wvd(dpoaets,Fs,'smoothedPseudo',kaiser(511,10),kaiser(511,10),'NumFrequencyPoints',4000,'NumTimePoints',3990)

耳音響放射の詳細については、CWT による時間-周波数解析 (Wavelet Toolbox)の「解析 CWT による正確な周波数の決定」を参照してください。

再割り当てとシンクロスクイージング

説明

  • "再割り当て" では、スペクトル推定の局所化が鮮明になり、読み取りと解釈の容易なスペクトログラムが作成されます。この手法では、各スペクトル推定はビンの幾何学的中心ではなく、そのビンのエネルギー中心に移動されます。これにより、チャープとインパルスの厳密な局所化が行われます。

  • "フーリエ シンクロスクイーズド変換" は、短時間フーリエ変換から開始し、時間-周波数平面で瞬時周波数の曲線の周囲に集中するようにその値を "押し込み" ます。

  • "ウェーブレット シンクロスクイーズド変換" は、信号エネルギーを周波数で再割り当てします。

  • フーリエ シンクロスクイーズド変換とウェーブレット シンクロスクイーズド変換は両方とも可逆変換です。

  • 再割り当て手法とシンクロスクイージング手法は、時間-周波数 "リッジ" の追跡と抽出に特に適しています。

考えられる用途

この時間-周波数手法の用途には以下のものが含まれますが、これらに限定されません。

  • 音声信号処理: シンクロスクイージング変換 (SST) は、元はオーディオ信号解析のコンテキストで導入されました。

  • 地震データ: 石油トラップやガス トラップを見つけるための地震データの解析。シンクロスクイージングにより、地震データ内の通常は不鮮明な深層の弱い信号も検出できます。

  • 電力系統の振動: 蒸気タービンと発電機は、さまざまなタービン段と発電機の間に機械的準同期振動 (SSO) モードをもつことがあります。SSO の周波数は一般に 5 Hz ~ 45 Hz であり、多くの場合、モード周波数は互いに近くなります。WSST のノイズ低減機能と時間-周波数分解能により、時間-周波数の表示の可読性が向上します。

  • 深層学習: シンクロスクイーズド変換を使用して、時間-周波数の特徴を抽出し、時系列データを分類するネットワークに送ることができます。深層学習を使用した波形セグメンテーションは、ECG 信号を分類する LSTM ネットワークに fsst 出力を送る方法を示しています。

使用方法

  • spectrogram'reassigned' オプションを使用して、pspectrum の引数 'Reassigned'true に設定するか、信号アナライザーのスペクトログラム表示内の [再割り当て] ボックスをオンにして、再割り当てされたスペクトログラムを計算します。

  • fsst はフーリエ シンクロスクイーズド変換を計算します。フーリエ シンクロスクイーズド変換の逆変換を行うには、関数 ifsst を使用します。(ifsst を使用した音声信号の再構成については、音声信号のフーリエ シンクロスクイーズド変換を参照してください)。

  • wsst (Wavelet Toolbox) はウェーブレット シンクロスクイーズド変換を計算します。ウェーブレット シンクロスクイーズド変換の逆変換を行うには、関数 iwsst (Wavelet Toolbox) を使用します。(iwsst (Wavelet Toolbox) を使用した二次チャープの再構成については、Inverse Synchrosqueezed Transform of Chirp (Wavelet Toolbox)を参照してください)。

例:反響定位パルス

オオクビワコウモリ (Eptesicus Fuscus) の発する反響定位パルスを読み込みます。サンプリング間隔は 7 マイクロ秒です。

load batsignal
Fs = 1/DT;

信号の再割り当てされたスペクトログラムを計算します。

subplot(2,1,1)
pspectrum(batsignal,Fs,'spectrogram','TimeResolution',280e-6, ...
    'OverlapPercent',85,'MinThreshold',-45,'Leakage',0.9)
subplot(2,1,2)
pspectrum(batsignal,Fs,'spectrogram','TimeResolution',280e-6, ...
    'OverlapPercent',85,'MinThreshold',-45,'Leakage',0.9,'Reassign',true)

この例では、イリノイ大学 Beckman Center の Curtis Condon 氏、Ken White 氏、Al Feng 氏にコウモリのデータの提供および使用許可をいただきました。ご協力に謝意を申し上げます [3]。

例:音声信号

女性と男性が発声する「strong」という単語が含まれているファイルを読み込みます。この信号は 8 kHz でサンプリングされています。これらを 1 つの信号に連結します。

load Strong
x = [her' him'];

信号のフーリエ シンクロスクイーズド変換を計算します。形状係数 β=20 のカイザー ウィンドウを使用して、信号にウィンドウを適用します。

fsst(x,Fs,kaiser(256,20),'yaxis')

例:合成地震データ

100 Hz で 1 秒間サンプリングされた合成地震データを読み込みます。

load SyntheticSeismicData

Bump ウェーブレットとオクターブあたり 30 の音の数を使用して、地震データのウェーブレット シンクロスクイーズド変換を計算します。

wsst(x,Fs,'bump','VoicesPerOctave',30,'ExtendSignal',true)

地震信号は、『Time-Frequency Analysis of Seismic Data Using Synchrosqueezing Transform』(Ping Wang、Jinghuai Gao、Zhiguo Wang 著) で説明されている 2 つの正弦波を使用して生成されます [4]。

例:地震振動

地震の条件下で 3 階建ての試験構造物の 1 階で記録された加速度の測定値を読み込みます。この測定値は 1 kHz でサンプリングされています。

load quakevib
Fs = 1e3;

加速度の測定値のウェーブレット シンクロスクイーズド変換を計算します。循環動作を示す振動データを解析します。シンクロスクイーズド変換を使用すると、約 11 Hz で区切られた 3 つの周波数成分を分離できます。主要な振動周波数は 5.86 Hz で、等間隔の周波数ピークはこれらが調和的に関連していることを示しています。振動の循環動作も見られます。

wsst(gfloor1OL,Fs,'bump','VoicesPerOctave',48)
ylim([0 35])

例:阪神淡路大震災のデータ

1995 年の阪神淡路大震災の発生時に記録された地震計データを読み込みます。このデータのサンプルレートは 1 Hz です。

load kobe
Fs = 1;

地震データのさまざまな周波数成分を分離するウェーブレット シンクロスクイーズド変換を計算します。

wsst(kobe,Fs,'bump','VoicesPerOctave',48)
ylim([0 300])

このデータは、オーストラリアのホバートにあるタスマニア大学の地震計で 1995 年 1 月 16 日 20:56:51 (GMT) から 51 分間にわたって 1 秒間隔で記録された測定値 (垂直加速度、nm/sq.sec) です [5]。

例:電力系統の準同期振動

電力系統の準同期振動データを読み込みます。

load OscillationData

Bump ウェーブレットとオクターブあたり 48 の音の数を使用して、ウェーブレット シンクロスクイーズド変換を計算します。4 つのモード周波数は、15 Hz、20 Hz、25 Hz、および 32 Hz です。15 Hz および 20 Hz でのモードのエネルギーは時間と共に減少するのに対し、25 Hz および 32 Hz でのモードのエネルギーは時間と共に徐々に増加することに注目してください。

wsst(x,Fs,'bump','VoicesPerOctave',48)
ylim([10 50])

この合成準同期振動データは、『Application of Synchrosqueezed Wavelet Transforms for Extraction of the Oscillatory Parameters of Subsynchronous Oscillation in Power Systems』で Zhao らにより定義されている方程式を使用して生成されました [6]。

Q ガボール変換

説明

  • "定 Q 非定常ガボール変換" はさまざまな中心周波数と帯域幅のウィンドウを使用し、帯域幅に対する中心周波数の比率 (Q 係数) は一定のままです。

  • Q ガボール変換では、安定した逆変換の構成が可能になり、信号の完全再構成ができます。

  • 周波数空間では、ウィンドウは対数的に等間隔の中心周波数を中心とします。

考えられる用途

この時間-周波数手法の用途には以下のものが含まれますが、これらに限定されません。

音声信号処理: 音楽のトーンの基本周波数は、幾何学的に配置されています。人間の聴覚系の周波数分解能はほぼ定 Q であるため、この手法は音楽信号処理に適しています。

使用方法

  • cqt (Wavelet Toolbox) は定 Q ガボール変換を計算します。

  • icqt (Wavelet Toolbox) は定 Q ガボール変換の逆変換を行います。

例:ロック音楽

ボーカル、ドラム、およびギターによるロック音楽の断章を含むオーディオ ファイルを読み込みます。この信号のサンプルレートは 44.1 kHz です。

load drums

CQT が対数周波数応答を持つ周波数範囲を、最小許容周波数が 2 kHz になるように設定します。オクターブごとに 20 のビンを使用して、信号の CQT を実行します。

minFreq = fs/length(audio);
maxFreq = 2000;
cqt(audio,'SamplingFrequency',fs,'BinsPerOctave',20,'FrequencyLimits',[minFreq maxFreq])

経験的モード分解とヒルベルト・ファン変換

説明

  • "経験的モード分解" は、信号を "固有モード関数" に分解し、元の信号の完全でほぼ直交する基底を形成します。

  • "ヒルベルト・ファン変換" は、各固有モード関数の瞬時周波数を計算します。

  • これらの 2 つの手法を組み合わせると、非線形信号と非定常信号の解析に役立ちます。

考えられる用途

この時間-周波数手法の用途には以下のものが含まれますが、これらに限定されません。

  • 生体信号処理: 大脳皮質の経頭蓋磁気刺激 (TMS) に対する人間の EEG 応答を解析します。

  • 構造物に関する用途: 梁や板にきれつ、層間剥離、剛性損失として現れる異常を見つけます。

  • システム同定: 近接した間隔のモーダル周波数をもつ構造のモーダル減衰比を分離します。

  • 海洋工学: 水中の電磁環境で人間によって引き起こされる過渡電磁による外乱を特定します。

  • 太陽物理学: 黒点データの周期的な成分を抽出します。

  • 大気乱流: 安定した境界層を観察して、乱流運動と非乱流運動を分離します。

  • 疫学: デング熱などの伝染病の移動速度を評価します。

使用方法

  • emd は経験的モード分解を計算します。

  • hht は経験的モード分解のヒルベルト・ファンスペクトルを計算します。

例:ベアリング振動

振動信号のヒルベルト スペクトルの計算の例で生成された欠陥ベアリングの振動信号を読み込みます。この信号は 10 kHz のレートでサンプリングされています。

load bearingVibration

信号の最初の 5 つの固有モード関数 (IMF) を計算します。最初と 3 番目の経験的モードのヒルベルト スペクトルをプロットします。最初のモードでは、ベアリングの外輪への高周波数の影響による摩耗の増加が明らかになります。3 番目のモードは、ベアリングの欠陥の原因となった、測定プロセス中に発生している共振を示しています。

imf = emd(y,'MaxNumIMF',5,'Display',0);
subplot(2,1,1)
hht(imf(:,1),fs)
subplot(2,1,2)
hht(imf(:,3),fs,'FrequencyLimits',[0 100])

参照

[1] The Pacific blue whale file is obtained from the library of animal vocalizations maintained by the Cornell University Bioacoustics Research Program.

[2] Moody G. B, Mark R. G. The impact of the MIT-BIH Arrhythmia Database. IEEE Eng in Med and Biol 20(3):45-50 (May-June 2001). (PMID: 11446209)

[3] Thanks to Curtis Condon, Ken White, and Al Feng of the Beckman Center at the University of Illinois for the bat echolocation data.

[4] Wang, Ping, Gao, J., and Wang, Z. Time-Frequency Analysis of Seismic Data Using Synchrosqueezing Transform, IEEE Geoscience and Remote Sensing Letters, Vol 12, Issue 11, Dec. 2014.

[5] Seismograph (vertical acceleration, nm/sq.sec) of the Kobe earthquake, recorded at Tasmania University, Hobart, Australia on 16 January 1995 beginning at 20:56:51 (GMTRUE) and continuing for 51 minutes at 1 second intervals.

[6] Zhao et al. Application of Synchrosqueezed Wavelet Transforms for Extraction of the Oscillatory Parameters of Subsynchronous Oscillation in Power Systems MDPI Energies; Published 12 June 2018.

[7] Boashash, Boualem. Time-Frequency Signal Analysis and Processing: A Comprehensive Reference Elsevier, 2016.

参考

アプリ

関数