instrumentation 1

clear all;
close all;
clc;
warning off MATLAB:colon:nonIntegerIndex;
% ----- Question 1: Oversample signal and perform DFT. -----
load('ECG.mat');
% Extract time-domain signal
t_over = ECG(:, 1);
ecg_t_over = ECG(:, 2);
% Define sampling characteristics
fs_over = 100; % Oversampling frequency of 100Hz
n_over = length(ecg_t_over); % Number of samples
tmax = t_over(n_over); % Should be 2.49 seconds
% Calculate dft
ecg_f_over = fft(ecg_t_over);
% ----- Plot Question 1 -----
figure('NumberTitle', 'off', 'Name', 'Question 1 (20%)');
subplot(211);
  plot(t_over, ecg_t_over);
  title('Oversampled ECG Signal in Time Domain');
  xlabel('Time (s)');
  ylabel('Amplitude (mV)');
subplot(212);
  stem(abs(ecg_f_over));
  title('Raw DFT of Oversampled ECG Signal');
  xlabel('Sample #');
  ylabel('Amplitude (mV)');
% ----- Question 2: Scale both the frequency and magnitude axes. -----
% Calculate frequency axis
f_over = linspace(0, fs_over, n_over);
% Apply Parseval's theorem. Frequency spectrum needs to be divided by
% the square-root of the number of samples. This is because the POWER of
% the returned fft is a factor of N higher than it should be (fft is not
% unitary). Therefore the amplitude is incorrect by a factor of sqrt(n).
ecg_f_over_scaled = ecg_f_over / sqrt(n_over);
% ----- Plot Question 2 -----
figure('NumberTitle', 'off', 'Name', 'Question 2 (20%)');
% Plot only half of the spectrum as the other half is just a mirror image
% of the first half.
stem(f_over(1:n_over/2), abs(ecg_f_over_scaled(1:n_over/2)));
title('Properly Scaled DFT of Oversampled ECG Signal');
xlabel('Frequency (Hz)');
ylabel('Amplitude (mV)');
% ----- Question 3: Determine fm and use Nyquist theorem to find fs. -----
% Determining fm is done by hand and requires justification. No
% justification resulted in a loss of 5%. fm should be somewhere between
% 15Hz and 25Hz. This also needs to be plotted somewhere, failing to plot
% it resulted in a further loss of 5%.
fm = 22.5;
fs_min = fm * 2;
% ----- Plot Question 3 -----
figure('NumberTitle', 'off', 'Name', 'Question 3 (20%)');
hold on;
  stem(f_over(1:n_over/2), abs(ecg_f_over_scaled(1:n_over/2)));
  plot([fm fm], [0 0.8], 'r');
hold off;
title('Properly Scaled DFT of Oversampled ECG Signal');
xlabel('Frequency (Hz)');
ylabel('Amplitude (mV)');
% ----- Question 4: Resample and re-DFT signal with fs_min -----
% Resample the signal (done here using linear interpolation).
t_min = 0:1/fs_min:tmax;
ecg_t_min = interp1(t_over, ecg_t_over, t_min);
n_min = length(ecg_t_min);
% Calculate DFT
ecg_f_min = fft(ecg_t_min);
% Scale frequency and amplitude axes
f_min = linspace(0, fs_min, n_min);
ecg_f_min_scaled = ecg_f_min / sqrt(n_min);
% ----- Plot Question 4 -----
figure('NumberTitle', 'off', 'Name', 'Question 4 (20%)');
subplot(211);
  stem(f_over(1:n_over/2), abs(ecg_f_over_scaled(1:n_over/2)));
  title('DFT of Oversampled ECG Signal')
  xlabel('Frequency (Hz)');
  ylabel('Amplitude (mV)');
subplot(212);
  stem(f_min(1:n_min/2), abs(ecg_f_min_scaled(1:n_min/2)));
  title('DFT of Minimalistically sampled ECG Signal')
  xlabel('Frequency (Hz)');
  ylabel('Amplitude (mV)');
% ----- Question 5: Determine whether measurement sets are the same. -----
% Begin by loading and parsing the data
load('Measurements.mat');
meas1 = Measurements(:, 1);
meas2 = Measurements(:, 2);
n1=length(meas1);
n2=length(meas2);
% Now calculate mean and std. dev. for each set...
mean1 = mean(meas1);
mean2 = mean(meas2);
sd1 = std(meas1);
sd2 = std(meas2);
% ...and display them
disp(['Mean of set 1:      ' num2str(mean1)]);
disp(['Mean of set 2:      ' num2str(mean2)]);
disp(['Std. dev. of set 1: ' num2str(sd1)]);
disp(['Std. dev. of set 2: ' num2str(sd2)]);
disp(char(10));
% Calculate the Z-Score...
ZScore = abs(mean1-mean2)/sqrt(sd1^2/n1 + sd2^2/n2);
%...and display it too!
disp(['Z = ' num2str(ZScore)]);
disp(char(10));
% Now test the hypothesese
disp('First Test of Reproducibility');
disp('-----------------------------');
   
if (abs(ZScore) < 2.575)
    disp('Z < 2.575, thus Ha is confirmed (both sets of measurements are similar).');
    disp(char(10));
   
    disp('Second Test of Reproducibility');
    disp('------------------------------');
    if (abs(ZScore) < 0.0125)
        disp('Z < 0.0125, thus Ha is confirmed (both sets of measurements are the same).');
    else
        disp('Z > 0.0125, thus H0 is confirmed (both sets of measurements are not the same).');
    end
else
    disp('Z > 2.575, thus H0 is confirmed (both sets of measurements are different).');
end