Power spectrum demo

powerspec.matlab

%% Power spectrum demo

clear all;

dt = 0.01;

% Ten seconds at dt=0.01
t = 0:dt:10;

% Frequencies to be included in the signal. NB: Due
% to Shannon-Nyquist, 400Hz won't be "rendered" here.
freqs = [2 3 5 10 20 400];

% Generate the signal by weightless addition:
signal = zeros(size(t));
for idx = 1:length(freqs)
    signal = signal + sin(2 * pi * freqs(idx) * t);
end

% Add normally distributed noise with arbitrary weight
% of 3:
signal = signal + 3 * randn(size(t));

% Plot
figure;

plot(t, signal);

xlabel('Time (s)');
ylabel('Amplitude (a.u.)');
title('Signal');

%% Do DFT (Discrete Fourier Transform)

ft = fft(signal);

sf = 1 / dt; % Sampling frequency

% We need everything up to "the middle" (due to
% FT symmetry properties)
n = round(length(signal) / 2 + 1);
f_spectrum = (sf / 2) * linspace(0, 1, n);

% Calculate power spectrum
ps = (ft .* conj(ft)) / length(signal);

% Plot
figure;

plot(f_spectrum, ps(1:n));

xlabel('Frequency (Hz)');
ylabel('Energy of component (a.u.)');
title('Power spectrum');