実践に即した周波数領域解析の紹介

ライブスクリプトを開く

この例では、基本的な周波数領域信号解析の実行と解釈の方法を示します。ここでは、信号の周波数領域表現の利点と時間領域表現の利点について説明します。また、シミュレーションデータと実際のデータを使用して基本的な概念を説明します。この例は、次のような基本的な質問に答えます。FFT の振幅や位相とは何を意味しているのか、利用する信号が周期的であるか、パワーはどのように測定するのか、この帯域にある信号が 1 つかそれとも 2 つ以上か。

周波数領域解析は、信号処理アプリケーションにおける最も重要なツールです。周波数領域解析は、通信、地質学、リモートセンシング、イメージ処理などの領域で広く使用されています。時間領域解析では信号が時間の経過に伴ってどのように変化するかが示されますが、周波数領域解析では、信号のエネルギーが周波数のある範囲にわたってどのように分布しているかが示されます。周波数領域表現には位相シフトに関する情報も含まれます。この情報は、それぞれの周波数成分のすべてを組み合わせて元の時間信号を復元するために、各周波数成分に適用されなければなりません。

信号は、変換と呼ばれる数学演算子のペアを使用して時間領域、周波数領域間で変換できます。例として、フーリエ変換は、関数を多数の正弦波周波数成分の (無限) 和に分解します。周波数成分のスペクトルは信号の周波数領域表現です。逆フーリエ変換は周波数領域関数を時間の関数に戻します。MATLAB® では、fft関数およびifft関数によって、信号の離散フーリエ変換 (DFT) とその逆変換をそれぞれ計算できます。

FFT の振幅および位相の情報

信号の周波数領域表現は、信号の各周波数での振幅および位相の情報をもっています。このため、FFT 演算の出力は複素数になります。複素数 $x$ には実数部 $x_{r}$ と虚数部 $x_{i}$ があり、 $x = x_{r} + j x_{i}$ で表されます。 $x$ の振幅は、 $\sqrt{(x_{r}^{2} + x_{i}^{2})}$ で計算され、 $x$ の位相は、 $\arctan (x_{i} / x_{r})$ で計算されます。複素数の振幅と位相はそれぞれ MATLAB 関数absおよびangleを使用して求められます。

オーディオ信号の例を使って、信号の振幅と位相がどのような情報を保持しているかを理解します。まず、15 秒間のアコースティックギター音楽を含むオーディオファイルを読み込みます。オーディオ信号のサンプルレートは 44.1 kHz です。

Fs = 44100;
y = audioread("guitartune.wav");

fftを使用して、信号の周波数成分を調べます。

NFFT = length(y);
Y = fft(y,NFFT);
F = ((0:1/NFFT:1-1/NFFT)*Fs).';

FFT 出力は複素数ベクトルで、信号の周波数成分の情報が含まれます。振幅は、他の周波数成分に対する成分の強度を示します。位相は、周波数成分のすべてが時間領域でどのように並んでいるかを示します。

信号の周波数スペクトルの振幅成分と位相成分をプロットします。振幅は、便宜上、対数スケール (dB) でプロットします。位相は、unwrap関数を使用してアンラップし、周波数の連続関数を見ることができるようにします。

magnitudeY = abs(Y);        % FFT magnitude
phaseY = unwrap(angle(Y));  % FFT phase

helperFrequencyAnalysisPlot1(F,magnitudeY,phaseY,NFFT)

Figure contains 2 axes objects. Axes object 1 with title Magnitude response of the audio signal, xlabel Frequency in kHz, ylabel dB contains an object of type line. Axes object 2 with title Phase response of the audio signal, xlabel Frequency in kHz, ylabel radians contains an object of type line.

周波数領域ベクトル Y を逆フーリエ変換して時間信号を復元することができます。"symmetric" フラグは、ifftに対して実数値の時間信号を扱っているということを示し、計算時の数値誤差によって逆変換に現れる小さな虚数部の値は 0 に設定されます。元の時間信号 y と復元された信号 y1 が実質的には同じであることがわかります (差のノルムは約 1e-14 です)。この 2 つの値の間の非常に小さい差も上述の数値誤差によるものです。逆フーリエ変換された信号 y1 を再生して聴きます。

y1 = ifft(Y,NFFT,"symmetric");
norm(y-y1)

ans = 
3.7900e-14

hplayer = audioplayer(y1,Fs);
play(hplayer);

信号の振幅応答の変化の影響を確認するため、振幅を 0 にして、1 kHz を超える周波数成分を直接 FFT 出力から除去します。それから、音を聴いてオーディオファイルの音に与える影響を確認します。信号の高周波数成分を取り除くことをローパスフィルター処理と呼びます。

Ylp = Y;
Ylp(F>=1000 & F<=Fs-1000) = 0;

helperFrequencyAnalysisPlot1(F,abs(Ylp),unwrap(angle(Ylp)), ...
    NFFT,"Frequency components above 1 kHz have been zeroed")

Figure contains 2 axes objects. Axes object 1 with title Magnitude response of the audio signal Frequency components above 1 kHz have been zeroed, xlabel Frequency in kHz, ylabel dB contains an object of type line. Axes object 2 with title Phase response of the audio signal, xlabel Frequency in kHz, ylabel radians contains an object of type line.

ifftを使用して、フィルター処理を行った信号を時間領域に戻します。

ylp = ifft(Ylp,"symmetric");

信号を再生します。メロディーは聞こえますが、耳をふさいだときのように聞こえます (耳をふさぐと高い周波数の音は排除されます)。弦で音を出すとき、ギターは 400 Hz ～ 1 kHzの音を出しますが、弦は基本周波数の倍数でも振動します。これらの高い周波数成分は高調波と呼ばれ、これがギターにその特有の音色を与えています。高調波を除去すると、音が "こもり" ます。

hplayer = audioplayer(ylp, Fs);
play(hplayer);

信号の位相は、時間領域で曲の音が現れるタイミングに関する重要な情報をもっています。オーディオ信号の位相の重要性を説明するために、各周波数成分の振幅を取ることによって位相情報を完全に除去します。これによって振幅応答は変化しないことに注目してください。

% Take the magnitude of each FFT component of the signal
Yzp = abs(Y);
helperFrequencyAnalysisPlot1(F,abs(Yzp), ...
    unwrap(angle(Yzp)),NFFT,[],"Phase has been set to zero")

Figure contains 2 axes objects. Axes object 1 with title Magnitude response of the audio signal, xlabel Frequency in kHz, ylabel dB contains an object of type line. Axes object 2 with title Phase response of the audio signal Phase has been set to zero, xlabel Frequency in kHz, ylabel radians contains an object of type line.

信号を時間領域に戻し、オーディオ信号を再生します。元の音とはまったく異なる音になっています。今回は振幅応答が同じで、周波数成分は何も除去されていません。しかし、音の順序がすっかり消えています。ここでの信号は、時刻 0 ですべてが整列した正弦波のグループで構成されています。通常、フィルター処理によって生じる位相歪みは、それと認識できない状態になるほど信号を損ないます。

yzp = ifft(Yzp,"symmetric");
hplayer = audioplayer(yzp, Fs);
play(hplayer);

信号の周期性の検出

信号の周波数領域表現によって、時間領域で信号を見たときには確認するのが容易でないか、まったく見えないいくつかの信号特性を調べることができます。たとえば、周波数領域解析は信号の周期的な性質を見つけるときに役立ちます。

オフィスビルにおける温度の周期的な性質の解析

冬期のオフィスビルにおける温度測定の結果を見てみましょう。測定は約 16.5 週間、30 分間隔で行われました。時間軸が週にスケーリングされた時間領域データを見てみましょう。このデータに何らかの周期的な性質はあるでしょうか。

load officetemp.mat
Fs = 1/(60*30);   % Sample rate is 1 sample every 30 minutes
t = (0:length(temp)-1)/Fs;

helperFrequencyAnalysisPlot2(t/(60*60*24*7),temp, ...
  "Time (weeks)","Temperature (\circF)")

Figure contains an axes object. The axes object with xlabel Time (weeks), ylabel Temperature ( degree F) contains an object of type line.

時間領域信号を見て、オフィスの温度に何らかの周期的な性質があるどうかを知ることはほとんど不可能です。しかし、温度の周波数領域表現を見れば、その周期的な性質は明らかです。

信号の周波数領域表現を求めます。周波数軸が回/週にスケーリングされた FFT 出力の振幅をプロットすると、他の周波数成分より明らかに大きな 2 つのスペクトル線があることがわかります。1 つのスペクトル線は 1 回/週、もう 1 つのスペクトル線は 7 回/週にあります。7 日単位の日程表により温度制御されているビルのデータであるということを考えれば、これは理にかなっています。最初のスペクトル線は、ビルの温度が週末は温度が低く、平日は高いという週次サイクルに従うことを示しています。2 番目のスペクトル線は、1 日ごとの周期性もあり、夜間は温度が低く、日中は高いということを示しています。

NFFT = length(temp);              % Number of FFT points
F = (0 : 1/NFFT : 1/2-1/NFFT)*Fs; % Frequency vector

TEMP = fft(temp,NFFT);            
TEMP(1) = 0; % remove the DC component for better visualization

helperFrequencyAnalysisPlot2(F*60*60*24*7,abs(TEMP(1:NFFT/2)), ...
  "Frequency (cycles/week)","Magnitude",[],[],[0 10])

Figure contains an axes object. The axes object with xlabel Frequency (cycles/week), ylabel Magnitude contains an object of type line.

パワーの測定

periodogram関数は信号の FFT を計算し、その出力を正規化してパワースペクトル密度 PSD、またはパワーの測定元となるパワースペクトルを求めます。PSD は、時間信号のパワーが周波数についてどのように分布しているかを示します。単位はワット/Hz です。PSD の各点を、その点が定義された周波数範囲にわたって (つまり、PSD の分解能帯域幅にわたって) 積分して、パワースペクトルを計算します。パワースペクトルの単位はワットです。パワーの値は、周波数範囲にわたって積分せずに、直接パワースペクトルから読み取ることができます。PSD もパワースペクトルも実数なので、位相情報は何も含まれていないことに注意してください。

非線形パワー増幅器の出力での高調波測定

パワー増幅器の出力で測定したデータを読み込みます。出力は 3 次の歪みをもつ $v_{o} = v_{i} + 0.75 v_{i}^{2} + 0.5 v_{i}^{3}$ という形で表されます。ここで、 $v_{o}$ は出力電圧で、 $v_{i}$ は入力電圧です。データは、サンプルレート 3.6 kHz で測定されました。入力 $v_{i}$ は大きさ 1 の振幅をもつ 60 Hz の正弦波です。非線形歪みの性質から、増幅器の出力信号は DC 成分、60 Hz 成分および 120 Hz と 180 Hz の第 2、第 3 高調波を含みます。

増幅器出力を 3600 サンプル読み込み、パワースペクトルを計算し、結果を対数スケール (dBW) でプロットします。

load ampoutput1.mat
Fs = 3600;
NFFT = length(y);

% Power spectrum is computed when you pass a "power" flag input
[P,F] = periodogram(y,[],NFFT,Fs,"power");   

helperFrequencyAnalysisPlot2(F,10*log10(P),"Frequency (Hz)", ...
  "Power spectrum (dBW)",[],[],[-0.5 200])

Figure contains an axes object. The axes object with xlabel Frequency (Hz), ylabel Power spectrum (dBW) contains an object of type line.

パワースペクトルのプロットには、予想される 4 つのピークのうちの 3 つである、DC 成分、60 Hz 成分、120 Hz 成分があります。他にも、信号に含まれるノイズによるスプリアスピークがいくつかあります。180 Hz の高調波はすっかりノイズに埋もれていることに注目してください。

観測できるピークのパワーを測定します。

PdBW = 10*log10(P);
power_at_DC_dBW = PdBW(F==0)

power_at_DC_dBW = 
-7.8873

[peakPowers_dBW, peakFreqIdx] = findpeaks(PdBW,MinPeakHeight=-11);
peakFreqs_Hz = F(peakFreqIdx)

peakFreqs_Hz = 2×1

    60
   120

peakPowers_dBW

peakPowers_dBW = 2×1

   -0.3175
  -10.2547

ノイズを含む信号のパワー測定の改良

上のプロットで見られるように、ピリオドグラムには目的の信号と関係のない周波数ピークがいくつか表示されています。スペクトルにはノイズがかなり含まれているように見えます。これは、ノイズを含む信号について短時間の実現を 1 回しか解析していないからです。数回実験を繰り返して平均をとると、スプリアススペクトルピークが除去され、より精度の高いパワー測定値になります。pwelch関数を使用すると、この平均化を行うことができます。この関数は、大きなデータベクトルを指定した長さの小さなセグメントに分けて、セグメント数分のピリオドグラムを計算し、それらの平均値を求めます。有効セグメント数が多いほど、pwelch関数はパワースペクトルをより平滑化 (分散を小さく) し、パワーの値は期待値に近づきます。

500e3 点の増幅器出力から成る、より大きな測定データを読み込みます。FFT 実行に使用するデータ点の数を同じ 3600 として、floor(500e3/3600) = 138 点の FFT を平均してパワースペクトルを求めます。

load ampoutput2.mat
SegmentLength = NFFT; 

% Power spectrum is computed when you pass a "power" flag input
[P,F] = pwelch(y,ones(SegmentLength,1),0,NFFT,Fs,"power");

helperFrequencyAnalysisPlot2(F,10*log10(P),"Frequency (Hz)", ...
  "Power spectrum (dBW)",[],[],[-0.5 200])

Figure contains an axes object. The axes object with xlabel Frequency (Hz), ylabel Power spectrum (dBW) contains an object of type line.

プロットを見ると、pwelchはノイズに起因するすべてのスプリアス周波数ピークを効果的に除去しています。ノイズに埋もれていた 180 Hz のスペクトル成分が今は見えています。平均化によってスペクトルの分散をなくし、より精度の高いパワー測定値を得ることができます。

平均パワー合計および周波数帯域パワーの測定

時間領域信号の平均パワーの合計の測定は、簡単な、通常タスクです。増幅器の出力信号 y の平均パワー合計は時間領域で次のように計算されます。

pwr = sum(y.^2)/length(y) % in watts

pwr = 
8.1697

周波数領域では、平均パワー合計は信号のすべての周波数成分のパワーの和として計算されます。pwr1 の値は、信号のパワースペクトル内において有効な周波数成分すべての和で構成されます。この値は、時間領域の信号を使用して計算された上述の pwr の値と一致します。

pwr1 = sum(P) % in watts

pwr1 = 
8.1698

では、周波数帯域にわたって有効なパワーの合計を測定する場合はどうでしょうか。この場合は、bandpower関数を使用して、目的の周波数帯域にわたるパワーを計算することができます。時間領域信号を入力として直接この関数に渡すことで、指定した周波数帯域全体のパワーが求められます。この場合、関数はピリオドグラム法でパワースペクトルを推定します。

50 Hz ～ 70 Hz 帯域のパワーを計算します。この結果は、60 Hz のパワーに対象の帯域にわたるノイズパワーを加算したものになります。

pwr_band = bandpower(y,Fs,[50 70]);
pwr_band_dBW = 10*log10(pwr_band) % dBW

pwr_band_dBW = 
0.0341

帯域内のパワーの測定に使用されるパワースペクトルの計算を調整する場合は、PSD ベクトルをbandpower関数に渡します。たとえば、前と同じようにpwelch関数を使用して PSD を計算し、ノイズの影響を確実に平均化することができます。

% Power spectral density is computed when you specify the "psd" option
[PSD,F]  = pwelch(y,ones(SegmentLength,1),0,NFFT,Fs,"psd");
pwr_band1 = bandpower(PSD,F,[50 70],"psd");
pwr_band_dBW1 = 10*log10(pwr_band1) % dBW

pwr_band_dBW1 = 
0.0798

スペクトル成分の検出

信号が 1 以上の周波数成分を含んでいることがあります。すべてのスペクトル成分を観測できるかどうかは、解析の周波数分解能によって決まります。パワースペクトルの周波数分解能または分解能帯域幅は、R = Fs/N で定義されます。ここで、N は信号観測の長さです。周波数分解能より大きい周波数で隔てられているスペクトル成分だけが分解されます。

建物の地震振動制御システムの解析

アクティブマスドライバー (AMD) 制御システムは、地震発生時の建物の振動を減らすために使用されます。アクティブマスドライバーは建物の最上階に設置され、建物の床の変位と加速度の測定値に基づいて動作します。制御システムは信号をドライバーに送り、地面の揺れを吸収するように錘を動かします。加速度の測定値は、地震の条件下で 3 階建ての試験構造物の 1 階での記録です。値は、アクティブマスドライバー制御システムなし (開ループ条件) およびアクティブマスドライバー制御システムあり (閉ループ条件) で測定されました。

加速度データを読み込み、1 階の床における加速度のパワースペクトルを計算します。データベクトルの長さは 10e3、サンプルレートは 1 kHz です。長さが 64 データ点のセグメントにpwelchを使用し、floor(10e3/64) = 156 点の FFT の平均値と、分解能帯域幅 Fs/64 = 15.625 Hz を求めます。前述のように、平均化によってノイズの影響が減り、より精度の高いパワー測定値が得られます。FFT 点の数を 512 にします。NFFT > N を使用すると、効果的に周波数点が内挿され、より細かいスペクトルプロットが得られます (これは、時間信号の最後に NFFT-N 個の 0 を追加し、ゼロで埋められたベクトルの NFFT 点の FFT を実行することによって実現できます)。

開ループ時および閉ループ時の加速度のパワースペクトルから、制御システムありの場合、加速度のパワースペクトルが 4 ～ 11 dB 減少していることが分かります。約 23.44 kHz における減衰が最大です。11 dB の減少は、振動のパワーを 12.6 の倍率で減少していることを意味します。合計パワーは 0.1670 ワットから 0.059 ワットに減り、これは 2.83 の倍率での減少となります。

load quakevibration.mat

Fs = 1e3;                 % sample rate
NFFT = 512;               % number of FFT points
segmentLength = 64;       % segment length

% open loop acceleration power spectrum
[P1_OL,F] = pwelch(gfloor1OL,ones(segmentLength,1),0,NFFT,Fs,"power"); 

% closed loop acceleration power spectrum
P1_CL     = pwelch(gfloor1CL,ones(segmentLength,1),0,NFFT,Fs,"power"); 

helperFrequencyAnalysisPlot2(F,10*log10([(P1_OL) (P1_CL)]),...
  "Frequency (Hz)","Acceleration Power Spectrum (dB)",...
  "Resolution bandwidth = 15.625 Hz",["Open" "Closed"]+" loop",[0 100])

Figure contains an axes object. The axes object with title Resolution bandwidth = 15.625 Hz, xlabel Frequency (Hz), ylabel Acceleration Power Spectrum (dB) contains 2 objects of type line. These objects represent Open loop, Closed loop.

現在、振動データの解析をしていて、振動が周期的な性質をもつことはわかっています。では、周期的動作特有のシャープなスペクトル線が、上図のスペクトルプロットに 1 つもないのはどういうことでしょうか。これらの線が 64 点のセグメント長で得られる分解能では分解できず、見逃されている可能性があります。周波数分解能を増やし、前に分解できなかったスペクトル線があるかどうか確認します。そのために、pwelch関数で使用するデータセグメント長を 512 点に増やします。これにより、新しい分解能は Fs/512 = 1.9531 Hz になります。この場合、FFT 平均の数は floor(10e3/512) = 19 に減ります。明らかに、pwelch使用時には FFT 平均の数と周波数分解能とのトレードオフが存在します。FFT 点の数は同じ 512 にします。

NFFT = 512;            % number of FFT points
segmentLength = 512;   % segment length

[P1_OL,F] = pwelch(gfloor1OL,ones(segmentLength,1),0,NFFT,Fs,"power"); 
P1_CL     = pwelch(gfloor1CL,ones(segmentLength,1),0,NFFT,Fs,"power"); 

helperFrequencyAnalysisPlot2(F,10*log10([(P1_OL) (P1_CL)]),  ...
  "Frequency in Hz","Acceleration Power Spectrum in dB", ...
  "Resolution bandwidth = 1.95 Hz",{"Open loop", "Closed loop"},[0 100])

Figure contains an axes object. The axes object with title Resolution bandwidth = 1.95 Hz, xlabel Frequency in Hz, ylabel Acceleration Power Spectrum in dB contains 2 objects of type line. These objects represent Open loop, Closed loop.

周波数分解能を増やすことによって、開ループスペクトルで 3 つのピークを、閉ループスペクトルで 2 つのピークを観測できることに注目してください。これらのピークは、以前は分解できませんでした。開ループスペクトルにおけるピーク間の間隔は約 11 Hz で、セグメントの長さ 64 で得られる周波数分解能より小さく、セグメントの長さ 512 で得られる周波数分解能より大きくなっています。ここでは、振動の周期的な振る舞いが見られます。主要な振動周波数は 5.86 Hz で、等間隔の周波数ピークはこれらが調和関係にあることを示しています。制御システムが振動の全体的なパワーを減らすことは既に確認済みですが、分解能の高いスペクトルを見るとわかるように、制御システムのもう 1 つの効果は 17.58 Hz で高調波成分を記録することです。制御システムは振動を減らすだけでなく、振動を正弦波に近づけています。

ここで重要な点は、周波数分解能が FFT 点の数ではなく、信号点の数で決まることです。FFT 点の数を増やすと、周波数データが内挿され、より細かいスペクトルが得られますが、分解能を改善するわけではありません。

まとめ

この例では、fft、ifft、periodogram、pwelch、およびbandpowerの各関数を使用して、信号の周波数領域解析を行う方法を説明しました。FFT の複雑な性質および周波数スペクトルの振幅および位相に含まれる情報が何であるかを示しました。信号の周期性を解析するときに周波数領域データを使用する利点を確認しました。合計パワーおよびノイズを含む信号の特定の周波数帯域全体のパワーの計算方法を紹介しました。スペクトルの周波数分解能を増やすことによって間隔が狭い周波数成分を観測できること、周波数分解能とスペクトルの平均化との間にトレードオフがあることを学びました。

参考情報

周波数領域解析の詳細については、スペクトル解析を参照してください。

参考文献: J.G. Proakis and D. G. Manolakis, "Digital Signal Processing. Principles, Algorithms, and Applications", Prentice Hall, 1996.

付録

この例では、以下の補助関数を使用します。

参考

fft | periodogram | pspectrum | pwelch