タグ: 信号処理

  • FFTのゼロ埋め、周波数の刻みと分解能の違い

    FFTのゼロ埋め、周波数の刻みと分解能の違い

    音や振動の記録をFFTにかけると、どの周波数が強いかを山の形で見ることができます。計算する点を増やしたら、その山が滑らかになり、頂点の位置も細かく読めるようになった。けれど、それだけで近接した二つの振動を見分ける力まで増えたのでしょうか。

    この違いを整理する鍵が、ゼロ埋め(zero-padding)です。記録済みのデータの後ろに0を足し、FFTの出力を細かな周波数間隔で得る方法です。SciPy 1.14.1の公式解説と関数仕様をもとに、計算の刻み、観測時間、窓関数がそれぞれ何を変えるかを見ていきます。

    簡単にまとめると?

    • ゼロ埋めは、同じ記録からスペクトルを読む周波数の刻みを細かくします。新しい観測値は増えません。
    • 近い周波数の成分を見分けるには、観測時間と窓関数が決める山の幅、雑音、成分の強さの差を確かめます。
    • ピーク位置を読むためのFFT点数と、振幅・パワーをそろえるための正規化は、別々に設定します。

    FFTの点数と、実際に記録した点数を分ける

    FFTは高速フーリエ変換の略で、離散フーリエ変換(DFT)を効率よく計算する方法です。時間順に並んだ信号を、異なる周波数の波がどれだけ含まれるかという表現へ変換します。ここでは等間隔で取得したデータを扱います。[1][2]

    1秒あたりに記録する点数をサンプリング周波数 fs、実データの点数を N、FFTの計算点数を M と置きます。ゼロ埋めをしなければ M=N です。M>N にすると、足りない部分に0を並べてから変換します。SciPyの fft(x, n=M) はこの動作をします。逆に入力より小さい n を指定すると、データが切り詰められる点には注意が必要です。[2]

    出力の周波数の刻みは Δf=fs/M です。一方、実データに対応する観測区間の長さは τ=N/fs。この二つは、ゼロ埋めを始めると同じ点数では決まりません。ここでの区間長は公式解説と同じく「点数×サンプル間隔」で表し、最初と最後のサンプル時刻の差 (N−1)/fs とは区別します。[1]

    fs=1,000 Hzとした計算例。Cだけは、同じ現象を実際に4秒間記録し直す。
    条件 実データ N 観測区間 τ FFT点数 M 刻み Δf
    A:そのまま変換 1,000点 1秒 1,000点 1 Hz
    B:0を3,000点追加 1,000点 1秒 4,000点 0.25 Hz
    C:観測を4秒へ延長 4,000点 4秒 4,000点 0.25 Hz

    BとCは、どちらも0.25 Hz刻みの結果を返します。しかしBの後半は作った0であり、Cの後半は新たな観測値です。FFT点数が同じというだけでは、同じ情報を持つ測定にはなりません。この表は関係式を用いた説明例で、センサーによる実測結果ではありません。

    山の間を細かく読んでも、山そのものは狭くならない

    記録を1秒で切り出すことは、その外側を見ない「窓」を信号に掛けることに相当します。区間内をそのまま使う矩形窓では、単一の周波数成分も、周波数軸上で中央の山と周辺の小さな山を持つ形になります。中央を主ローブ、周辺をサイドローブと呼びます。有限時間で切り出したために他の周波数へ広がって見える現象が、スペクトル漏れです。[1]

    ゼロ埋めは、その広がりを持つスペクトルを、より多くの周波数位置で計算し直します。追加した0の項は和の値に寄与しないので、元の記録が作る周波数ごとの形自体を変えるわけではありません。粗い刻みでは飛び越えていた頂点付近を読めることが、見た目が滑らかになる理由です。

    幅を考えるための単純な例として、雑音のない一つの複素正弦波、つまり振動の位相を複素数で表した理想的な単一周波数の信号を、矩形窓で等間隔に切り出します。この場合、中央の頂点から最初の零点までの距離は fs/N=1/τ、左右の零点間の幅は 2/τ です。先ほどのAとBでは幅がともに2 Hz、Cでは0.5 Hzになります。Bの刻みは細かくても、窓が決める幅はAと同じです。これは窓の周波数応答から得られる例であり、実信号のあらゆるピーク形状をこの幅に固定する説明ではありません。

    「分解能」という言葉を、ソフトウェアによっては周波数の刻みの意味で使うこともあります。ここでは、Δfは結果を並べる刻み、二つの近い成分を見分ける能力は別の判断として扱います。主ローブ幅はその判断材料ですが、一定の間隔を超えれば必ず二つに分かれるという万能の境界ではありません。窓の種類、二成分の強さと位相、雑音、採用する推定方法によって見え方は変わります。

    公式図が示すのは、二つの別々の信号の比較

    次の原図は、SciPyの公式解説にある、20 Hzと20.5 Hzの正弦波をそれぞれ別に調べた例です。両方とも100 Hzで100点、1秒間の記録を使い、Hann窓を掛けています。菱形の点はゼロ埋めなし、細かな点をつないだ線はFFT点数を3倍にした結果です。[1][4]

    SciPy公式図。左は20 Hz、右は20.5 Hzの別々の正弦波。Hann窓を掛けた同じ100点の記録を、100点FFTの菱形と300点FFTの細かな点で比較している。
    図1 Hann窓を使った振幅スペクトルと3倍のゼロ埋め。左右は合成した二成分ではなく、別々の入力。原図の縦軸は片側を2倍しない表現で、単位振幅の正弦波の頂点は約0.5となる。SciPy 1.14.1 User Guide「Spectral Analysis」より、原図を無改変で掲載。Copyright (c) 2001–2002 Enthought, Inc.; 2003–2024 SciPy Developers. BSD 3-Clause、末尾のライセンス表示を参照。

    右側では、1 Hz刻みの菱形だけを見ると、20.5 Hzの頂点を直接読めません。1/3 Hz刻みにすると頂点の近くを捉えやすくなります。ただし20.5 Hzはその刻みにも厳密には乗らないため、細かくしたから頂点そのものを必ず計算したことになるわけでもありません。

    この図は「20 Hzと20.5 Hzが同時に含まれる信号を分離できた」という実験ではありません。また、Hann窓はサイドローブを抑える一方で主ローブを広げるので、矩形窓と比べた山の幅の変化を、ゼロ埋めの効果だけに帰すこともできません。図の入力、窓、FFT点数を分けて読むと、滑らかな山を過大評価しにくくなります。

    振幅を比べるときは、0の数で割らない

    FFTの点数を増やした後、振幅まで小さくなったように見えることがあります。そこで確認したいのが正規化、つまり計算結果を何で割って表示するかです。SciPyの fft の既定値 norm="backward" は、順方向の変換では点数による割り算を行いません。[2]

    矩形窓で得た同じ記録の振幅を比べるなら、ゼロ埋め後の M を無条件に分母にしないことが大切です。N 点の元の記録を4倍のFFTにしたとき、同じ周波数の変換値を M=4N で割れば、N で割った表示の4分の1になります。これは信号の振幅が減った証拠にはなりません。

    窓を掛けた振幅表示では、公式解説は窓の重みの和の絶対値 |Σw| での正規化を説明しています。図1の描画コードも、ゼロ埋め前後で同じ値を分母にしています。ただし、有限の刻みが真の頂点を外すことによる誤差まで、正規化だけで消えるわけではありません。実信号の片側スペクトルを2倍するかどうかも、別の表示上の約束です。[1][4]

    さらに、振幅の表示とパワースペクトル密度(PSD)は異なります。SciPyの periodogram では、入力がV、サンプリング周波数がHzなら、scaling="density" はV²/Hz、scaling="spectrum" はV²となります。窓や表示方式を変えたグラフの縦軸を、そのまま同じ「強さ」として比較しないようにします。[3]

    何を改善したいかで、変える設定を選ぶ

    一つの山の頂点付近を細かく読みたいなら、ゼロ埋めは役に立ちます。一方、隣り合う成分が重なって見えるなら、まず観測区間を長くできるか、窓の主ローブ幅とサイドローブが目的に合うかを確認します。長く記録する場合も、その間に回転数や振動状態が変われば、一つの定常的なスペクトルとしてまとめる判断そのものを見直す必要があります。

    サンプリング周波数を上げるときも、何を固定するかで結果が変わります。ゼロ埋めなしで、1秒間の記録を1,000 Hzから4,000 Hzへ変えるなら、実データは1,000点から4,000点へ増えますが、周波数の刻みはどちらも1 Hzです。反対に1,000点という点数だけを固定して4,000 Hzで記録すると、区間は0.25秒、刻みは4 Hzになります。時間軸を細かく刻むことと、周波数軸を細かく読むことは、同じ操作ではありません。

    たとえば測定ソフトの画面に「4,000点」とだけ表示されていても、それが実データ数なのか、0を追加した後の計算点数なのかは区別できません。取得条件と解析条件を別々に確認してから比較すると、設定値だけが増えた変更を、測定能力の向上と取り違えずに済みます。

    高い周波数まで観測したい場合には、FFT点数よりもサンプリング周波数と取得前の帯域制限が関係します。実信号の折り返しを避ける条件として、公式解説はナイキスト周波数 fs/2 を示しています。ゼロ埋めをしても、取得時に失った周波数情報や折り返しを元に戻すことはできません。[1]

    測定結果を比較するときは、実データ点数、サンプリング周波数、観測区間、窓、FFT点数、縦軸の単位と正規化を一緒に残す。この組み合わせが分かれば、「細かく見えるようになった」のか、「新しい観測によって見分けられるようになった」のかを、同じ条件で考えられます。

    参考資料

    参照版はSciPy 1.14.1。原図は公式資料の計算例であり、本文の比較表・零点幅・正規化の説明は関係式を用いた例です。実機測定やSciPyの性能ベンチマークは行っていません。

    1. SciPy 1.14.1 User Guide「Signal Processing — Spectral Analysis」。Continuous-time Sine Signal、Sampled Sine Signal、窓と正規化の節。
    2. SciPy 1.14.1「scipy.fft.fft」。n、norm、DFTの定義。
    3. SciPy 1.14.1「scipy.signal.periodogram」。nfft、scaling、単位。
    4. SciPy v1.14.1、図1の描画ソース。同一Hann窓・同一の振幅正規化と、100点/300点の計算条件。
    図1の著作権・ライセンス表示(BSD 3-Clause)

    SciPy v1.14.1 LICENSE.txt

    Copyright (c) 2001-2002 Enthought, Inc. 2003-2024, SciPy Developers.
    All rights reserved.
    
    Redistribution and use in source and binary forms, with or without
    modification, are permitted provided that the following conditions
    are met:
    
    1. Redistributions of source code must retain the above copyright
       notice, this list of conditions and the following disclaimer.
    
    2. Redistributions in binary form must reproduce the above
       copyright notice, this list of conditions and the following
       disclaimer in the documentation and/or other materials provided
       with the distribution.
    
    3. Neither the name of the copyright holder nor the names of its
       contributors may be used to endorse or promote products derived
       from this software without specific prior written permission.
    
    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.