TNJ-082: 窓関数という傍流の裏方こそ深い趣きがある(後編)

2021年12月06日
myAnalogに追加

myAnalog のリソース セクション、既存のプロジェクト、または新しいプロジェクトに記事を追加します。

新規プロジェクトを作成

はじめに

窓関数を用いる必要性は、[1]でも示したとおり、FFT におけるスペクトル・リークを低減させるためです。FFT では必須の裏方といえるでしょう。前回(TNJ-081)と今回は、窓関数を用いた FFT における補正係数の算出方法について検討しています。

今回はこれを離散信号にしたときにどうなるかを考えます。しかしこれが離散信号で考えるとつじつまが合わない状態に陥り、かなり難儀した結果としてこの技術ノートを書き上げることができました。

これらの補正係数は、ネットを探しても数値自体や、補正係数算出手順の詳しい説明をあまりみることはありません。しかしこのような傍流のネタ、トピックこそ、『趣き』があるもんだなあと、いたく(難儀もしたので「痛く」)感心するものでありました。

 

『趣き』を『味わう』のが「趣味」なのではないか

前回、今回の技術ノートでは『趣き』というものをキーワードとして、ストーリーを構成しています。漢字の「趣」について想い巡らせると、「趣味」という単語が思い浮かびます。なるほど「趣味」とは、「『趣き』を『味わう』もの」ではないか!と独り納得したものでした(笑)。

「さて、そういう自分の趣味は何かな」と思うと、最近忙しいこともあり、意外と大したものがないなと思うところです。「WEB ラボの執筆が趣味か(汗)」とか…。

「あらためて、そういう自分の趣味は何かな」と思うと、気がつくことがありました。「オレの趣味は電子機器の修理ではないか?」。このところはヒマもないので、「やりたい趣味」というところかもしれません。これまでの測定器あたりの実績(?)で思いつく主なものを挙げてみると、

・オシロ TDS784D の RUN/STOP ボタン接触不良(ボタン・ラバー交換)

・SSG 8460B の周波数表示異常(古い…^o^;。ロータリ・スイッチの接触ワイパ交換)

・ロジアナ 16500C の HDD を SSD 化(HDD メカ的故障の心配がなくなった。といっても簡単ではなかった)

・オーディオ・アナライザ VP-7723A/B(A はヤフオク購入時には完全に不動で、二次電池液漏れでパターンの相当数が切れていた!全部繋ぎなおし、ちゃんと動作した。その A はすでに売却)

・89410A CRT 交換と A35 モジュールのリレー修理(TNJ080 で紹介したもの)

そしてヒマができたらやってみたいのが、図 1 の BCLラジオ [2]、クーガー2200 の再生です。たしか中学のころ購入したものです

図 1. 大掃除で出てきたクーガー2200(修理したいなとワクワク しているが、いつになるやら…)
図 1. 大掃除で出てきたクーガー2200(修理したいなとワクワクしているが、いつになるやら…)

が、最近のコロナ禍での大掃除で出てきました(笑)。電源は入り、音は出ますが、放送を受信できません。「そのうちに修理を」とワクワクしながら(?)、慌ただしい日常を過ごしています。

しかし、「そのどこが『趣き』を『味わっている』のだろうか」と自分でも思いますが(笑)。

 

FFT で抑えておく基本ポイント(再掲)

前回と今回は、離散信号における窓関数の補正係数の『趣き』を『味わって』みましょうというお話です。

前々回の TNJ-080 で FFT について、その一番基本的な抑えておくべきポイントをご紹介しました。今回も窓関数を「離散信号」として仔細に調べていくため、再掲とはなりますが、TNJ-080 の説明を再度ここで行います。

図 2 は 16 サンプルで 2 周期となるピーク電圧 2V の波形を、16 サンプル長で AD 変換したものです。サンプリング・レートが 16sps(sample per sec)、サンプル全長時間は(1 サンプルの中央でサンプリング動作が行われているとして)1 sec、波形自体の周期は 0.5sec になり周波数は 2Hz です。

これを FFT したスペクトルが図 3 です。時間波形 16 サンプル (16 ポイント)を FFT すると、16の周波数ポイント(これを bin といいます)が得られます。スペクトル・リーク(前回示した)もなく、綺麗に 2Hz が検出できています。

なお今回は「1Hz 帯域幅」を考えていきますので、それに合うように、図 2、図 3 は横軸をそれぞれ全長 1sec まで、15Hz まで (サンプリング・レート 16 sps)に修正してあります。

図 3 の最大周波数は 16Hz-1Hz = 15Hz です。これから、

● FFT した離散周波数の最大周波数 bin はサンプリング・レート-周波数ステップ(1 bin の帯域幅)になる

● FFT した離散周波数の周波数ステップはサンプリング周期×サンプル・ポイント数に相当する時間の逆数になる

つまり図 3 の FFT したスペクトルの周波数ステップ(1 bin の帯 域幅)は 1Hz になっています。

ピーク 2V の正弦波を 16 ポイントで FFT していますが、2Hz に 16、14Hz(= 16Hz - 2Hz。折り返し成分)に 16 と、2 本のスペクトルが得られています。この大きさをポイント数 16 で割ると、 16/16 = 1V となります。本来の周波数成分 1V と、折り返しの周 波数成分 1V の合計が、信号のピーク値 2V を形成します。

 

窓関数ノイズ等価帯域幅補正係数を離散信号で考える(まずは理論的なしくみ)

今回は窓関数のノイズ等価帯域幅補正係数を離散信号で検討していきます。図 4 に窓関数補正係数と窓関数ノイズ等価帯域幅を窓関数ごとにまとめたものを示します([3]から転記。TNJ080, TNJ-081 でも示したもの)。

図 2. 2 周期の波形を 16 サンプル AD 変換した時間軸サンプル (TNJ-080 の図からは全長を 1sec に変えている)
図 2. 2 周期の波形を 16 サンプル AD 変換した時間軸サンプル(TNJ-080 の図からは全長を 1sec に変えている)
図 3. 図 2 を FFT した周波数スペクトル (最大周波数は 15Hz、周波数 bin ステップは 1Hz。 TNJ-080 の図からは全長を 15Hz に変えている)
図 3. 図 2 を FFT した周波数スペクトル(最大周波数は 15Hz、周波数 bin ステップは 1Hz。 TNJ-080 の図からは全長を 15Hz に変えている)
図 4. いろいろな窓関数の補正係数([3]の Table 3 より転記)
図 4. いろいろな窓関数の補正係数([3]の Table 3 より転記)

 

ノイズ等価帯域幅補正係数の計算の考え方

ノイズ等価帯域幅補正係数𝛾(本技術ノートでは記号𝛾を用います)の考え方[4, 5]を改めて示します。図 5(TNJ-081 の図 15 再掲)のような、周波数によらずレベルが一定な(一様に分布する)ホワイト・ノイズを入力信号とし、

① これに窓関数𝑤(𝑛)を掛けたもの、つまり窓関数によりフィルタリングしたものを FFTし、そこから得られた𝑚番目の 1 bin の電力Ω𝑃(𝑚)(この 1 bin の帯域幅を𝑓𝐵INとします)を

② 窓関数フィルタの周波数特性の中央(bin の中央)における通過電力密度を基準値とした1Hzの矩形フィルタ内に存在する電力量Ω𝑃SD

に変換する比が「ノイズ等価帯域幅𝑓BIN_𝐴DJ」です。そして𝛾はこの係数です。つまり[前回の TNJ-081 の式(4)再掲]

数式1

さらに𝑓BIN = 1Hzとすれば電力Ω𝑃は 1Hz bin の電力Ω𝑃@1となり、変形させると、

数式2

となり、ノイズ等価帯域幅補正係数𝛾(𝑓BIN= 1Hz時のノイズ等価帯域幅)自体を得ることもできます(本技術ノートのポイントです)。

16 ポイントで FFT するとします。16 ポイント FFT 用の Hann 窓の窓関数𝑤(𝑛)を、離散値で表してみたものを図 6 に示します。以降では、この窓関数の定義を『正しい定義となるように』変えて考えていきますので[このままの考え方では間違っているので(汗)]、その様子(違い)が分かりやすくなるように、縦軸は対数軸でプロットしてみました。

サンプリング周波数𝑓𝑆 = 16Hz、サンプル数 16 ポイントなので、サンプル/窓関数全長は 1sec、図 3 のとおり周波数ステップ(1 bin の帯域幅)は 1Hz になります。この技術ノートでは 1Hz 帯域幅について考えていますが、ここでは 1 bin = 1Hz としていますので、とくに「bin の帯域幅で割る」という操作をする必要はありません。

しかし以降に示していくように、このままの窓関数の考え方では「なんだかおかしいぞ!」という状態になっており、「おお窓関数!かなり『趣き』あり」というストーリー展開になっています…。

その以降では、途中で「間違った計算/考え方」が入り込んでいますのでご注意ください。なお「間違った計算/考え方」のあたりでは、その部分を赤でハイライトしますから、おかしい部分は分かるようになっていますので、ご心配なく。

図 5. 周波数によらずレベルが一定な(一様に分布する)ホワイ ト・ノイズが入力信号となり、これを窓関数もしくは矩形フィ ルタでフィルタリングする(TNJ-081 の図 15 再掲)
図 5. 周波数によらずレベルが一定な(一様に分布する)ホワイト・ノイズが入力信号となり、これを窓関数もしくは矩形フィルタでフィルタリングする(TNJ-081 の図 15 再掲)
図 6. Hann 窓を 16 ポイントの離散値で表してみた(サンプリン グ周波数 16Hz。縦軸は対数表示)
図 6. Hann 窓を 16 ポイントの離散値で表してみた(サンプリング周波数 16Hz。縦軸は対数表示)

しかし信号処理やスペクトル解析ご専門の方は、「おぬし、そんなことも知らないのか(鉄拳!)」とおっしゃるネタかもしれません…(汗)。

 

離散値の窓関数を掛けたときの電力を計算する

図 5 のような、周波数によらず一様レベルの入力信号(ホワイト・ノイズ)のサンプル列に窓関数を掛けたとき、FFT 出力でどのような電力が得られるかを考えてみます。

ここで「電力」と表現していますが、本来、サンプリング・データ列の全長は任意長です。その全体は「エネルギ」であり、 1sec が基準となる「電力」ではありません。しかしここではサンプリング・データ列の全長を 1sec として、話しを分かりやすくしていますので、「1sec のエネルギ」=「電力」となります。

また入力されるホワイト・ノイズは電圧密度として 1V/√Hz だと仮定します。そうすれば窓関数の周波数特性自体が、窓関数でフィルタされ出力に現れるノイズ量(ノイズ密度)に相当するものとして計算できます。

窓関数を時間軸で𝑤(𝑛)(ただし𝑛 = 0, 1, ⋯ , 15)とし、𝑊(𝑓) を𝑤(𝑛)のフーリエ変換、つまり周波数軸表示とします。𝑊(𝑓) は𝑤(𝑛)の窓関数により得られるフィルタの周波数軸特性を表していることになります。これを図 7 に示します。

「離散信号の𝑤(𝑛)から連続スペクトルの𝑊(𝑓)???が得られる???」と思われるかもしれませんが、入力信号の周波数は bin のステップ間隔ぴったりではない場合もあると考えれば、「連続スペクトル」というのも腹落ち感があるのではないでしょうか。なおこのお話しは、[6]を再度ご覧いただくとご理解できると思いますし、この図 7も[6]の考え方を使って LTspiceで計算してみたものです。

入力信号は周波数によらずレベルが一定(一様分布)、電圧密度 1V/√Hz と仮定し、またサンプリング・データ列の全長を 1sec として、以降の話しを分かりやすくしています。

入力信号(上記の仮定)に窓関数を掛け合わせた出力(窓関数でフィルタした出力)で得られる電力𝑃𝑊_𝐴llは、𝑊(𝑓)を自乗し、周波数領域で積分したものになります。またプラス・マイナスの周波数がありますから、図 7 を折り返したプラス・マイナスの周波数範囲全長として

数式3

で−𝑓𝑆/2~+𝑓𝑆/2の積分で計算します。ここで𝑓𝑆はサンプリング周波数です(この例では 16Hz)。入力信号が上記仮定により、周波数によらずレベルが一定(一様分布)でしたから、結局図 7 のスペクトル全体のエネルギが 1 bin の帯域落ち込んでくることになります。つまり式(3)の答え𝑃𝑊_𝐴llが 1Hz 帯域幅の bin の電力Ω𝑃@1となるわけです。

なお単に自乗し積分しているだけで「電力」としているのは、1秒あたりとして規定していること、また電力の式で考えれば 𝑃 = 𝑉2/𝑅で𝑅 = 1としているからです。

この式を「計算できれば」、窓関数を通して 1Hz 帯域幅の 1 bin に現れる電力Ω𝑃@1を計算できます。しかしこの計算は簡単そうではありません…。

図 7. 窓関数により得られるフィルタの周波数軸特性。サンプリ ング周波数 16Hz(LTspice を用いて計算。図 6 の Hann 窓係数を 使用。以降の説明のように、本当は図 11 による Hann 窓係数を 使用したほうが正しい)
図 7. 窓関数により得られるフィルタの周波数軸特性。サンプリング周波数 16Hz(LTspice を用いて計算。図 6 の Hann 窓係数を使用。以降の説明のように、本当は図 11 による Hann 窓係数を使用したほうが正しい)

 

パーセバルの定理の考えを使って時間軸で計算すればよい

この簡単そうではない計算は、違う道筋からその答えを得ることができます。「時間軸で等価的に求めてみる」というアプローチです。

ここで離散信号についてのパーセバルの定理[7]を取り上げます。

数式4

𝑋は𝑥のフーリエ変換、𝜔は角周波数です。この式ではサンプリング周波数𝑓𝑆は 1Hz として表現しています(先の「エネルギ」と「電力」の話しと同様、一般的にもサンプリング全長は 1sec ではありませんので、この両辺の単位は電力ではなく、エネルギになります)。またサンプル・ポイント数を 16 で表記しています。右辺は変数が角周波数𝜔 [rad/sec]なので、この右辺は連続関数です[式(3)と同じ考えです]。ここで𝜔 = 2𝜋𝑓(𝑓は周波 数 [Hz])として右辺を置換積分します。まず

数式5

を得ます。これより上記の式(4)の右辺は

数式6

ここで𝑥(𝑛)を窓関数の時間軸表現𝑤(𝑛)とすると

数式7

𝑊は𝑤のフーリエ変換です。

この式はサンプリング周波数𝑓𝑆 = 1Hz と仮定しています。積分範囲が±1/2 になっていますが、式(3)においては、

● サンプリング・データ列の全長を 1sec

● サンプリング周波数𝑓𝑆は 16Hz

であったため、𝑓𝑆が 1Hz 以外であれば、その全体±𝑓𝑆/2で積分します。

数式8

そうするとこれは、式(3)の右辺と同じ表現になるわけですね。

さらにサンプリング・データ列の全長は 1sec でしたから電力単位となり、この式(8)の右辺は式(3)そのものとして、「窓関数を掛け合わせた(窓関数でフィルタした)FFT 出力で得られる 1Hz 帯域幅の bin の電力Ω𝑃@1」を求めていることになるわけです。

さらにこの式(8)の左辺は時間サンプルを表しています。つまり式(3)の𝑃𝑊_𝐴ll = Ω𝑃@1を求めたいのであれば、周波数軸での積分ではなく、式(8)の左辺、窓関数の時間サンプル𝑤(𝑛)において

数式9

を計算すればよいという結論に辿りつきます。時間サンプル全長 1sec の𝑤(𝑛)を自乗して足し合わせることで、窓関数でフィルタされ出力に現れる(電圧密度 1V/√Hz と仮定した)ノイズの電力を計算できるということです。

 

やりたいことは「ノイズ等価帯域幅補正係数𝜸を求める」

ここでやりたいことは、周波数によらずレベルが一定(一様分布)で、電圧密度 1V/√Hz の入力信号の、1Hz 帯域幅の bin の電力Ω𝑃@1自体を求めることではありません。ノイズ等価帯域幅補正係数[式(2)再掲]

数式2-2

つまり「比」が求められればいいわけで、入力信号𝑥(𝑛)の時間形状には影響されないことにも気がつきます。

 

実際に MATLAB/Octave で窓関数ノイズ等価帯域幅補正係数を計算してみる(でもこの計算は何 だかおかしいぞ!)

MATLAB/Octave を使って、窓関数ノイズ等価帯域幅補正係数を計算してみましょう(以下は視覚的な理由で、一部の表示を修正してあります)。しかしこの考え方は間違っているようなのです…。多くの記事や書籍でこのようなかたちで窓関数の設定を説明や計算をしているにも関わらず…。

 

Hann 窓の離散値を得ることと基本的な設定条件

ではやってみましょう。まず 16 サンプルの Hann 窓を得ます。

数式9-2

なお w は行ベクトルです(視覚的な理由で表示を修正)。実はこの部分の考え方(Hann 窓の得かた)が間違っているのです! この Hann窓は、サンプリング周期を 1/16 sec(サンプリング周波数𝑓𝑆 = 16Hz)、サンプル数を 16 サンプルとすることで全長が 1sec となります。この条件で1 bin が 1Hz 帯域幅となります。

上記に説明したように、窓関数ノイズ等価帯域幅補正係数𝛾は式(2)のように、1Hz 帯域幅の bin の電力Ω𝑃@1と 1Hz の矩形フィルタ内に存在する電力量Ω𝑃SDとの「比」ですから、入力信号の形状には影響されず、計算にも𝑥(𝑛)は含まれていません。

なお、この計算において周波数軸で一様分布の信号を時間軸で作りたいなら、デジタル信号処理としては「インパルス信号」になりますから、この計算手順が応用できなくなります(笑)。

 

Hann 窓により形成されたフィルタ帯域内のノイズ電力を計算する

窓関数𝑤(𝑛)というフィルタを通したとき 1 bin帯域に現れる電力 Ω𝑃を計算するには、示してきたとおり、式(9)のように電力相当 として Hann 窓𝑤(𝑛)を自乗して足し算します。

ここでは 1 サンプルにおける占有時間は 1/16 sec なので、個々のサンプルに相当する電力は、その値を 1/16 にする必要があります。実効値の計算方法と同じです。この考え方から 1 秒間の電力値に相当する大きさ(1Hz 帯域幅の bin の電力)Ω𝑃@1は MATLAB/ Octave では、

数式9-3

としてオール 1 との内積で計算でき、

数式10

が得られます(実はこの値が間違っているのです!)。全長が 1sec ですから、そのまま 1sec のエネルギ = Watt(J/s = Joule per second)と計算できます。

図 8. サンプリング周波数が 16Hz、サンプル数が 16 ポイント(全長 1sec)で、bin の帯域幅が 1Hz と なるときの 1Hz 矩形フィルタの考え方
図 8. サンプリング周波数が 16Hz、サンプル数が 16 ポイント(全長 1sec)で、bin の帯域幅が 1Hz と なるときの 1Hz 矩形フィルタの考え方

 

サンプリング周波数 16Hz、全長 1sec の条件で帯域幅 1Hz の 矩形フィルタをどう考えるか

つづいてこのサンプリング条件において、帯域幅 1Hz の矩形フィルタで得られる電力を計算してみます。

ここまで検討してきた条件では、図 8 のように 1 bin の帯域幅が 1Hz でした。「帯域幅 1Hz の矩形フィルタ」というのは、この場合では 1 bin 単体の通過特性を矩形と見立てたものです。そこで得られた電力が、帯域幅 1Hz の矩形フィルタ出力での電力、つまり電力密度になります。

 

各 bin で帯域幅 1Hz の矩形フィルタの大きさを得る

ここで「帯域幅 1Hz の矩形フィルタ」を通過してきた大きさは、窓関数を掛けた操作によるフィルタ中央(各 bin の中央)の通 過量(図 7 で示したもの)と同じと考えます。

この各 bin の中央における通過量を計算するには、0Hz に相当する binでの通過量を考えれば一番簡単です。窓関数に加わる入力信号を 0Hz つまり DCとし、その大きさを 1[𝑥(𝑛) = 1]とすれば、窓関数の形状自体が 0Hz に相当する成分の通過量と考えることができるわけです。そしてその通過量の計算は、窓関数の形状自体(時間軸の)の平均値を得ればよいことになります。

他の(DC 以外の)bin がどうなるかは、この考え方を拡張して考えれば簡単です。想定している入力信号は、「周波数によらず一定で一様に分布するホワイト・ノイズ」としています。離散フーリエ変換結果の各 binの中心に相当する周波数において、大きさ 1 の信号(一様分布のノイズ)の変換結果として得られたものも、DC の場合と同じだと考えることができます。

 

窓関数離散値の平均値を計算する

さて、それでは窓関数の離散値𝑤(𝑛)の平均値を得てみます。本題に立ち返って、もう一度ご説明させていただきます。

𝑥(𝑛) = 1に窓関数𝑤(𝑛)を掛け合わせた時間サンプル、𝑥(𝑛) ∙ 𝑤(𝑛) = 𝑤(𝑛)の FFT 結果としては、周波数ゼロの bin のところの大きさに相当します。窓関数の離散値𝑤(𝑛)の平均値はこの大きさと等しくなります。

ところでこの(以降の)算出の考え方は、ひとつまえの TNJ081 の図 9、図 10、図 11 で示した、窓関数減衰補正係数の計算手順と近い(実際は同じ)ものとなっています。

平均値(伝達関数)を求める計算は MATLAB/Octave では、

数式10-2

として先と同様、オール 1 との内積計算です。サンプル・ポイント数が 16 なので 16 で割っています。図 6 の𝑤𝑤(0), …, 𝑤𝑤(15) を内積計算した答えは

数式11

が得られます。

同じく図 6 の𝑤(0), …, 𝑤(15)を FFT してみると、図 9 のように周波数要素𝑊(0), …, 𝑊(15)となり、この要素 0、つまり DC 成 分𝑊(0)は 7.5 になっています。この大きさabs[𝑊(0)] = 7.5を図 3 に関する解説や、前々回の TNJ-080 の式(2)で 1/65536 としたように、16 で割ることで DC 成分(この bin)の大きさが上記の平均値計算の答え 7.5000/16 と同じになります。なるほど…。

 

帯域幅 1Hz の矩形フィルタ出力の電力値を得る

上記の数値、7.5000/16 は各 bin の中央における通過量です。これを 1Hz 帯域で考えますので、自乗して電力相当にし、(1Hz 帯域なので)×1倍したものが、帯域 1Hzの矩形フィルタで得られる電力になります(再度、図 7、図 8 を参照してください)。これを MATLAB/Octaveでは

数式11-2

とすることで

数式12

が得られます。なお単位は Watt としていますが、ここでは実際は 1Hz 帯域で考えていますので、Watt/Hz が正しいというか、本来の単位です。

 

窓関数ノイズ等価帯域幅補正係数を計算すると「答えが合わないぞ!」

それではこれまで得られた結果から、ノイズ等価帯域幅補正係数𝛾を計算してみましょう。式(10)と式(12)を用いて

図 9. 図 6 を 16 ポイントで離散フーリエ変換すると DC 成分 𝑋(0)の大きさは 7.5 になっていることが分かる
図 9. 図 6 を 16 ポイントで離散フーリエ変換すると DC 成分 𝑋(0)の大きさは 7.5 になっていることが分かる

数式13

あれ?答えが合いません!1.5 になっていません!これはなぜでしょう。またまた『趣き』がでてきましたね!(笑)。

 

正しい係数を得るためにはどう考えるのか

なぜ答えが合わないのでしょうか。ここが窓関数の趣き深いところです。そういう私もこの技術ノートの執筆ネタとして、簡単に答えが得られるだろう(執筆が進むだろう)と思っていましたが、上記の計算をしていくなかで、「あれ?なぜ答えが合わないんだ?おかしいな!」と思い立ったのが、このようなストーリーが形成されている原因なのでした。

「執筆しながら考える『回路設計 WEB ラボ』」それをそのまま地で行くモノカキ状態なのでありました。そのためだいぶ(相当…)執筆に時間がかかっているという事実もあるのでした(笑)。そして「窓関数だなんて、傍流の裏方みたいなものだが、とても『趣き』があるもんなのだな」といたく感心したのでした。

離散値では窓関数のサンプル・ポイント数をどう取り扱うのか

[8]にこっそりこんなことが書いてあります([8]より引用)。

数式13-2

'periodic'のところを見てみると([8]より引用)、

数式13-3

簡単にいえば、「スペクトル解析にはこちらを用いる。DFT を考慮したとき完全周期性が得られる」ということです(正確な日本語翻訳は日本語版の[9]をご覧ください)。

図 2、図 3 で示したように、離散フーリエ変換(実用的には Fast Fourier Transform; FFT ですが)で完全な周期性をもつ 16 サンプルで 1 周期となる波形(コヒーレント・サンプリング)であれば、スペクトル・リークなしの周波数スペクトルを得ることができます([1]でも示しました)。

この条件は図 2 を見ていただいても分かるように、17 サンプル目(0 から数えると 16 ポイント目)が次の 1 周期の波形スタートとなり、図 2 の一番右側の 16 サンプル目(0 から数えると 15 ポイント目)がその波形の最後となります。しかしここは「ゼロ値(つまり波形スタート)」ではありません。これが「完全な周期性」ということの裏返しでもあるわけです。

また[10]の 696 ページの式(10.2)の上にある記述と、その式(10.2)自体として以下の記載があります([10]より引用)

数式13-4

数式14

これが何を示そうとしているかというと、「周波数領域での周期的畳み込みで用いられる信号𝑥、そしてそれをフーリエ変換したもの𝑋と、窓関数𝑤、そしてそれをフーリエ変換したもの𝑊は、

① それぞれ同じ長さであり

② 𝑥, 𝑋が「完全な周期性」を要求しているのであれば、 𝑤, 𝑊も「完全な周期性」を要求している(ように読める)

ということに気がつきました。また[11]でも離散窓関数の最終 要素を𝑛 = 𝑁 − 1にしてあること、symmetry, periodic という記述も見つけました。

 

離散値での窓関数の実際の取り扱い方

結論として正しくは、𝑤(𝑛)の 17 サンプル目(0 から数えると 16 ポイント目)が次の 1 周期の波形スタート、w(17) = 0 となるように設定すべき、ということになるわけです。

つまり

数式14-2

で得られた、図 6 でも示した w(16) = 0(0 から数えると 15 ポ イント目)となる状態は、「完全な周期性」を持っておらず、

数式14-3

を得て、ここから

数式14-4

で最初の 16 要素の窓関数にすべき[w(16) ≠ 0 で w(17) = 0] ということです。もしくは

数式14-5

とします。つまり図 10 のように設定すべきということです。こ れは図 6 とは最後のゼロ位置が異なる状態です。

完全周期性を得られる窓関数の設定で計算しなおしてみる

それではここまでの議論の流れを繰り返して、完全周期性のある 16 サンプルの窓関数を用いて計算しなおしてみましょう。まず 17 サンプルの Hann 窓を得ます。

数式14-6

なお w は行ベクトルです。このうち最初の 16 要素を取り出します。

数式14-7

図 10. 16 ポイントの Hann 窓を完全周期性としてみた  (サンプリング周波数 16Hz。縦軸は対数表示)
図 10. 16 ポイントの Hann 窓を完全周期性としてみた(サンプリング周波数 16Hz。縦軸は対数表示)

もしくは(先に示したように)

数式14-8

とします。さきの式(10)の計算と同じように、𝑤(𝑛)という窓関数フィルタを通したとき、1 bin に現れる電力Ω𝑃@1を計算するため、式(9)のように電力相当として w を自乗して足し算します。

数式14-9

こたえは

数式15

が得られます。式(10)で得られた数字とは異なります。

つづいて帯域幅 1Hz の矩形フィルタで得られる電力を計算するため、𝑤(𝑛)の平均値を得ます。その平均値(伝達関数)計算は、

数式15-1

サンプル・ポイント数が 16[hann(17)のうちの前半 16 要素]なので、16 で割っています。こたえは

数式15-2

が得られます。ここでも式(11)で得られた数字とは異なりますね(笑)。これにより

数式16

なお単位は Watt としていますが、ここでは実際は 1Hz 帯域で考えていますので、Watt/Hz が正しいというか、本来の単位です[先の式(12)の説明のコピペのままでスイマセン…]。

式(15)と式(16)からノイズ等価帯域幅補正係数𝛾を計算してみると、

数式17

なんと!これで 1.5 になりました!

式(13)の 1.6 と、この 1.5 の違い(誤差)が 7%程度になっているのは、サンプル数(FFT ポイント数)を少ない 16 で計算しているからです。サンプル数(FFT ポイント数)が大きくなれば、というか一般的な FFT 計算であれば、無視できる誤差になります。そのため、多くのケースで見過ごされてきたのではないかと思います。

しかし窓関数、「趣き深い」ですね!

 

まだまだ趣き深い FFT。解決できていない疑問

ここまで FFT、もう少し原理的な言い方をすると「離散フーリエ変換」における窓関数というものを考えてきました。しかし私は昔から離散フーリエ変換に関する「疑問」があるのです。

疑問

DC 1V を FFT すると 1 が得られる。AC(各 bin 中央の周波数と仮定して)のピーク 1V の時間信号を FFT すると 1 が得られる。

しかしこれまでの複数の本技術ノートで議論してきたように、電力の計算をする場合には、実効値で考える必要がある。DC 1V の実効値は 1 となるが、AC でピーク値 1V というのは実効値 1/√2V であり、「電力の計算」、「実効値」という点では、FFT 結果はつじつまが合わないのではないか?

「そんなことも知らないの?」とか言われそうですが(汗)…。

結局これは、離散フーリエ変換の本質に由来するものと考えることができそうです。結局離散フーリエ変換は各周波数 binにおける、複素周波数参照/基準信号(参照源/基準源)と観測信号を乗算し積分する、つまり 2 信号同士の内積をとる、もっと簡単にいえば、「複素周波数信号と観測信号の二つの信号が、どれだけ似通っているかを示す指標」なわけですから、DC で 1Vは 1V(100%似通っている)だし、ACでも(くりかえしますが bin 中央の周波数では)ピーク 1V の波形は 1V(100%似通っている)と言えることになります。

つまり「離散フーリエ変換では実効値という概念は蚊帳の外」ということなのですね。

私はだいぶ前に「FFT では自乗することを power と書く。ただし power といっても『べき乗』という意味であり、電力のパワーを示すものではない」というのを何かの書籍で読んだことがあります。当時は「同じことだろう(自乗すれば電力になる)」と思っていいました。

しかし上記のように、今回あらためて熟考してみると、DC では離散フーリエ変換結果を自乗すればそのまま電力(抵抗は 1Ω として)になりますが、AC では変換結果を自乗してもそのままでは電力にならない(実効値の考えではない)というわけですから、「powerといっても『電力のパワー』を示すものではない」という説明も頷(うなず)けることになるわけですね。

 

それでも本技術ノートではその疑問点を回避してあいまいさを排除して説明している

式(3)で連続スペクトル𝑊(𝑓)を積分する計算を示しました。この𝑊(𝑓)は𝑊(𝑛)として離散フーリエ変換した値で置き換えて考えることもできます。

そうするとここでご提示した疑問、つまり「DC では実効値がそのまま得られているが、DC 以外の bin ではピーク値であり、実効値を得るにはそれから 1/√2をしなくてはならない」という点から、「おぬしの説明の展開でも、DC 以外では 1/√2 をする必要があるだろう?」というご指摘があろうかと思います。

しかしそれでもここまでの検討では、この疑問はスルーできていたのです。その式(3)以降の展開をご覧いただくと、周波数軸で計算をしているのではなく、式(4)から式(9)の流れのようにパーセバルの定理を用いて、時間軸で計算していることが分かります。

つまりこの疑問、あいまいさを回避して、説明が組み立てられていたわけなのです。

 

まとめ

2 回に亘って(わたって)窓関数ネタをご紹介してきました。あらためて考えてみると、本当に『趣き深い』、「奥深い」ものであることに気が付きます。極論となりますが「窓関数操作」自体は、時間サンプル全長と同じ全長の窓関数を「無理やり」適用するから無理があるわけで、それが『趣き深さ』を作り上げている理由なのかもしれません(フィルタのインパルス応答が時間サンプル全長より長ければ、より良好なフィルタが実現できるため)。

前回では枕草子を取り上げて、「あとから気が付く『趣き』の深さ」という切り出しで技術ノートを書き始めました。書き始めて、いろいろと検討をしてみると、「こたえが合わないぞ!」と自分の知識が足りないことをあらためて思い知らされたわけでした。

仕上がった技術ノートはすらすらと執筆が進んでいるようにも見えるかもしれませんが、複数の壁にぶつかり、実はかなり難儀して、その『趣き』を堪能(?)した、これら技術ノートでありました。さらには「『趣き』を堪能」などと涼しい顔したふりして書いていますが、他業務で相当忙しいなか、「分からん!時間がない!困ったぞー」と焦りながらの執筆だったことを白状させていただきます(笑)。

 

ラジオと離散フーリエ変換

今回の技術ノートの最初に BCL ラジオ、クーガー2200 というものをご紹介しました。

ラジオは「ある周波数」の信号の「大きさ」(とくに AM の場合は大きさ)を検出します。これは離散フーリエ変換で特定の bin の大きさを考えることと同じなんですね。離散フーリエ変換のイメージが難しいとすれば、「ラジオと同様」と考えてもらえると良いかもしれません。

そして「ラジオとアマチュア無線が私の原点なんですよ」ということをお話しし、AD 変換から始まった長い、フーリエ変換ネタに関係する複数の技術ノートを終わりにさせていただきます。

さて次のネタは何にしましょうか…。

著者について

石井 聡
1963年千葉県生まれ。1985年第1級無線技術士合格。1986年東京農工大学電気工学科卒業、同年電子機器メーカ入社、長く電子回路設計業務に従事。1994年技術士(電気・電子部門)合格。2002年横浜国立大学大学院博士課程後期(電子情報工学専攻・社会人特別選抜)修了。博士(工学)。2009年アナログ・デバイセズ株式会社入社、現在に至る。2018年中小企業診断士登録。
デジタル回路(FPGAやASIC)からアナログ、高周波回路まで多...

最新メディア 21

Subtitle
さらに詳しく
myAnalogに追加

myAnalog のリソース セクション、既存のプロジェクト、または新しいプロジェクトに記事を追加します。

新規プロジェクトを作成