FFT tutorial

Fs = 4096;             % set sampling frequency
T = 1/Fs;              % sample time
segdur = 4;            % segment duration is the duration of the measurement
L = Fs * segdur;       % number of data points in one segment
NFFT = 2^nextpow2(L);  % calculate nearest power of 2
t = (0:L-1)*T;         % array of times for each sample
% array of frequencies associated with one-sided FFT
f = Fs/2*linspace(0, 1, NFFT/2+1);

f0 = 1000;             % define a calibration line
a = 1;
cal = a*sin(2*pi*f0*t);

% notice we are see seeding randn and not rand
randn('seed', mod(sum(100*clock),2^32-1));

% simulate Gaussian white noise with calibration line
y1 = cal + 1*randn(size(t));
y2 = cal + 1*randn(size(t));

% perform FFTs -- note the normalization by sqrt(1/L)
y1_tilde = fft(y1, NFFT)*sqrt(1/L);
y2_tilde = fft(y2, NFFT)*sqrt(1/L);

% calculate one-sided power spectra -- note (2/L) normalization
CSD = real(conj(y2_tilde(1,1:NFFT/2+1)).*y1_tilde(1,1:NFFT/2+1)) * (2/L);
P1 = abs(y1_tilde(1:NFFT/2+1).^2) * (2/L);
P2 = abs(y2_tilde(1:NFFT/2+1).^2) * (2/L);

figure;
semilogy(f, abs(CSD), 'r');
print('-djpeg', 'csd.jpg');

Back to Resources/matlab