Main Content

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

フーリエ変換

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

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

yk+1=j=0n-1ωjkxj+1.

ω=e-2πi/nn 個ある 1 の複素根の 1 つで、i は虚数単位です。xy に対して、インデックス j とインデックス k の範囲は 0 から n-1 までです。

MATLAB® の関数 fft では、高速フーリエ変換アルゴリズムを使用してデータのフーリエ変換を計算します。時間 t の関数であり、周波数成分が 15 Hz と 20 Hz の正弦波信号 x について考えます。150 秒ごとに 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')

Figure contains an axes. The axes contains an object of type line.

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

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

信号の振幅を周波数の関数としてプロットすると、振幅のスパイクは、信号の周波数成分である 15 Hz および 20 Hz に対応します。

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

Figure contains an axes. The axes with title Magnitude contains an object of type line.

この変換では、スパイクのミラー コピーも生成され、これは信号の負の周波数に対応します。この周期性をさらにわかりやすく可視化するには、関数 fftshift を使用して、変換でゼロを中央にした循環シフトを実行します。

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

Figure contains an axes. The axes contains an object of type line.

ノイズを含む信号

科学的な応用では、ランダム ノイズによって信号が破損し、周波数成分が隠されることが多くあります。フーリエ変換によってランダム ノイズを排除し、周波数を明示することができます。たとえば、元の信号 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')

Figure contains an axes. The axes with title Power contains an object of type line.

計算の効率性

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

カリフォルニア沿岸沖で水中マイクから収集されたオーディオ データについて考えます。このデータは、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)])

Figure contains an axes. The axes contains an object of type line.

元の長さよりも大きい最小の 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')

Figure contains an axes. The axes contains an object of type line.

参考

| | | | | |

関連するトピック