Main Content

このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。

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

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

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

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

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

信号の周波数領域表現は、信号の各周波数での振幅および位相の情報をもっています。このため、FFT 演算の出力は複素数になります。複素数 $x$ には実数部 $x_r$ と虚数部 $x_i$ があり、$x = x_r + ix_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);        % Magnitude of the FFT
phaseY = unwrap(angle(Y));  % Phase of the FFT

helperFrequencyAnalysisPlot1(F,magnitudeY,phaseY,NFFT)

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

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

   3.7851e-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')

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')

信号を時間領域に戻し、オーディオ信号を再生します。元の音とはまったく異なる音になっています。今回は振幅応答が同じで、周波数成分は何も除去されていません。しかし、音の順序がすっかり消えています。ここでの信号は、時刻 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 in weeks','Temperature (Fahrenheit)')

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

信号の周波数領域表現を求めます。周波数軸が回/週にスケーリングされた 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])

パワーの測定

関数 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 in Hz',...
  'Power spectrum (dBW)',[],[],[-0.5 200])

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

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

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

[peakPowers_dBW, peakFreqIdx] = findpeaks(PdBW,'minpeakheight',-11);
peakFreqs_Hz = F(peakFreqIdx)
peakPowers_dBW
power_at_DC_dBW =

   -7.8873


peakFreqs_Hz =

    60
   120


peakPowers_dBW =

   -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 in Hz',...
  'Power spectrum (dBW)',[],[],[-0.5 200])

プロットを見ると、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 in Hz','Acceleration Power Spectrum in dB',...
  'Resolution bandwidth = 15.625 Hz',{'Open loop', 'Closed loop'},[0 100])

現在、振動データの解析をしていて、振動が周期的な性質をもつことはわかっています。では、周期的動作特有のシャープなスペクトル線が、上図のスペクトル プロットに 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])

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

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

まとめ

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

参考情報

周波数領域解析の詳細は、Signal Processing Toolbox を参照してください。

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

付録

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

参考

| | |