MatlabでFFT

LSI

FFTで得られる周波数データにどのような性質があるのかを、簡単な波形を例にして調べてみました。Sampling Rateが1[Hz]で、周波数が 1/16, 2/16, 4/16, 8/16[Hz]の４通りで、コサインの波形をN=16点のFFTをしました。下のグラフで、上側にある４個の波形が時間軸でのコサインの波形、下の４個がそれぞれのFFTの結果です

  • 16個の周波数データの中の、non-zeroの周波数スペクトラムの位置
    • ０番目のデータはDCです。下のグラフでは常にゼロです。
    • 入力周波数が1/16 (= SamplingRate/N)のときには、non-zeroは１番目と１５番目に現れます。
    • 入力周波数が2/16 (= SamplingRate/N*2)のときには、non-zeroは２番目と１４番目に現れます。
    • 入力周波数が8/16 (= Nyquist周波数)のときには、non-zeroは８番目だけに現れます。
  • どの入力周波数でもパーセバルの定理が成立しています。
    • 入力周波数が1/16, 2/16, 4/16の場合は、時間軸の波形の２乗の総和が 8 (=N/2)となっていて、周波数スペクトラムの２乗の総和は (N/2)^2*2=N^2/2
    • 時間軸の波形の２乗の総和が入力周波数にかかわらずN/2となる理由ー時間軸の波形のある１点に着目して、そのデータと、それから９０度だけ位相がずれたデータのペアに着目すると、それらの２乗の和は必ず１になる（コサインとサインの関係のため）。そのようなペアがN/2個だけ存在するため。
    • 入力周波数が8/16の場合は、時間軸の波形の２乗の総和が 16 (=N)となっていて、周波数スペクトラムの２乗の総和は (N)^2=N^2

まとめ

  • 周波数スペクトラムの大きさの定義を、時間波形で振幅が１の場合に周波数スペクトラムの大きさが１と定義する場合には、FFTで得られた周波数データをN/2で割れば良い。
  • ただし、DCとNyquist周波数の場合はNで割る必要がある。ただし、実際の信号処理ではDCもNyquist周波数の成分も存在しないことが多いのであまり問題にはならない。(DCは無視するし、Nyquistはsystemのlowpass特性で減衰される)


clear;
N = 16;
fs = 1;
fx = [1/16, 2/16, 4/16, 8/16]';

x = cos(2*pi*fx/fs*[0:N-1]);
y = abs(fft(x)');

figure;
for n=1:4
   subplot(2,4,n);
   stem(1/fs*[0:N-1], x(n,:));
   grid;
   set(gca, 'XLIM', [0,N-1]);
   title(strcat('F_i_n = ', num2str(fx(n)),  ' [Hz]'));
   xlabel('Time [s]');

   subplot(2,4,n+4);
   stem(1/fs*[0:N-1], abs(fft(x(n,:))));
   grid;
   xlabel('Frequency');
   set(gca, 'XLIM', [0,N-1]);
   set(gca, 'YLIM', [0,N]);
end