フーリエ変換

R2026b

フーリエ変換は、時間または空間でサンプリングされた信号を、時間周波数または空間周波数でサンプリングされた同じ信号に変換する公式です。信号処理では、フーリエ変換により、信号の重要な特性、つまり周波数成分が明らかになる場合があります。

ベクトル $$x$$ が $$n$$ 個の等間隔でサンプリングされた点をもつ場合、フーリエ変換は次の式で定義されます。

$$y_{k+1}=\sum _{j=0}^{n-1}{\omega }^{jk}{x}_{j+1}.$$

$$\omega =e^{-2\pi i/n}$$ は $$n$$ 個ある 1 の複素根の 1 つで、$$i$$ は虚数単位です。$$x$$ と $$y$$ に対して、インデックス $$j$$ とインデックス $$k$$ の範囲は $$0$$ から $$n-1$$ までです。

MATLAB® の関数 fft では、高速フーリエ変換アルゴリズムを使用してデータのフーリエ変換を計算します。時間 t の関数であり、周波数成分が 15 Hz と 20 Hz の正弦波信号 x について考えます。1/50 秒ごとに 10 秒間にわたってサンプリングされた時間ベクトルを使用します。

Ts = 1/50;
t = 0:Ts:10-Ts;
x = sin(2*pi*15*t) + sin(2*pi*20*t);
plot(t,x)
xlabel('Time (seconds)')
ylabel('Amplitude')

信号のフーリエ変換を計算し、周波数空間での信号のサンプリングに対応するベクトル f を作成します。

y = fft(x);
fs = 1/Ts;
f = (0:length(y)-1)*fs/length(y);

変換後の信号の振幅を周波数の関数としてプロットします。

plot(f,abs(y))
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Magnitude')

このプロットには周波数ピークが 4 つ示されていますが、信号は 15 Hz と 20 Hz の 2 つの周波数ピークをもつと予想されます。このプロットの後半は、前半と左右対称になっています。時間領域信号の離散フーリエ変換は周期的な性質をもっており、スペクトルの前半が正の周波数、後半が負の周波数となります。このプロットの 30 Hz と 35 Hz の周波数成分は、–20 Hz と –15 Hz の周波数成分に対応します。この周期性をさらにわかりやすく可視化するには、関数 fftshift を使用して、変換でゼロを中央にした循環シフトを実行します。

n = length(x);
fshift = (-n/2:n/2-1)*(fs/n);
yshift = fftshift(y);
plot(fshift,abs(yshift))
xlabel('Frequency (Hz)')
ylabel('Magnitude')

ノイズを含む信号

科学的な応用では、ランダム ノイズによって信号が破損し、周波数成分が隠されることが多くあります。フーリエ変換によってランダム ノイズを排除し、周波数を明示することができます。たとえば、元の信号 x にガウス ノイズを挿入することにより、新しい信号 xnoise を作成します。

rng('default')
xnoise = x + 2.5*randn(size(t));

周波数の関数としての信号パワーは、信号処理で多く使用される測定基準です。パワーは、信号のフーリエ変換の振幅を二乗して、周波数サンプルの数で正規化したものです。中央を周波数ゼロに揃えて、ノイズを含む信号のパワー スペクトルを計算してプロットします。ノイズがあっても、パワーのスパイクにより、信号の周波数を判別できます。

ynoise = fft(xnoise);
ynoiseshift = fftshift(ynoise);
power = abs(ynoiseshift).^2/n; 
plot(fshift,power)
title('Power')
xlabel('Frequency (Hz)')
ylabel('Power')

計算の効率性

フーリエ変換式を直接使用して $$y$$ の $$n$$ 個の要素をそれぞれ計算するには、約 $$n^2$$ 回の浮動小数点演算が必要です。高速フーリエ変換アルゴリズムでは、計算に必要な演算は約 $n\log n$ 回のみです。この計算の効率性は、数百万ものデータ点をもつデータを処理する場合は大きなメリットになります。特化された高速フーリエ変換アルゴリズムを多数実装すると、$$n$$ が 2 のべき乗である場合など、$$n$$ に小さい素因数がある場合にさらに効率的になります。

カリフォルニア沿岸沖で水中マイクから収集されたオーディオ データについて考えます。このデータは、Cornell University Bioacoustics Research Program で保管されているライブラリにあります。bluewhale.au のデータのサブセットを読み込んでフォーマットします。このファイルには太平洋のシロナガスクジラの鳴き声が収録されています。シロナガスクジラの鳴き声は低い周波数の音であるため、人間はほとんど聞き取れません。データの時間スケールは、音程を上げて鳴き声がよりはっきりと聞こえるように 10 倍に圧縮されています。コマンド sound(x,fs) を使用すると、オーディオ ファイル全体を聴くことができます。

whaleFile = 'bluewhale.au';
[x,fs] = audioread(whaleFile);
whaleMoan = x(2.45e4:3.10e4);
t = 10*(0:1/fs:(length(whaleMoan)-1)/fs);

plot(t,whaleMoan)
xlabel('Time (seconds)')
ylabel('Amplitude')
xlim([0 t(end)])

元の長さよりも大きい最小の 2 のべき乗である新しい信号長を指定します。次に、fft を使用して、新しい信号長を使ってフーリエ変換を計算します。fft では、サンプル サイズを拡大するために、データにゼロが自動的に付加 (パディング) されます。このパディングにより、特にサンプル サイズに大きな素因数が含まれている場合に、変換の計算が大幅に速くなります。

m = length(whaleMoan);
n = pow2(nextpow2(m));
y = fft(whaleMoan,n);

信号のパワー スペクトルをプロットします。プロットは、うめき声が 17 Hz 前後の基本周波数と、2 番目の調波が強い高調波群から構成されていることを示しています。

f = (0:n-1)*(fs/n)/10; % frequency vector
power = abs(y).^2/n;   % power spectrum
plot(f(1:floor(n/2)),power(1:floor(n/2)))
xlabel('Frequency (Hz)')
ylabel('Power')

正弦波の位相

フーリエ変換を使用して、元の信号の位相スペクトルを抽出することもできます。たとえば、周波数が 15 Hz と 40 Hz の 2 つの正弦波で構成される信号を作成します。最初の正弦波は位相が $-\pi /4$ の余弦波で、2 番目の正弦波は位相が $\pi /2$ の余弦波です。信号を 100 Hz で 1 秒間サンプリングします。

fs = 100;
t = 0:1/fs:1-1/fs;
x = cos(2*pi*15*t - pi/4) - sin(2*pi*40*t);

信号のフーリエ変換を計算します。変換の振幅を周波数の関数としてプロットします。

y = fft(x);
z = fftshift(y);

ly = length(y);
f = (-ly/2:ly/2-1)/ly*fs;

stem(f,abs(z))
xlabel("Frequency (Hz)")
ylabel("|y|")
grid

振幅の小さい変換の値を除去して、変換の位相を計算します。位相を周波数の関数としてプロットします。

tol = 1e-6;
z(abs(z) < tol) = 0;

theta = angle(z);

stem(f,theta/pi)
xlabel("Frequency (Hz)")
ylabel("Phase / \pi")
grid

参考

fft | fftshift | nextpow2 | ifft | fft2 | fftn | fftw

トピック

• 2 次元フーリエ変換

授業用リソース

• Fourier Analysis (MathWorks Teaching Resources)