読者です 読者をやめる 読者になる 読者になる

粗大メモ置き場

個人用,たまーに来訪者を意識する雑記メモ

Matlab program for dft,fft (離散フーリエ変換)

メモ メモ-MATLAB フーリエ変換

突然ですがmaltabのreferenceページって結構見づらくないですか?

特にClassとして実装されてる奴らはどういう処理をするのかかなりわかりづらいと私は思っています。
どの操作が必要だったか忘れがちなので記しておきます。アホらしいですが...

fft in MATLAB

fftコマンドは非常に便利ですが,
ナイキスト周波数で折り返しが生じること
・ 周波数のスケーリングを自分で行う必要があること
の2点が実用上面倒です。

DFT用のコード(FFTも同じ)

以下のように関数化しました。パワースペクトル密度の表示まで意外とやることがあるのでこれで結構快適になったと思われます。

11/08 追記:直流成分の除去も中でするようにしました。単に平均を引いてやっています。

function [f power] = DFTandShow(t,x)
if size(t) == 1
    sampling = t;
else
    sampling = t(2) - t(1);
end

% transpose so that row become longer
if size(x,1) < size(x,2)
    x = x.';
end

sf = 1/sampling % Sample frequency (Hz)
L = max(size(x)) % data length(window length)

% calc average
average = ones(1,L)*x(:,1) / L;

x = x - average;

% fft
N = pow2(nextpow2(L));  % Transform length
N = 512;

y = fft(x,N); %fft

f = (0:N-1)*(sf/N);     % Frequency range ( sampling fq / fft length )
power =abs(y/N) *2 ;   % Power of the DFT

figure;
plot(f,power)
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('frequency (Hz)')
ylabel('Power')
grid on;
xlim([0 f(floor(N/2))]) % cut off faster frequency than nyquist freq
end