instrumentation 2

%clear all;
close all;
clc;
warning off MATLAB:colon:nonIntegerIndex;
% ----- Question 1.1: Create a moving average filter with n points. -----
load('ECG.mat');
% Extract time-domain signal
t = ECG(:, 1);
ecg = ECG(:, 2);
n = length(ecg);
fs = 100;
% Calculate 5-point moving average
ecg_avg5 = moving_average(ecg, 5);
t_avg5 = t(3:(n-2));
% Calculate 15-point moving average
ecg_avg15 = moving_average(ecg, 15);
t_avg15 = t(8:(n-7));
% ----- Plot Question 1.1 -----
figure('NumberTitle', 'off', 'Name', 'Question 1.1 (20%)');
title('Moving Averages of ECG Signal');
xlabel('Time (s)');
ylabel('Amplitude (mV)');
hold on;
  l(1) = plot(t, ecg, 'b');
  l(2) = plot(t_avg5, ecg_avg5, 'r');
  l(3) = plot(t_avg15, ecg_avg15, 'g');
hold off;
legend('No Averaging', '5-Point Averaging', '15-Point Averaging');
%saveas(gcf, 'c:\temp\test.pdf', 'pdf');
% ----- Question 1.2: Calculate time-constant of filters. -----
% Model 5-point filter as IIR
B5pt = ones(1, 5) * (1/5);
A5pt = 1;
[H5pt, W5pt] = freqz(B5pt, A5pt, n, fs);
% Convert to dB magnitude
H5pt_dB = 10*log(abs(H5pt));
% This finds the point closest to -3dB... a reasonable approximation for
% us.
[min5pt, arrpos5pt] = min(abs(H5pt_dB - (-3)));
minus3dB_5pt = W5pt(arrpos5pt);
% Finally calculate the time constant
Tau5pt = 1/(2*pi*minus3dB_5pt);
% Repeat for 15-point filter. Model as IIR
B15pt = ones(1, 15) * (1/15);
A15pt = 1;
[H15pt, W15pt] = freqz(B15pt, A15pt, n, fs);
% Convert to dB magnitude
H15pt_dB = 10*log(abs(H15pt));
% And find the point closest to -3dB... reasonable enough approximation for
% us.
[min15pt, arrpos15pt] = min(abs(H15pt_dB - (-3)));
minus3dB_15pt = W15pt(arrpos15pt);
% Finally calculate the time constant
Tau15pt = 1/(2*pi*minus3dB_15pt);
% ----- Plot/Output Question 1.2 -----
figure('NumberTitle', 'off', 'Name', 'Question 1.2 (10%)');
subplot(2, 1, 1);
title('5-Point Moving Average Filter Impulse Response');
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
hold on;
  plot(W5pt, H5pt_dB, 'b');
  plot([0 minus3dB_5pt], [H5pt_dB(arrpos5pt) H5pt_dB(arrpos5pt)], 'r:');
  plot([minus3dB_5pt minus3dB_5pt], [H5pt_dB(arrpos5pt) -60], 'r:');
  text(9, -36, ['-3dB = ' num2str(minus3dB_5pt) ' Hz']);
  text(9, -44, ['tau = ' num2str(Tau5pt, 3) ' s']);
hold off;
subplot(2, 1, 2);
title('15-Point Moving Average Filter Impulse Response');
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
hold on;
  plot(W15pt, H15pt_dB, 'b');
  plot([0 minus3dB_15pt], [H15pt_dB(arrpos15pt) H15pt_dB(arrpos15pt)], 'r:');
  plot([minus3dB_15pt minus3dB_15pt], [H15pt_dB(arrpos15pt) -80], 'r:');
  text(5, -56, ['-3dB = ' num2str(minus3dB_15pt) ' Hz']);
  text(5, -65, ['tau = ' num2str(Tau15pt, 3) ' s']);
hold off;
% ----- Question 1.3: Generate a gradient plot. -----
N = 125; % Number of samples
gradientplot = zeros(n-N-1, 1);
% Loop through the full signal the number of times allowed by a 125 sample
% window
for i = 1:(n - N - 1)
    % Create a vector that is shows the difference between a sample and the
    % subsequent samples.
    diff_vec = ecg(i+1:(i+N+1)) - ecg(i:(i+N));
   
    % Since a logical true is also a numerical one, we can count the number
    % of elements above 0 and below 0 as:
    increase_count = sum(diff_vec > 0);
    decrease_count = sum(diff_vec < 0);
   
    gradientplot(i) = increase_count/decrease_count;
end
% ----- Plot/Output Question 1.3 -----
figure('NumberTitle', 'off', 'Name', 'Question 1.3 (20%)');
hold on;
plot(gradientplot);
plot([0 125], [1 1], 'k:');
hold off;
title('Gradient Plot of Original ECG Signal');
xlabel('Window Number');
ylabel('Ratio (Increases/Decreases)');
axis([0 125 0 1.5]);
% ----- Question 2.1: Determine if the transducer is suitable. -----
load 'Step';
% Split into two vectors
t_step = Step(:,1);
data_step = Step(:,2);
n_step = length(data_step);
fs_step = 1000;
% Subtract off the baseline 0.05V
data_step = data_step - data_step(1);
% Find the point where the value of the step function is 63% of
% steady-state
[step_63_min, step_63_index] = min(abs(data_step - (data_step(n_step)*0.63)));
step_tau_rough = t_step(step_63_index);
% Calculate a rough bandwidth
step_fc_rough = 1/(2*pi*(step_tau_rough/1000));
% ----- Plot/Output Question 2.1 -----
figure('NumberTitle', 'off', 'Name', 'Question 2.1 (20%)');
title('Pressure Sensor Step Response');
xlabel('Amplitude (V)');
ylabel('Time (ms)');
hold on;
  plot(t_step, data_step, 'b');
  plot([0 step_tau_rough], [data_step(step_63_index) data_step(step_63_index)], 'r:');
  plot([step_tau_rough step_tau_rough], [data_step(step_63_index) 0], 'r:');
  text(40, 0.3, ['tau = ' num2str(step_tau_rough) ' ms']);
  text(40, 0.25, ['fc = ' num2str(step_fc_rough, 3) ' Hz']);
  text(40, 0.2, 'Not suitable for 100 Hz sensor');
hold off;
% ----- Question 2.2: Determine bandwidth of sensor. -----
% Calculate approximate derivative of the step response
data_step_diff = gradient(data_step);
% Now take fft of the differentiated step response
data_step_diff_fft = abs(fft(data_step_diff));
data_step_diff_fft = data_step_diff_fft(1:n_step/2);
f_step = linspace(0, fs_step/2, n_step/2);
% Use interp1(...) to linearly interpolate the cutoff frequency. 56
% frequency 'buckets' over 500 Hz means that we only have a frequency
% resolution of roughly 9Hz, and we expect to lie somewhere in between the
% first two 'buckets'. Linear interpolation should give us a reasonable
% approximation. Note the reversal of 'x' and 'y' data in calling the
% method, as the function is designed to find a certain YI value given an
% XI, and we need to find XI (the frequency) from a given YI (0.707).
step_fc_precise = interp1(data_step_diff_fft, f_step, sqrt(2)/2);
% ----- Plot/Output Question 2.2 -----
figure('NumberTitle', 'off', 'Name', 'Question 2.2 (20%)');
title('Pressure Sensor Impulse Response Frequency Spectrum');
xlabel('Amplitude');
ylabel('Frequency (Hz)');
hold on;
  plot(f_step, data_step_diff_fft, 'b');
  plot([0 step_fc_precise], [sqrt(2)/2 sqrt(2)/2], 'r:');
  plot([step_fc_precise step_fc_precise], [sqrt(2)/2 0], 'r:');
  text(50, 0.3, ['fc = ' num2str(step_fc_precise) ' Hz']);
  text(50, 0.25, 'Not suitable for 100 Hz sensor');
hold off;
% ----- Question 2.3: Determine the sensor's static sensitivity. -----
% Determine Kss
Kss_step = data_step(n_step) / 1.0; % Newtons
% Display it!
disp('Question 2.3 (10%)');
disp('------------------');
disp(' ');
disp(['Kss = ' num2str(Kss_step) ' V/N with offset compensation']);