LabOneの窓関数でスペクトルのリークをコントロールする

November 25, 2021 by Mehdi Alem

はじめに

信号やシステムのスペクトル特性は、物理学や工学の幅広い分野で必要とされています。信号を周波数領域で分析することにより、時間領域での分析では得られない、信号の挙動の有益な情報を引き出すことができます。そのため、信号のスペクトル解析を、高い信頼性で効率的に行えることは、すべての測定機器にとって非常に重要です。測定された時間信号を、周波数領域の表示に変換するための手段として、フーリエ変換 (FT) が一般的に用いられています。そのデジタル版は離散フーリエ変換 (DFT) と呼ばれ、高速フーリエ変換 (FFT) アルゴリズムを使用してデジタル測定システムで広く使われてされています。しかし、この手法にはいくつかの制限があり、その多くはスペクトルのリークにより生じます。この記事では、信号のスペクトルを取得する際のリークについて検討し、窓関数を使用してどのようにそれを抑制できるかを説明します。また、測定精度の低下を防ぐために、窓関数が信号のスペクトル密度に与える影響を補正する方法についても説明します。

スペクトルのリーク

フーリエ変換は、信号が無限の時間長であり、そのすべてをフーリエ変換に使うことを仮定しています。しかし、信号全体を考慮することは、不可能であったり、膨大な計算量が必要になります。そこで、信号の一部分のみを選択したスペクトル解析、すなわち図1に示すように、信号を長方形の窓で切り出すように行います。

図1: 窓関数処理によるスペクトル漏れのデモ。上図は1 kHzのサイン信号と長さ40 msの窓を表示したプロット。中央のプロットは、窓関数を適用した後の信号を表示。この結果には、40サイクルの正弦信号が含まれています。下の図は、実際の信号 (オレンジ) と窓付き信号 (青) のスペクトルを比較したもので、窓関数による明らかなリークが見られます。

図1の上図に示すような正弦信号から、スペクトル及びその周波数を求める場合を考えます。実際には、図1中段に示すように、信号に長方形の窓をかけるように、信号を切り出し、そのフーリエ変換を計算します。図1の下段は、理想的な正弦波信号と窓関数を掛けた信号のスペクトルを比較したものです。正弦波のスペクトルには、狭い幅のスペクトル線 (理想的にはインパルス関数) となりますが、窓関数によって信号の長さが制限されているため、信号周波数付近に複数のピークを持つ広がったスペクトルとなっています。これは、ある周波数の信号成分が、元の信号には存在しない別の周波数に漏れて生じるため、スペクトル・リークと呼ばれています。スペクトル・リークは、振幅の精度、周波数分解能を低下させます。スペクトル・リークを回避する方法の1つは、このブログ記事に示したように、信号特性に基づいてサンプリングレートと収録長を適切に選択することです。しかし、スペクトル・リークを回避できる状態にパラメータを設定できることは、簡単でない場合が多いです。そこで一般的には、後述するように、信号に合わせた窓関数を適用します。

窓関数

フーリエ変換を計算するために信号の区間を選択することは、図1に示すように、信号に長方形の窓をかけることと同意です。ただしこのような窓関数は、スペクトルのリークという好ましくない結果を引き起こします。ここで、窓の形状を緩やかにすることで、スペクトルのリークを調整し、信号セグメントのDFTを計算する際にその影響を抑制することが可能です[1]。窓関数には2つのタイプがあり、それぞれ異なる目的をもっています。

対称窓
非対称窓

対称窓関数は、中点を中心に対称な形状をしており、周波数領域の測定におけるスペクトル・リークの対策に広く用いられています。一方、非対称窓は中点に対して対称でなく、リングダウン法などの時間領域の測定技術を採用した場合のスペクトル・リークの緩和に主に用いられていています。図2に、Zurich Instruments測定器およびデータ収録ソフトウェアであるLabOne®で利用できる、対称型および非対称型の窓関数をあげました。

図2：LabOneでスペクトル測定と時間測定にそれぞれ使用した対称・非対称の窓関数。

ここでは、スペクトル測定に注目し、対称的な窓関数だけを対象にします。図2の上段に示したように、信号区間のテーパー付けに利用できる窓関数がいくつかあります。ここで問題となるのは、どの窓関数を選択するのがベストなのかということです。もちろん、この答えはアプリケーションの要件により異なりますが、2つの重要な基準があり、それによって、2つの窓関数のセットが存在します。

強い周波数成分から微弱な周波数成分を分離するための高いダイナミックレンジ
同じような強度の2つの近接した周波数成分を分解するための高い周波数分解能

この2つの基準はトレードオフの関係にあり、1つの条件に優れた窓関数は、通常、もう1つの条件には劣っています。LabOneで利用できる高いダイナミックレンジを持つ窓関数として用いられるのはBlackman Harrisであり、高い周波数分解能の典型例はRectangularです。この2つの基準のバランスの取れたものとして、HannとHammingの窓関数があります。LabOneではHannをスペクトル測定のデフォルトの窓関数として使用しています。表1は、LabOneでスペクトラム(対称)とリングダウン(非対称)の測定に使用される窓関数の数学的定義を示したものです。ここで、\(N\)は、\(0 \leq n \leq N\)となる窓の長さを表します。

表1:LabOneに実装されている窓関数の一覧
ウィンドウ	式	特性
Rectangular(ボックスカー)	\(1\)	対称
Hamming	\(0.54 - 0.46\cos(\frac{2\pi n}{N})\)	対称
Hann	\(0.5 - 0.5\cos(\frac{2\pi n}{N})\)	対称
Blackman-Harris	\(0.35875 - 0.48829\cos(\frac{2\pi n}{N}) + 0.14128\cos(\frac{4\pi n}{N}) - 0.01168\cos(\frac{6\pi n}{N})\)	対称
Exponential	\(\exp(\frac{n}{N})\)	非対称
Cosine	\(\cos(\frac{\pi n}{2N})\)	非対称
Cosine squared	\(\cos^2(\frac{\pi n}{2N})\)	非対称

これらの窓関数は、LabOneのスコープ、データ収集（DAQ）、スペクトルの各ツールの設定項目で選べます。実信号に対する窓関数の影響を調べるために、LabOneのスコープモジュールを使って100 kHzの正弦波のスペクトルを測定してみます。このモジュールでは、時間ドメインに加えて、周波数ドメインの表示を、さまざまな窓関数を使って行うことができます。図3は、LabOneで使える4つの対称的な窓関数を使った場合の測定結果を示しています。

図3：MFLIロックインアンプで測定した100kHzの信号のスペクトルをLabOneのScopeモジュールで取得し、4つの窓関数(Rectangular, Hamming, Hann, and Blackman-Harris)で表示したもの。

図3から、窓関数の信号スペクトルへの影響が明確に読み取れ、窓関数の種類を変えると、スペクトルリークが大きく変わることがわかります。信号周波数の近傍では、Rectangular窓からHamming窓、Hann窓、Blackman-Harris窓に変化させると、それぞれスペクトルの拡がりが大きくなります。一方、中心周波数から離れると、Blackman-Harris、Hann、Hamming、Rectangularと逆の順番でリーク量が増えることがわかります。

スペクトル密度

多くの周波数領域測定技術では、被測定システムに含まれる信号やノイズ源のスペクトル密度を求めることが重要な目的です。振幅スペクトル密度であれ、パワースペクトル密度であれ、測定方法の多くは、時間領域の信号を取得し、そのフーリエ変換を計算し、それを測定器の周波数分解能で割るという手法で行います。窓関数を適用しない場合、あるいは矩形窓を実際に使用する場合は、サンプリング周波数\(F_s\)と取得サンプル数\(N\)から、周波数分解能は単純に\(\Delta f = \frac{F_s}{N}\)で求められます。ただし、窓を適用する場合は、周波数分解能は対応する窓関数の等価雑音帯域幅 (ENBW）により変わります。すなわち、等価周波数分解能は\(c_w\Delta f\)で表されます。ここでENBW補正係\(c_w\)は以下の式から求められます。

\[c_w = \frac{\int{|w(t)|^2 dt}}{|\int{w(t)dt}|^2}\]

ここで表2にLabOneで使える4つの対称窓関数のENBW補正係数を示します。

表2: 対称窓関数の等価ノイズ帯域幅
窓	\(c_w\)
Rectangular	\(1.00000\)
Hamming	\(1.36283\)
Hann	\(1.50000\)
Blackman-Harris	\(2.00435\)

ここで、電圧信号\(x\)のパワースペクトル密度 (PSD) \(\frac{\text{V}^2}{\text{Hz}}\)を次のように定義できます。

\[S_x = \frac{|\text{FFT}(x)|^2}{c_w\Delta f}\]

\(\frac{\text{V}}{\sqrt{\text{Hz}}}\)で与えられる振幅スペクトル密度は、PSDの平方根をとることで簡単に得られます。その結果、周波数差分\(\Delta f\)を取得し、FFT計算前に窓関数\(w\)を乗じた信号\(x\)の振幅スペクトル密度は、次のように表されます。

\[\sqrt{S_x} = \frac{|\text{FFT}(x)|}{\sqrt{c_w\Delta f}}\]

LabOneのグラフィカルユーザーインターフェイスでは、スペクトル密度を表示する際にはENBW補正係数\(c_w\)が考慮されます。しかし、信号スペクトルをモジュールヒストリーから記録する場合や、LabOneのAPIを使用する場合は、信号スペクトルしか記録できないため、そのスペクトル密度は得られません。信号のスペクトル密度を求めるには、記録した信号に対して、後処理の段階で上記の式を適用する必要があります。記録されたファイルには、信号のスペクトル、周波数分解能、選択した窓関数の補正係数が含まれています。この情報を使って、次のMATLAB®のコードで示すように、信号スペクトルから信号のスペクトル密度を得ることができます。同様のコードは、PythonやC言語など他の言語でも書くことができます。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 測定器のシリアル番号
device = 'dev3305';

% ファイル 'spectrum_00000.mat' から測定したスペクトルを読み込む
struct = load('spectrum_00000'); 
data = struct.(device).demods.sample_xiy_fft_abs_avg{1,1};

% 周波数情報を抽出する 
freq_offset = data.header.gridcoloffset;
freq_delta = data.header.gridcoldelta;
freq_num = data.header.gridcols;  

% 周波数範囲を計算する
freq = freq_delta*(0:double(freq_num)-1) + freq_offset;

% 読み込んだデータから信号のスペクトルを抽出する
spectrum = data.value(1,:);

% 窓関数の補正係数を求める  
factor_nep_window = data.header.nenbw;

% スペクトル密度の計算
spectral_density = spectrum / sqrt(freq_delta*factor_nep_window);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

まとめ

窓関数は、信号のフーリエ変換を計算する際にスペクトルのリークを抑えるために広く使用されています。LabOneソフトウェアでは、Rectangular, Hamming, Hann, Blackman-Harrisの4つの対称窓関数を使うことができます。スペクトル測定時のスペクトルリークを抑制するために、用途に応じて使い分けることを推奨します。窓関数によって等価雑音帯域幅が変化するため、フーリエ変換から信号のスペクトル密度を求める場合は注意が必要です。

参考文献

Wikipedia: Window function.
Mathuranathan Viswanathan, Equivalent noise bandwidth (ENBW) of window functions, GaussianWaves, 2020.