% Physiobank data 
% normal data with Pan-Tompkins QRS detection
% and either peak detection or simple threshold crossing

clear all
figure(1); clf;

fname = 'e0103m.mat';
pname = 'Z:\Matlab\5030\PhysioBank\';

load([pname, fname])
Fs = 250 ; % Hz
t = linspace(0, length(val)/Fs, length(val));
noise = 100*sin(2*pi*60*t)+ 50*sin(2*pi*0.5*t) + 100*randn(1,length(t));
trunc = 400;
 val = val(1:fix(end/trunc))+ noise(1:fix(end/trunc))  ;
 t = t(1:fix(end/trunc));

np = 5;

subplot(np,1,1)
plot(t, val)
a = gca ;

%lowpass
 [b,a] = butter(2, 12/(Fs/2)) ;
 val_filtered = filter(b, a, val);
 subplot(np,1,2)
plot(t, val_filtered)

%highpass
[b,a] = butter(2, 5/(Fs/2), 'high') ;
val_filtered = filter(b, a, val_filtered);
subplot(np,1,3)
plot(t, val_filtered)

%differentiator --> squared
a = 1 ;
b = [0.25, 0.125, -0.125, -0.25];
val_filtered = filter(b, a, val_filtered);
val_filtered = val_filtered .^2 ;
subplot(np,1,4)
plot(t, val_filtered)

%smoothed
a = 1;
b = ones(1,20) / 20 ;
val_filtered = filter(b, a, val_filtered);

subplot(np,1,5)
plot(t, val_filtered);
hold on

% get the peaks
threshold = max(val_filtered)/3;
val_binary = diff(val_filtered > threshold) ;
start_index = find(val_binary == 1);
end_index =  find(val_binary == -1);
for i = 1:length(start_index)
    [local_max, max_index] = max(val_filtered(start_index:end_index)) ;
    local_max_time(i) = t(max_index + start_index(i)-1) ;
end

% get simple threshold crossing 
val_binary = diff(val_filtered > threshold) ;
start_index = find(val_binary == 1);
start_times = t(start_index) ;
plot(start_times, max(val_filtered)*0.5, 'r.','markersize', 20)
plot(local_max_time, max(val_filtered)*0.6, 'g.','markersize', 20)
subplot(np,1,1)
hold on
plot(start_times-0.1, 100+max(val_filtered)*0.5, 'r.','markersize', 20)
plot(local_max_time-0.1, 200+max(val_filtered)*0.5, 'g.','markersize', 20)

%!!!!!!!!!!!!!!!
return

figure(2)
subplot(2,1,1)
rr_interval_start = diff(start_times);
rr_interval_start = rr_interval_start(find(rr_interval_start < 1.5)) ;
hist(rr_interval_start, 100)
subplot(2,1,2)
rr_interval_peak = diff(local_max_time);
rr_interval_peak = rr_interval_peak(find(rr_interval_peak < 1.5)) ;
hist(rr_interval_peak, 100)
disp(['beats/min start=', num2str(60 / mean(rr_interval_start),'%2.2f')])
disp(['beats/min peak=', num2str(60 / mean(rr_interval_peak),'%2.2f')])

figure(3)
pwelch(val_filtered, [], [], [], Fs)




