Exercise 5 (I+II)

part1.m

% EXERCISE SET 5 -- Part I

close all;
clear all;

%% Parameters

dt = 0.001; % s
dur = 1; % s
t = 0:dt:dur;

% Noise component
snr = 1;

% Frequency components of signal in format:
%  f1 weight1;
%  f2 weight2;
%  ...

freqs = ...
   [30 0.8;
    60 0.5;
    90 1;
    120 1;
    130 0.5;
    400 0.4;
    410 0.4;
    420 0.4];

%% Signal generation

signal = zeros(size(t));

for j = 1:size(freqs, 1)
    signal = signal + freqs(j,2) * sin(2 * pi * freqs(j,1) * t);
end

signal = signal + snr * randn(size(signal));

%% FT

[f, pfc] = powerspectrum(signal, dt);

%% Plot

figure;

subplot(2,1,1);
plot(t, signal);
xlabel('Time (s)');
ylabel('Amplitude (arb.)');

subplot(2,1,2);
plot(f, pfc / max(pfc));
xlabel('Frequency (Hz)');
ylabel('Contribution (relative to maximum)');

axis tight;

part2.m

% EXERCISE SET 5 -- Part II

close all;
clear all;

%% Parameters

dt = 0.001; % s
dur = 1.5; % s
t = 0:dt:dur;

% Noise component
snr = 2.5;

% Filter characteristics
cutoff = 12; % Hz

% Frequency components of signal in format:
%  f1 weight1;
%  f2 weight2;
%  ...

freqs = ...
   [10 1];

%% Signal generation

signal = zeros(size(t));

for j = 1:size(freqs, 1)
    signal = signal + freqs(j,2) * sin(2 * pi * freqs(j,1) * t);
end

nsignal = signal + snr * randn(size(signal));

%% Filter

fc = fft(nsignal);

nfc = fc;
cutoff_n = round(dur * cutoff);
nfc(cutoff_n:end-cutoff_n) = 0;

reconstructed = real(ifft(nfc));

%% Plot

figure;

a = subplot(3, 2, 1); plot(t, signal); xlabel('Time (s)'); ylabel('Amplitude (arb.)');
b = subplot(3, 2, 3); plot(t, nsignal); xlabel('Time (s)'); ylabel('Amplitude (arb.)');
c = subplot(3, 2, 5); plot(t, reconstructed); xlabel('Time (s)'); ylabel('Amplitude (arb.)');

linkaxes([b a c]);

% Power spectra
[f1, p1] = powerspectrum(signal, dt);
[f2, p2] = powerspectrum(nsignal, dt);
[f3, p3] = powerspectrum(reconstructed, dt);

d = subplot(3, 2, 2); plot(f1, p1); xlabel('Frequency (Hz)'); ylabel('Contribution');
e = subplot(3, 2, 4); plot(f2, p2); xlabel('Frequency (Hz)'); ylabel('Contribution');
f = subplot(3, 2, 6); plot(f3, p3); xlabel('Frequency (Hz)'); ylabel('Contribution');

linkaxes([d e f]);

powerspectrum.m

function [f, spec] = powerspectrum(signal, dt)
% Calculates power spectrum for supplied signal (including
% correct frequency axis in Hz) given sampling frequency 1/dt

fc = fft(signal);
pfc = fc .* conj(fc);
spec = pfc(1:round(length(pfc) / 2 + 1));

f = (1 / (dt*2)) * linspace(0, 1, round(length(fc) / 2 + 1));

end