MATLABで波形データを周波数分析！フーリエ変換の基本をわかりやすく

By /

はい、承知いたしました。
MATLABを用いたフーリエ変換による周波数分析について、基本から応用までを網羅した詳細な解説記事を作成します。約5000語のボリュームで、初心者の方でもステップバイステップで理解を深められるように構成します。

MATLABで波形データを周波数分析！フーリエ変換の基本をわかりやすく徹底解説

1. はじめに

私たちの周りには、音、振動、電波、株価の変動など、時間とともに変化する「波形データ」が溢れています。これらのデータをただ眺めているだけでは、その背後に隠された本質的な特徴を見抜くことは困難です。例えば、ある機械から発生する異音の波形を見ただけでは、どの部品が異常な振動を引き起こしているのかを特定するのは難しいでしょう。

このような時間領域のデータを分析し、その中に「どのような周波数成分が、どれくらいの強さで含まれているか」を明らかにする強力なツールが周波数分析であり、その中核をなす数学的手法がフーリエ変換です。

この記事では、科学技術計算の分野で広く使われているソフトウェア「MATLAB」を使い、フーリエ変換による周波数分析の基本を、理論から実践まで徹底的に解説します。

対象読者

  • MATLABの基本的な操作（変数の作成、スクリプトの実行、簡単なプロットなど）ができる方
  • 信号処理やフーリエ変換については初心者、または知識を再確認したいエンジニア、学生、研究者の方
  • 手元の時系列データ（音、振動、センサーデータなど）を周波数分析してみたいと考えている方

なぜ周波数分析が重要なのか？

時間領域の波形は、いわば「素材そのままの料理」です。美味しいかもしれませんが、どのような食材（周波数成分）で構成されているかは分かりません。周波数分析は、この料理を「食材ごと」に分解し、どの食材がどれくらい使われているかを明らかにするプロセスに例えられます。

これにより、以下のようなことが可能になります。

  • 異音・異常振動の特定: 機械の正常時の周波数パターンと比較し、異常な周波数成分を検出することで故障診断を行う。
  • ノイズ除去: 信号に含まれる不要な周波数成分（ノイズ）を特定し、フィルタリングで除去する。
  • 音声認識: 人の声に含まれる周波数のパターン（フォルマント）を分析し、母音などを識別する。
  • 通信: 限られた周波数帯域を効率的に利用して、多くの情報を送受信する。

この記事を読むことで得られること

この記事を最後まで読めば、あなたは以下のスキルを習得できます。

  1. フーリエ変換の基本的な概念が、数式アレルギーの方でも直感的に理解できる。
  2. MATLABのfft関数を使い、具体的な手順に沿って周波数分析を実行できるようになる。
  3. 分析結果として得られる周波数スペクトルのグラフを正しく解釈し、意味のある情報を引き出せるようになる。
  4. 周波数分析で陥りがちな「エイリアシング」や「スペクトル漏れ」といった落とし穴を理解し、その対策を講じられるようになる。

さあ、MATLABと共に、波形データの隠された世界を探る旅に出かけましょう。

2. フーリエ変換とは？ (理論編)

実践的なコーディングに入る前に、フーリエ変換が一体何をしているのか、その基本的な考え方を理解しておくことが重要です。難しい数式は最小限にとどめ、直感的なイメージを掴むことを目指します。

直感的な理解：「波のレシピ」を暴く魔法

フーリエ変換の核心的なアイデアは、フランスの数学者ジョゼフ・フーリエが提唱した「あらゆる複雑な波形は、単純なサイン波（正弦波）の無限の足し合わせで表現できる」というものです。

これは、オーケストラの演奏に例えると分かりやすいかもしれません。

  • 時間領域の波形: 私たちが耳で聞くオーケストラの「演奏」そのものです。非常に複雑で豊かな音色に聞こえます。
  • 周波数領域のスペクトル: その演奏を「どの楽器（周波数）が、どれくらいの音量（振幅）で鳴っているか」という情報に分解したものです。「バイオリンのパートがこのくらいの音量で、チェロがこのくらい、トランペットが…」というように、構成要素を分析した結果が周波数スペクトルです。

フーリエ変換は、この「演奏（時間領域の波形）」から「各楽器の音量（周波数スペクトル）」を導き出す、いわば「波のレシピ」を暴く魔法なのです。


（概念図：左側に複雑な波形（時間領域）、右側にその波形を構成する複数のサイン波の周波数と振幅を示す棒グラフ（周波数領域）を描画）

図1: 時間領域と周波数領域の関係

左のグラフ（時間領域）では、横軸が時間、縦軸が信号の振幅です。一方、右のグラフ（周波数領域）では、横軸が周波数、縦軸がその周波数成分の強さ（振幅やパワー）を表します。フーリエ変換は、左のグラフを右のグラフに変換する操作です。

コンピュータが扱う離散フーリエ変換（DFT）と高速フーリエ変換（FFT）

本来のフーリエ変換（連続フーリエ変換）は、連続的な無限の長さの信号を対象としますが、コンピュータはアナログで連続的なデータを直接扱うことができません。コンピュータが扱うデータは、一定の時間間隔で値を測定（サンプリング）した、離散的で有限のデータです。

この離散的なデータに対してフーリエ変換を行うのが、離散フーリエ変換（Discrete Fourier Transform, DFT）です。

しかし、DFTはデータ点数が増えると計算量が爆発的に増加するという欠点がありました。そこで、このDFTの計算を劇的に高速化するアルゴリズムとして考案されたのが高速フーリエ変換（Fast Fourier Transform, FFT）です。

MATLABでfftという名前の関数を使うのはこのためです。私たちが実践で使うのは、このFFTアルゴリズムということになります。FFTはDFTと数学的に等価ですが、計算方法が非常に効率的である、と理解しておけば十分です。

なぜ結果が複素数になるのか？

fftを実行すると、結果が複素数の配列として返ってきます。これに戸惑う初心者の方は少なくありません。なぜ虚数単位 i が出てくるのでしょうか？

これは、フーリエ変換が各周波数成分の「振幅（強さ）」だけでなく「位相（タイミングのズレ）」の情報も同時に保持しているためです。

  • 振幅 (Amplitude): その周波数のサイン波がどれくらいの大きさか。
  • 位相 (Phase): そのサイン波が、基準となるタイミングからどれだけズレて始まっているか。

複素数は、a + bi の形で表され、複素平面上の点で表現できます。この点の原点からの距離が振幅に、実軸とのなす角度が位相に対応します。このように、一つの複素数で振幅と位相という2つの情報を同時に表現できるため、フーリエ変換の結果は複素数となるのです。

多くの周波数分析では、まずは「どの周波数成分が強いか」を知ることが目的なので、位相情報は一旦無視して振幅情報だけを取り出してプロットすることが一般的です。振幅は、複素数の絶対値 abs() を計算することで簡単に得られます。

理論はここまでです。次は、いよいよMATLABを使って実際に手を動かしていきましょう。

3. MATLABで周波数分析を実践してみよう (実践編)

ここからは、MATLABのコードを使いながら、周波数分析の具体的な手順をステップバイステップで学んでいきます。まずは、答えが分かっているシンプルな合成波形を分析対象とすることで、手法の正しさを確認しながら進めましょう。

Step 1: 分析対象の波形データを作成する

今回は、振幅0.7の50Hzのサイン波と、振幅1.0の120Hzのサイン波を足し合わせた合成波を作成します。最終的に、周波数分析の結果で50Hzと120Hzの2つのピークが正しく現れることを目指します。

まず、信号の基本的なパラメータを定義します。

matlab
% --- 信号のパラメータ設定 ---
Fs = 1000; % サンプリング周波数 (Hz)
T = 1/Fs; % サンプリング周期 (s)
L = 1500; % 信号の長さ (サンプル数)
t = (0:L-1)*T; % 時間ベクトル (s)

各パラメータの意味を解説します。

  • Fs (サンプリング周波数): 1秒間に何回データをサンプリング（測定）するかを示す値です。単位はHz。ここでは1000Hz、つまり1秒間に1000回のデータを取得することを意味します。この値は非常に重要で、分析したい周波数の2倍以上である必要があります（詳細は後述の「エイリアシング」で解説）。
  • T (サンプリング周期): 1回のサンプリングにかかる時間です。サンプリング周波数の逆数 (1/Fs) で計算できます。
  • L (信号の長さ): 合計で何個のデータ点を取得するか。データ長が長いほど、より細かい周波数分解能で分析できます。
  • t (時間ベクトル): 各データ点がいつの時点のものかを示す時間の配列です。0秒から始まり、T秒ずつ増えていきます。

次に、この時間ベクトル t を使って、2つのサイン波を合成した信号 S を作成します。

matlab
% --- 合成波形の作成 ---
% 振幅0.7, 周波数50Hzのサイン波と、振幅1.0, 周波数120Hzのサイン波を合成
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

sin(2*pi*f*t) が、周波数 f のサイン波を生成する基本的な数式です。

作成した波形がどのようなものか、プロットして確認してみましょう。

matlab
% --- 時間領域での波形プロット ---
figure; % 新しい図ウィンドウを作成
plot(t(1:100), S(1:100)); % データが多すぎると見づらいので最初の100点だけプロット
title('合成信号の時間領域波形');
xlabel('時間 (s)');
ylabel('振幅');
grid on;

これを実行すると、2つの波が混ざり合った複雑な波形が表示されるはずです。このグラフを見ただけでは、50Hzと120Hzの成分が含まれていることを見抜くのは困難です。


（生成されるグラフのイメージ：複雑に振動する波形）

図2: 作成した合成信号の時間領域波形

Step 2: FFTを実行する

いよいよ、この時間領域の信号 S をフーリエ変換します。MATLABでは、fft 関数を一行呼び出すだけで完了です。

matlab
% --- FFTの実行 ---
Y = fft(S);

これだけです。驚くほど簡単ですね。
先述の通り、この時点での Y は、L個の要素を持つ複素数の配列になっています。このままでは人間が解釈するのは難しいので、次のステップで意味のある形に整形していきます。

Step 3: FFTの結果を解釈できる形に整形する

ここが周波数分析で最も重要かつ、初心者がつまずきやすいポイントです。FFTの生の出力結果から、私たちが解釈したい「片側振幅スペクトル」を計算する手順を丁寧に見ていきましょう。

FFTの結果 Y は、両側スペクトル (Two-sided spectrum) と呼ばれる形式になっています。これは、正の周波数と負の周波数の両方の情報を含んでおり、通常は左右対称な形をしています。しかし、実信号を分析する場合、負の周波数の情報は正の周波数の情報と冗長なので、私たちは通常、正の周波数領域のみをプロットする片側スペクトル (Single-sided spectrum) を使います。

片側振幅スペクトル P1 を計算するコードは以下のようになります。

“`matlab
% — FFT結果の整形 (片側振幅スペクトルの計算) —

% 1. 複素数の絶対値を取って振幅を計算
P2 = abs(Y/L);

% 2. 片側スペクトルを計算 (データ点Lの半分まで)
P1 = P2(1:L/2+1);

% 3. 直流成分(0Hz)とナイキスト周波数成分を除き、振幅を2倍にする
P1(2:end-1) = 2*P1(2:end-1);
“`

このコードブロックを分解して解説します。

  1. P2 = abs(Y/L);

    • abs(Y): fft の結果 Y は複素数なので、abs() 関数でその絶対値（大きさ）を計算します。これにより、各周波数成分の振幅が得られます。
    • /L: なぜ信号長 L で割るのでしょうか？これは、FFTの出力の大きさが信号長に比例してしまうため、信号長に依存しない正しい振幅を得るために正規化する処理です。これにより、元のサイン波の振幅（0.7と1.0）とスケールが合います。この P2 が両側振幅スペクトルです。

  2. P1 = P2(1:L/2+1);

    • 両側スペクトル P2 は、0Hzからサンプリング周波数 Fs までの情報を持ちますが、Fs/2（ナイキスト周波数）を境に左右対称な形になっています。
    • そのため、意味のある情報は前半部分（0Hzから Fs/2 まで）に集約されています。ここでは、P2 の最初から L/2+1 番目の要素までを取り出して、片側スペクトル P1 を作成しています。

  3. P1(2:end-1) = 2*P1(2:end-1);

    • 両側スペクトルでは、各周波数成分のエネルギーが正の周波数と負の周波数に半分ずつ分配されています。
    • 片側スペクトルにする際に負の周波数側を捨てたため、失われたエネルギー分を補うために、振幅を2倍にする必要があります。
    • なぜ 2:end-1 なのか？:
      • 最初の要素 P1(1) は直流成分 (DC成分、0Hz) に対応します。これには負の周波数の対になる成分が存在しないため、2倍にはしません。
      • 最後の要素 P1(end) はナイキスト周波数 (Fs/2) に対応します。これも同様に対になる成分がない特別な点なので、2倍にはしません。（データ長 L が偶数の場合）
    • そのため、直流成分とナイキスト周波数成分を除く、間のすべての要素を2倍にしています。

この3つのステップを経て、ようやく解釈可能な片側振幅スペクトル P1 が完成しました。

Step 4: 周波数軸を作成し、結果をプロットする

スペクトルの縦軸（振幅 P1）は準備できましたが、横軸である周波数軸がまだありません。P1 の各要素が、具体的に何Hzに対応するのかを計算する必要があります。

周波数軸 f は以下のように作成できます。

matlab
% --- 周波数軸の作成 ---
f = Fs*(0:(L/2))/L;

この式の意味を考えてみましょう。
* FFTの結果は、Fs/L という周波数分解能で区切られています。つまり、スペクトル上の隣り合う点の間隔は Fs/L [Hz] です。
* 0:(L/2) は、0から L/2 までの整数のベクトルを作成します。
* これに Fs/L を掛けることで、0, Fs/L, 2*Fs/L, ..., Fs/2 という、P1 の各点に対応する周波数のベクトル f が得られます。

これで、横軸 f と縦軸 P1 が揃いました。いよいよプロットして結果を確認します。

matlab
% --- 片側振幅スペクトルのプロット ---
figure;
plot(f, P1);
title('片側振幅スペクトル');
xlabel('周波数 (Hz)');
ylabel('振幅 |P1(f)|');
grid on;
xlim([0 200]); % 見やすいように表示範囲を200Hzまでに制限

これを実行すると、以下のようなグラフが得られるはずです。


（生成されるグラフのイメージ：横軸が周波数、縦軸が振幅。50Hzの位置に高さ0.7のピーク、120Hzの位置に高さ1.0のピークが鋭く立っている）

図3: 合成信号の片側振幅スペクトル

見事に、50Hzの位置に振幅0.7のピーク、そして120Hzの位置に振幅1.0のピークがはっきりと現れました！
時間領域の波形（図2）では分からなかった信号の構成要素が、周波数分析によって一目瞭然となりました。これで、基本的な周波数分析の手順はマスターです。

全体のコードまとめ

ここまでの処理を一つのMATLABスクリプトにまとめます。

“`matlab
% MATLABでフーリエ変換による周波数分析 基本編

%% 1. パラメータ設定と波形生成
% 信号のパラメータ設定
Fs = 1000; % サンプリング周波数 (Hz)
T = 1/Fs; % サンプリング周期 (s)
L = 1500; % 信号の長さ (サンプル数)
t = (0:L-1)*T; % 時間ベクトル (s)

% 合成波形の作成
% 振幅0.7, 周波数50Hzのサイン波と、振幅1.0, 周波数120Hzのサイン波を合成
S = 0.7sin(2pi50t) + sin(2pi120*t);

% 時間領域での波形プロット
figure;
plot(t(1:100), S(1:100));
title(‘合成信号の時間領域波形’);
xlabel(‘時間 (s)’);
ylabel(‘振幅’);
grid on;

%% 2. FFTの実行と結果の整形
% FFTの実行
Y = fft(S);

% 両側振幅スペクトルの計算
P2 = abs(Y/L);

% 片側振幅スペクトルの計算
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

%% 3. 周波数軸の作成とプロット
% 周波数軸の作成
f = Fs*(0:(L/2))/L;

% 片側振幅スペクトルのプロット
figure;
plot(f, P1);
title(‘S(t)の片側振幅スペクトル’);
xlabel(‘周波数 (f)’);
ylabel(‘振幅 |P1(f)|’);
grid on;
xlim([0 200]); % 見やすいように表示範囲を制限
“`

4. より実践的なトピック (応用・注意点編)

基本が理解できたところで、次はより現実世界のデータに近い状況を想定し、周波数分析を行う上での重要な注意点について学んでいきましょう。

4-1. ノイズが含まれる信号の分析

実際のセンサーデータなどには、必ずと言っていいほどノイズが含まれています。周波数分析は、ノイズに埋もれた信号の特徴を抽出するのに非常に有効です。

先ほど作成した信号 S に、正規分布に従うランダムノイズ（ホワイトノイズ）を加えてみましょう。

“`matlab
% — ノイズの追加 —
noise_power = 0.5;
X = S + noise_power * randn(size(t)); % Sにノイズを加える

% ノイズが加わった波形をプロット
figure;
plot(t(1:100), X(1:100));
title(‘ノイズが付加された信号’);
xlabel(‘時間 (s)’);
ylabel(‘振幅’);
grid on;
“`


（生成されるグラフのイメージ：図2の波形がギザギザになったような波形）

図4: ノイズが付加された信号の時間領域波形

時間領域の波形はギザギザになり、元の滑らかな波形はほとんど見えなくなりました。
では、このノイズありの信号 X を周波数分析してみましょう。手順は先ほどと全く同じです。

“`matlab
% — ノイズあり信号のFFT —
Y_noise = fft(X);
P2_noise = abs(Y_noise/L);
P1_noise = P2_noise(1:L/2+1);
P1_noise(2:end-1) = 2*P1_noise(2:end-1);

% プロット
figure;
plot(f, P1_noise);
title(‘ノイズあり信号の片側振幅スペクトル’);
xlabel(‘周波数 (f)’);
ylabel(‘振幅 |P1(f)|’);
grid on;
xlim([0 200]);
“`


（生成されるグラフのイメージ：図3のグラフ全体に細かいギザギザのノイズフロアが乗り、その上に50Hzと120Hzのピークが突き出ている）

図5: ノイズあり信号の片側振幅スペクトル

結果を見ると、スペクトル全体にノイズ成分が乗っている（ノイズフロアが上昇している）ものの、50Hzと120Hzのピークは依然として明確に識別できます。
このように、時間領域ではノイズに埋もれて見えなかった信号成分を、周波数領域ではっきりと分離して検出できるのが、周波数分析の大きな利点の一つです。

4-2. 知っておくべき重要概念と注意点

FFTは非常に強力なツールですが、使い方を誤ると全く見当違いな結果を導いてしまいます。ここでは、絶対に知っておくべき2つの重要な注意点、「エイリアシング」と「スペクトル漏れ」について解説します。

サンプリング定理とエイリアシング（折り返し誤差）

「サンプリング周波数 Fs は、分析したい信号に含まれる最大周波数の2倍以上に設定する必要がある」というルールがあります。これをサンプリング定理と呼び、Fs/2 の周波数をナイキスト周波数と呼びます。

もし、このルールを破ってしまうと何が起こるでしょうか？
例えば、120Hzの信号を、サンプリング定理を満たさない Fs = 100Hz でサンプリングしてみます。ナイキスト周波数は 50Hz です。

“`matlab
Fs_low = 100; % 低いサンプリング周波数
T_low = 1/Fs_low;
L_low = 150;
t_low = (0:L_low-1)*T_low;

% 120Hzの信号
S_high_freq = sin(2pi120*t_low);

% FFT実行
Y_alias = fft(S_high_freq);
P2_alias = abs(Y_alias/L_low);
P1_alias = P2_alias(1:L_low/2+1);
P1_alias(2:end-1) = 2*P1_alias(2:end-1);

f_alias = Fs_low*(0:(L_low/2))/L_low;

% プロット
figure;
plot(f_alias, P1_alias);
title(‘エイリアシングの例 (120Hz信号をFs=100Hzでサンプリング)’);
xlabel(‘周波数 (f)’);
ylabel(‘振幅 |P1(f)|’);
grid on;
“`
結果のグラフを見ると、驚くべきことに、本来存在するはずのない 20Hz にピークが現れます。これがエイリアシング（折り返し誤差）です。


（生成されるグラフのイメージ：横軸0-50Hzの範囲で、20Hzの位置にピークが立つ）

図6: エイリアシング現象

なぜこうなるのでしょうか？
ナイキスト周波数(50Hz)を超える120Hzの信号は、サンプリング点が粗すぎるため、まるで低い周波数の信号であるかのように誤認識されてしまいます。計算上は |120 - 100| = 20Hz のように、ナイキスト周波数で折り返した位置に偽のスペクトル（エイリアス）が現れるのです。

対策:
1. 十分に高いサンプリング周波数を設定する: 分析したい信号の周波数成分を予測し、その最大周波数の2倍以上の Fs を選びます。安全マージンを見て、5倍～10倍程度に設定することも多いです。
2. アンチエイリアシングフィルタを使用する: サンプリングする前に、ナイキスト周波数以上の周波数成分を物理的なフィルタ（ローパスフィルタ）で予め除去しておきます。これにより、折り返してくる不要な高周波成分が存在しなくなるため、エイリアシングを根本的に防ぐことができます。

スペクトル漏れと窓関数

FFTは、切り出した有限長のデータが、無限に繰り返される周期信号であると仮定して計算を行います。

もし、切り出したデータの始まりと終わりの値が滑らかに繋がらない（つまり、切り出した区間の長さが信号の周期の整数倍でない）場合、その不連続な繋ぎ目が、本来存在しないはずの様々な周波数成分を生み出してしまいます。これにより、本来シャープであるべきスペクトルのピークが広がり、裾野が持ち上がってしまいます。この現象をスペクトル漏れ (Spectral Leakage) と呼びます。


（概念図：切り出した波形の始点と終点が一致せず、不連続になっている様子）

図7: スペクトル漏れの原因となる不連続点

このスペクトル漏れを軽減するために用いられるのが窓関数 (Window Function) です。
窓関数は、切り出したデータの両端を滑らかに0に近づけるような、山なりの形状をした関数です。元のデータにこの窓関数を掛け合わせることで、データの始まりと終わりの不連続性を緩和し、スペクトル漏れを抑えることができます。

MATLABで窓関数を適用するのは簡単です。代表的なハニング窓 (Hanning Window) を使ってみましょう。

“`matlab
% 周波数がFFTの周波数分解能の整数倍にならない信号を作成 (漏れを意図的に発生)
f_leak = 120.5; % 120Hzから少しずらす
S_leak = sin(2pif_leak*t);

% 窓関数なしでFFT
Y_leak = fft(S_leak);
P2_leak = abs(Y_leak/L);
P1_leak = P2_leak(1:L/2+1);
P1_leak(2:end-1) = 2*P1_leak(2:end-1);

% ハニング窓を作成
win = hanning(L)’; % L点のハニング窓、転置して行ベクトルにする

% 信号に窓関数を適用
S_win = S_leak .* win;

% 窓関数ありでFFT
Y_win = fft(S_win);
P2_win = abs(Y_win/L);
P1_win = P2_win(1:L/2+1);
P1_win(2:end-1) = 2*P1_win(2:end-1);

% 比較プロット
figure;
plot(f, P1_leak, ‘b’, ‘LineWidth’, 1.5);
hold on;
plot(f, P1_win, ‘r’, ‘LineWidth’, 1.5);
hold off;
title(‘窓関数の効果’);
xlabel(‘周波数 (f)’);
ylabel(‘振幅’);
legend(‘窓関数なし (矩形窓)’, ‘ハニング窓あり’);
grid on;
xlim([100 140]); % ピーク周辺を拡大
“`


（生成されるグラフのイメージ：青線（窓なし）はピークが太く、裾野が広く広がっている。赤線（窓あり）はピークは少し太くなるものの、裾野が劇的に低く抑えられている）

図8: 窓関数の適用によるスペクトル漏れの軽減効果

このグラフから、窓関数を適用しない場合（青線）はピークの裾野が広く、他の周波数帯にエネルギーが漏れ出しているのが分かります。一方、ハニング窓を適用した場合（赤線）は、裾野が劇的に低く抑えられ、スペクトル漏れが大幅に改善されていることが確認できます。

注意点: 窓関数を適用すると、ピークの振幅が少し低下したり、ピークが少し太くなる（周波数分解能がわずかに低下する）という副作用もあります。ハニング窓の他に、ハミング窓、ブラックマン窓など様々な種類があり、それぞれ特性が異なります。目的に応じて適切な窓関数を選択することが重要です。

4-3. パワースペクトルとパワースペクトル密度 (PSD)

これまで見てきたのは「振幅スペクトル」でしたが、信号のエネルギーやパワーに注目したい場合はパワースペクトルやパワースペクトル密度 (PSD) が用いられます。

  • パワースペクトル: 振幅の2乗に比例します。P_power = P1.^2 のように簡単に計算できます。信号のどの周波数にパワーが集中しているかを見るのに使います。
  • パワースペクトル密度 (PSD): パワースペクトルを周波数分解能で割ったもので、単位周波数あたりのパワーを表します。信号長やサンプリング周波数が異なる信号同士のパワーを比較する際に便利です。

MATLABには、より高度で安定したPSD推定を行うための関数 pwelch が用意されています。pwelch は、信号を短いセグメントに分割し、それぞれで窓関数を適用してFFTを行い、その結果を平均化する（ウェルチ法）ことで、ノイズの影響を抑えた滑らかなスペクトルを得ることができます。

matlab
% pwelch関数によるパワースペクトル密度推定
figure;
pwelch(X, hanning(500), 250, 512, Fs);
% pwelch(信号, 窓, オーバーラップ, FFT点数, サンプリング周波数)

fft を手動で計算するのに比べて、pwelch を使うとより手軽に高品質なスペクトルが得られる場合が多いので、特にノイズが多い信号を分析する際には試してみる価値があります。

5. まとめ

この記事では、MATLABを使ったフーリエ変換による周波数分析の基本を、理論的な背景から具体的な実装、そして実践的な注意点まで、詳細に解説してきました。

最後に、学んだことの要点を振り返りましょう。

  • フーリエ変換の核心: 「複雑な波は単純なサイン波の集まりである」という考え方に基づき、時間領域の信号を周波数ごとの成分に分解する強力なツールです。
  • MATLABでの実装:
    • fft 関数で簡単に高速フーリエ変換を実行できます。
    • FFTの出力は複素数であり、振幅と位相の情報を含んでいます。
    • abs(Y/L) で正規化された両側振幅スペクトルを計算し、そこから片側振幅スペクトルを求める手順を学びました。
    • f = Fs*(0:(L/2))/L のように、結果に対応する正しい周波数軸を作成することが不可欠です。
  • 実践上の重要事項:
    • エイリアシング: サンプリング周波数が不十分だと、高周波成分が低周波成分として誤認識される現象。サンプリング定理 (Fs > 2*f_max) を遵守することが対策の基本です。
    • スペクトル漏れ: 有限長のデータを切り出すことで生じる不連続性が原因で、スペクトルがぼやける現象。窓関数（ハニング窓など）を適用することで効果的に軽減できます。
  • 発展的なトピック:
    • 振幅の2乗であるパワースペクトルや、pwelch関数で推定できるパワースペクトル密度(PSD) も、信号のエネルギーを分析する上で重要です。

周波数分析は、信号処理の世界への入り口です。今回学んだFFTの知識をベースに、今後は時間と共に周波数が変化する信号を分析するための短時間フーリエ変換 (STFT) や、時間と周波数の両方で高い分解能を持つウェーブレット変換など、より高度な手法へと学びを進めていくことができるでしょう。MATLABには、これらの高度な分析をサポートする Signal Processing Toolbox も用意されています。

一見すると複雑な波形データも、フーリエ変換というメガネをかければ、その構造がはっきりと見えてきます。この記事が、あなたのデータ解析の旅における確かな第一歩となることを心から願っています。

コメントする

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です

上部へスクロール