時間——周波数解析を使用したDTMF信号内の番号識別

ほとんどの時変信号は,それがどういったものでも各セクションの信号が実質的に定常的となる程度に短い区間に分割できます。時間-周波数解析では、信号をそういった短い時間にセグメント化し、複数のスライディング ウィンドウにわたってスペクトルを推定する方法が最も一般的です。関数pspectrumを的谱图オプションと共に使用して,各スライディングウィンドウ全体に対するFFTベースのスペクトル推定を計算し,信号の周波数成分がどのように経時変化するかを可視化できます。

デジタル電話のダイヤルの信号伝達システムについて考えます。このようなシステムで生成される信号は,デュアルトーン多重周波数(DTMF)信号としてしられています。ダイヤル番号ごとに生成される音は,互いに排他な2つのグループから取得された周波数をもつ2つの正弦波 $-$またはトーン $-$の和で構成されます。トーンの各ペアは,低群(697 Hz, 770 Hz, 852 Hzまたは941 Hz)の1つの周波数と高群(1209 Hz, 1336 Hzまたは1477 Hz)の1つの周波数で構成され,固有のシンボルを表します。以下は電話パッドのボタンに割り当てられている周波数です。

DTMF信号を生成して聞きます。

[tone, Fs] = helperDTMFToneGenerator();p = audioplayer(音调、Fs、16);玩(p)

信号を聞くと3桁の番号がダイヤルされたことがわかります。ただし,どの番号かはわかりません。次に,650 ~ 1500 Hz帯域の時間領域と周波数領域で信号を可視化します。関数pspectrumの“漏”パラメーターを1に設定し,箱型ウィンドウを使用して周波数分解能を向上させます。

N =元素个数(音调);t = (0: n - 1) / Fs;次要情节(2,1,1)阴谋(1 e3 * t,音调)包含(“时间(ms)”) ylabel (“振幅”)标题(DTMF信号的次要情节(2,1,2)pspectrum(音调,Fs,“漏”, 1“FrequencyLimits”(650、1500))

信号の時間領域プロットでは,プッシュされた3つのボタンに対応して3つのエネルギーのバーストが存在することが確認できます。バーストの長さを測定するには,RMS包絡線のパルス幅を取得します。

env =信封(音调,80,“rms”）;脉冲宽度(env, Fs)

ans =3×10.1041 0.1042 0.1047

标题(“均方根包络线的脉冲宽度”）

このように,それぞれにおいて約100ミリ秒の長さのパルスが3つ確認できます。ただし,どの番号がダイヤルされたのかはわかりません。周波数領域プロットでは信号に存在する周波数が示されるため,ダイヤルされた番号の判別に役立ちます。

4つの異なる周波数帯域の平均周波数を推定することによって周波数ピークの位置を特定します。

f = [meanfreq(tone,Fs，[700 800])，．..meanfreq(音调,Fs 900 [800]),．..meanfreq(音调,Fs 1000 [900]),．..meanfreq(音调,Fs 1400 [1300])];轮(f)

ans =1×4770 852 941 1336

推定された周波数を電話パッドの図と突き合わせると,ダイヤルされた番号は' 5 ',' 8 'および' 0 'であることがわかります。ただし,周波数領域プロットではダイヤルされた順番を特定できるような時間情報は示されません。可能な組み合わせは' 580 ',' 508 ',' 805 ',' 850 ',' 085 'または058”です。このパズルを解くため,関数pspectrumを使用して,スペクトログラムを計算し,信号の周波数成分が時間と共にどのように変化しているかを観察します。

650 ~ 1500 Hz帯域のスペクトログラムを計算し, $-$10 dBパワーレベルに満たない成分を削除して,主な周波数成分のみを可視化します。トーンの持続時間と時間領域におけるそれらの位置を表示するには,0%のオーバーラップを使用します。

pspectrum(音调,Fs,的谱图，“漏”, 1“OverlapPercent”，0，．..“MinThreshold”, -10,“FrequencyLimits”(650、1500));

スペクトログラムの色は,周波数のパワーレベルを符号化しています。黄色は高いパワーの周波数成分を示し,青色は非常に低いパワーの周波数成分を示します。濃い黄色の水平線は,特定の周波数のトーンが存在することを示します。プロットにはダイヤルされた3つの数字すべてに1336 Hzのトーンが存在することが明確に示されており,そこからすべての数字がキーパッドの2列目にあることがわかります。最も低い周波数770 Hzが最初にダイヤルされたことがプロットからわかります。次は最も高い周波数941 Hzです。中間の周波数852 Hzが最後です。したがって,ダイヤルされた番号は508となります。

最良の信号表現を得るための時間と周波数分解能のトレードオフ

関数pspectrumは信号をセグメントに分割します。セグメントが長いほど周波数分解能が向上し,短いほど時間分解能が向上します。セグメントの長さは,“FrequencyResolution”と“TimeResolution”のパラメーターを使用して制御することができます。周波数分解能または時間分解能の値が指定されていないと,pspectrumは入力信号長に基づいて時間分解能と周波数分解能の間で適切なバランスを検出しようとします。

太平洋のシロナガスクジラの歌のふるえ声の部分を構成する,4 kHzでサンプリングされた次の信号について考えてみます。

负载whaleTrillp = audioplayer (whaleTrill Fs 16);玩(p)

ふるえ声の信号はトーンパルスの列で構成されます。分解能が指定されていない場合と時間分解能が10ミリ秒に設定されている場合に,pspectrumによって取得される時間信号とスペクトログラムを確認します。“漏”パラメーターを1に設定して,箱型ウィンドウを使用します。パルスの時間位置を決定するため,オーバーラップ率を0に設定します。最後に, $-$60 dBの“MinThreshold”を使用して,スペクトログラム表示からバックグラウンドノイズを削除します。

t =(0:长度(whaleTrill) 1) / Fs;图ax1 = subplot(3,1,1);plot(t,whaleTrill) ax2 = subplot(3,1,2);pspectrum (whaleTrill Fs,的谱图，“OverlapPercent”，0，．..“漏”, 1“MinThreshold”, -60) colorbar (ax2,“关闭”) ax3 = subplot(3,1,3);pspectrum (whaleTrill Fs,的谱图，“OverlapPercent”，0，．..“漏”, 1“MinThreshold”, -60,“TimeResolution”10 e - 3) colorbar (ax3“关闭”) linkaxes (ax₁,ax2 ax3],“x”）

pspectrumによって選択された47ミリ秒の時間分解能は,スペクトログラムのすべてのふるえ声パルスの位置を決定するには大きい値です。一方,10ミリ秒の時間分解能は時間領域での各ふるえ声パルスの位置を決定するのに十分です。これによって,少数のパルスにズームする場合でもより明確になります。

xlim ([0.3 - 0.55])

ここで,オオクビワコウモリ(“Eptesicus fuscus”)の発する反響定位パルスで構成される信号を読み込みます。信号は7マイクロ秒のサンプリング間隔で測定されます。信号のスペクトログラムを解析します。

负载batsignalFs = 1 / DT;图pspectrum (batsignal Fs,的谱图）

既定のパラメーター値を使用するスペクトログラムでは,4つの大まかな時間——周波数リッジが表示されます。周波数分解能値を3千赫に減らし,各リッジの周波数変動での詳細を取得します。

pspectrum (batsignal Fs,的谱图，“FrequencyResolution”, 3 e3)

周波数リッジの周波数での位置決めが改善していることを確認します。ただし,周波数と時間の分解能が反比例するため,スペクトログラムの時間分解能は大幅に小さくなります。オーバーラップを99%に設定して,時間枠を滑らかにします。 $-$60 dBの“MinThreshold”を使用して,不要なバックグラウンド成分を削除します。

pspectrum (batsignal Fs,的谱图，“FrequencyResolution”3 e3,．..“OverlapPercent”, 99,“MinTHreshold”, -60)

新しい設定では,反響定位信号の4つの周波数リッジを明確に示すスペクトログラムが得られます。

時間——周波数の再割り当て

4つの周波数リッジを特定できても,各リッジはまだ複数の隣接する周波数ビンに広がっていることがわかります。これは,時間と周波数の両方で使用されたウィンドウ処理法の漏れが原因です。

関数pspectrumは,時間と周波数の両方において各スペクトル推定のエネルギーの中心を推定することができます。各推定値のエネルギーを,この新しい時間および周波数の中心に最も近いビンに再度割り当てると,ウィンドウの漏れの一部を補正できます。“再分配”パラメーターを使用して,これを行うことができます。このパラメーターを真正的に設定すると,再代入された信号のスペクトログラムを計算します。

pspectrum (batsignal Fs,的谱图，“FrequencyResolution”3 e3,．..“OverlapPercent”, 99,“MinTHreshold”, -60,“再分配”,真正的)

周波数リッジは,さらにシャープになり,時間領域での位置決めが改善されました。関数fsstを使用して,信号エネルギーの位置を求めることもできます。これについては次の節で説明します。

時間——周波数リッジの再構成

周波数が時間とともに減少するチャープ信号と最後のスプラット音で構成されている,次の音声記録について考えます。

负载长条木板p = audioplayer (y, Fs, 16);玩(p) pspectrum (y, Fs,的谱图）

時間——周波数平面でリッジを抽出することで”スプラット”音の一部を再構成してみましょう。fsstを使用して,ノイズが含まれるスプラット信号のスペクトルをシャープにし,tfridgeを使用してチャープ音のリッジを識別し,ifsstを使用してチャープを再構成します。この処理で再構成後の信号のノイズが除去されます。

ガウスノイズを”スプラット”音のチャープ部分に追加します。追加されたノイズにより,安いマイクによって収録したオーディオ録音をシミュレートします。時間——周波数のスペクトル成分を確認します。

rng (“默认”) t =(0:长度(y)-1)/Fs;yNoise = y + 0.1*randn(size(y));yChirp = yNoise (t < 0.35);pspectrum (yChirp Fs,的谱图，“MinThreshold”, -70)

フーリエシンクロスクイーズド変換fsstを使用してスペクトルをシャープにします。fsstは,固定時間の周波数のエネルギーを再代入することで,時間——周波数平面でエネルギーの位置を求めます。ノイズを含むチャープのシンクロスクイーズド変換を計算し,プロットします。

fsst (yChirp Fs,“桠溪”）

チャープは時間——周波数平面で位置決めされたリッジとして表示されます。tfridgeを使用してリッジを識別します。変換と共にリッジをプロットします。

海温,[f] = fsst (yChirp Fs);[冰箱，iridge] = tfridge(sst,f,10);helperPlotRidge (yChirp、Fs、冰箱);

次に,リッジインデックスベクトルiridgeを使用してチャープ信号を再構成します。リッジの両側にビンを1つずつ含めます。再構成後の信号のスペクトログラムをプロットします。

yrec = ifsst (sst,皇帝(256年,10),iridge,“NumFrequencyBins”1);pspectrum (yrec Fs,的谱图，“MinThreshold”, -70)

リッジの再構成により,信号からノイズが除去されています。ノイズを含む信号とノイズ除去後の信号を連続して再生し,違いを聞き比べます。

p = audioplayer ([yChirp; 0(大小(yChirp)); yrec], Fs, 16);玩(p);

パワーの測定

複素線形周波数変調(lem)パルスについて考えます。これは一般的なレーダー波形です.1.27マイクロ秒の時間分解能と90%のオーバーラップを使用して信号のスペクトログラムを計算します。

Fs = 1 e8;bw = 60 e6;t = 0:1 / Fs: 10 e-6;IComp = chirp(t，-bw/2,t(end)， bw/2，“线性”, 90) + 0.15 * randn(大小(t));QComp = chirp(t，-bw/2,t(end)， bw/2，“线性”, 0) + 0.15 * randn(大小(t));IQData = IComp + 1i*QComp;segmentLength = 128;pspectrum (IQData Fs,的谱图，“TimeResolution”1.27 e-6“OverlapPercent”, 90)

スペクトログラムの計算に使用されるパラメーターは,线性调频信号の時間——周波数表現を明確にします。pspectrumはパワースペクトログラムを計算します。これは,カラー値がdB単位の真のパワーレベルに対応することを意味します。カラーバーは,信号のパワーレベルが約 $-$4 dBであることを示します。

対数周波数スケールによる可視化

特定の用途では,信号のスペクトログラムが対数周波数スケール上に可視化されている方が望ましい場合もあります。これを行うには,y軸のYScaleプロパティを変更します。たとえば1 kHzでサンプリングされた対数チャープについて考えます。チャープの周波数は10秒間で10 Hzから400 Hzに増加します。

Fs = 1 e3;t = 0:1 / Fs: 10;fo = 10;f1 = 400;y =唧唧声(f1 t, fo, 10日,“对数”）;pspectrum (y, Fs,的谱图，“FrequencyResolution”, 1．..“OverlapPercent”, 90,“漏”, 0.85,“FrequencyLimits”, 1 f / 2)

このチャープのスペクトログラムは,周波数スケールが対数のときに直線になります。

甘氨胆酸ax =;斧子。YScale =“日志”；

3次元ウォーターフォールによる可視化

视图コマンドを使用すると,信号のスペクトログラムを3次元ウォーターフォールプロットとして可視化できます。関数colormapによって表示色を変えることもできます。

Fs = 10 e3;t = 0:1 / Fs: 2;x1 = vco(锯齿(2*pi*t,0.5)，[0.1 0.4]*Fs,Fs);pspectrum (x1, Fs,的谱图，“漏”, 0.8)

colormap视图(-45、65)骨

パーシステンススペクトルを使用した干渉の検出

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

広帯域信号に組み込まれた狭帯域の干渉信号を考えてみます。1千赫で500秒間サンプリングされたチャープを生成します。チャープの周波数は,測定中に180 Hzから220 Hzに増加します。

fs = 1000;t = (0:1 / fs: 500)”;x =唧唧喳喳(220 t, 180 t(结束),)+ 0.15 * randn(大小(t));

0.05信号には振幅がの210 Hz干渉が含まれます。これは,信号の合計持続時間の1/6のみに存在します。

idx =地板(长度(x) / 6);(1: idx) = x (1: idx) + 0.05 * cos(2 *π* t (1: idx) * 210);

区間100 ~ 290 Hzの信号のパワースペクトルを計算します。弱い正弦波はチャープによって不明瞭になります。

pspectrum (x, fs,“FrequencyLimits”290年[100])

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

图colormapparulapspectrum (x, fs,“坚持不懈”，“FrequencyLimits”(100 290),“TimeResolution”, 1)

まとめ

この例では,関数pspectrumを使用して時間——周波数解析を行う方法とスペクトログラムデータとパワーレベルを解釈する方法を説明しました。また,信号の理解を高めるために時間と周波数の分解能を変更する方法とfsst、ifsst,およびtfridgeを使用して,スペクトルをシャープにし,時間——周波数リッジを抽出する方法について説明しました。対数周波数スケールおよび3次元可視化を使用するためのスペクトログラムプロットの設定方法を説明しました。最後に,パーシステンススペクトルを計算して,干渉信号を検出する方法を説明しました。

付録

この例では次の補助関数が使用されています。