snippets.m

% Matlab Snippets for SMB2

% Close all figures, clear workspace
close all; clear all; clc;

% Load ECG and set sampling frequency
    load e23;
    f0 = 250;
    signal = e23';
    t = 1/f0:1/f0:length(signal)/f0;
    
% Get number of samples
    N = size(e23, 1);

% Get impulse response from ecg
    h = flipud(x(1:21));
    
% Get length of impulse response
N = length(h);

% Add gaussian noise to signal
 x = 15*randn(length(x),1) + x;

% variables for Signals
    sR = 250;                                   % Sampling frequency
    T_A = 1 / sR;                               % Sampling interval
    N = 4000;                                   % Number of elements
    t = [0 : N-1] * T_A;                        % Time axis
    
% Breathing amplitude scaling
    scaling = 0.4;

% add ECG Signal synthesis: Noise
    signal = signal + scaling*randn(N,1)';      % Add gaussian noise

% Shift signal
    signal = signal';

% Plot ECG signal
    figure(1); clf; plot(t, signal);
    title('ECG Signal');
    xlabel('Time (sec)');
    ylabel('Amplitude');
    
% Compute spectrum
    F = fft(signal, N);                         % FFT
    F = abs( F ) / N;                           % Amplitude spectrum
    F = fftshift( F );                          % Shift FFT for convenience
    df = sR/N;                                  % FFT sampling interval
    frq = [0:N-1] * df;                         % Frequency axis
    frq = [-fliplr(frq(1:floor(end/2))),(frq(1:ceil(end/2)))];
    
% Plot spectrum
    figure(2); clf; stem(frq, F);
    title('Real ECG spectrum');
    xlabel('Frequency (Hz)');
    ylabel('Amplitude');
    xlim([floor(-sR/2),floor(sR/2)]);

% MATLAB derivative implementation / Ableitung 1 und 2
    y1_matlab = diff(signal);
    y2_matlab = diff(diff(signal));

    % For loop
    for i=3:length(signal)-2 %bordercontrol
    % do things here
  end
  
% Compute:  y = b*x + a
b = ( sum(x.*y) - 1/numel(x) * (sum(x) * sum(y)) ) / (sum(x.^2) - 1/numel(x)*(sum(x).^2) );
a = mean(y)-b*mean(x);

%Zerocrossings

diff_1 = [diff(e23)   ; 0];                         % First Derivation
diff_2 = [diff(diff_1) ; 0];                        % Second Derivation

local_maxima = zeros(1, N);                         % Matrix containing local maxima positions
local_minima = zeros(1, N);                         % Matrix containing local minima positions

% Find Zero Crossings
for k = 1:N-1
     if ((sign(diff_1(k)) * sign(diff_1(k+1)) < 0) | (sign(diff_1(k)) == 0))    % Detect Zero Crossing
          if (sign(diff_2(k)) < 0) 
               local_maxima(k+1) = 1;               % If second derivative is negative -> local maximum
          end
          if (sign(diff_2(k)) > 0)                  % If second derivative is positive -> local minimum
               local_minima(k+1) = 1; 
          end
     end
end

% Get Zero Crossings
idxMax = find(local_maxima);
idxMin = find(local_minima);

%% Create filter masks
mask_1 = zeros(1,3);                        % Flat structuring element
mask_2 = [zeros(1,2) 1 zeros(1,2)];         % peaky structuring element

%function and example dilation
function [ signal_output ] = dilate( f, k )
M = length(k);
N = length(f);
signal_output = zeros(1,length(f));

for m = 1 : 1 : N-M    
      values = [];    
   for n = 1:M;
        values = [values f(m+n)+k(n)];
   end    
    % Output signal adjusted to dilate_v2.m
    signal_output(m+round(M/2)) = max(values);
end

end

% Convolution Eq. (41)
y_sum = 0;
for n=(N+1):length(x)
   for i=1:N
       y_sum = y_sum + h(i) * x(n-i);
   end
   y(n-1) = y_sum;
   y_sum = 0;
end

% AMCD Eq. (42)
AMCD_sum = 0;
for n=(N+1):length(x)
   for i=1:N
       AMCD_sum = AMCD_sum + abs (  x(n-i) - h(i)  ) ;
   end
   y2(n-1) = AMCD_sum;
   AMCD_sum = 0;
end

%hilbert function
function a = hilbert_freq( x )

    % Reduce complex input to real only 
    if ~isreal(x)
        warning('Complex input is ignored.');
        x = real(x);
    end

    % Number of elements
    n = numel( x );
    % Compute FFT
    X = fft(x);
    % One-sided spectrum only
    XH = [X(1), 2*X(2:n/2), X(n/2+1), zeros(n/2-1,1)' ];
    % IFFT and return analytic signal
    a = ifft((XH)); 
    
end

% Normalize results
y = (y-min(y))/(max(y)-min(y));
y2 = (y2-min(y2))/(max(y2)-min(y2));
convolution = (conv(x,h)-min(conv(x,h)))/(max(conv(x,h))-min(conv(x,h)));

%interaction
    display('Reduce by hitting a key');
    pause