% programa Lista3_derivada.m  S Furuie/24/05/2012
% Exemplo de detecao de QRS 
% 1) Ler sinal de ECG do physionet e caracterizar
% 5) derivadas
% 6) resultados tabela
%
% Utiliza:  plotATM modificado=> plotATM_sf

clear all           % clears all active variable
close all

%% lendo arquivo physionet  %%
%fnam = input('Nome do arquivo physionet (sem .mat, sem .info):','s');
fnam='104m';
fnam_mat=strcat(fnam,'.mat');
fnam_info=strcat(fnam,'.info');
[val faHz]=plotATM_sf(fnam_mat, fnam_info);% carrega matriz de dados em 'val' 
                    % com num.sinais=size(val,1)
                    % cada sinal com num.amostras=size(val,2) e freq. de
                    % amostragem faHz
ecg=val(1,:);  fs=faHz;
sze = length(ecg);
time = (1 : sze)/fs;
fN=fs/2;  % Nyquist
N_amostras=size(val,2);
tempoDisplay= N_amostras/fs;  %s
NumPtosDisplay=N_amostras; 
deltaT=1/fs;  % em s

% caracterizando os sinais de entrada
fid = fopen('resultados.txt','wt');
fprintf(fid,'*** ECG %s ***\n',fnam);
media=mean(ecg);  dp=std(ecg); maximo=max(ecg); minimo=min(ecg);
s=sprintf('ecg de entrada   : (media, dp, minimo, maximo)=(%10.2f, %10.2f %10.2f %10.2f)', media, dp, minimo, maximo);
disp (s); fprintf(fid, '%s \n',s);  


ecg_f_n=zeros(N_amostras,1);
ecg_f_mod=zeros(N_amostras,1);
% derivada backward
ecg_f_n(1)=(ecg(2)-ecg(1))/deltaT;
for i=2:N_amostras
    ecg_f_n(i)=(ecg(i)-ecg(i-1))/(deltaT);
    ecg_f_mod(i)=abs(ecg_f_n(i));
end
media=mean(ecg_f_n);  dp=std(ecg_f_n); maximo=max(ecg_f_n); minimo=min(ecg_f_n);
s=sprintf('derivada backward: (media, dp, minimo, maximo)=(%10.2f, %10.2f %10.2f %10.2f)', media, dp, minimo, maximo);
disp (s);fprintf(fid, '%s \n',s);
media=mean(ecg_f_mod);  dp=std(ecg_f_mod); maximo=max(ecg_f_mod); minimo=min(ecg_f_mod);
s=sprintf('|derivada backward|: (media, dp, minimo, maximo)=(%10.2f, %10.2f %10.2f %10.2f)', media, dp, minimo, maximo);
disp (s);fprintf(fid, '%s \n',s);

% display resultados
figure;
subplot(4,1,1);  plot(time(1:NumPtosDisplay),ecg(1:NumPtosDisplay));
ylabel('ECG(mV)'); xlabel('tempo (s)');
subplot(4,1,2); plot(time(1:NumPtosDisplay), ecg_f_n(1:NumPtosDisplay)); 
ylabel('derivada backward (mV/s)'); xlabel('tempo (s)');
subplot(4,1,3); plot(time(1:NumPtosDisplay), ecg_f_mod(1:NumPtosDisplay)); 
ylabel('|derivada backward| (mV/s)'); xlabel('tempo (s)');

% derivada central
ecg_f_n(1)=(ecg(2)-ecg(1))/deltaT;
ecg_f_n(N_amostras)=(ecg(N_amostras)-ecg(N_amostras-1))/deltaT;
for i=2:N_amostras-1
    ecg_f_n(i)=(ecg(i+1)-ecg(i-1))/(2*deltaT);
    ecg_f_mod(i)=abs(ecg_f_n(i));
end
media=mean(ecg_f_n);  dp=std(ecg_f_n); maximo=max(ecg_f_n); minimo=min(ecg_f_n);
s=sprintf('derivada central: (media, dp, minimo, maximo)=(%10.2f, %10.2f %10.2f %10.2f)', media, dp, minimo, maximo);
disp (s);fprintf(fid, '%s \n',s);
media=mean(ecg_f_mod);  dp=std(ecg_f_mod); maximo=max(ecg_f_mod); minimo=min(ecg_f_mod);
s=sprintf('|derivada central|: (media, dp, minimo, maximo)=(%10.2f, %10.2f %10.2f %10.2f)', media, dp, minimo, maximo);
disp (s);fprintf(fid, '%s \n',s);

% display resultados
subplot(4,1,4); plot(time(1:NumPtosDisplay), ecg_f_mod(1:NumPtosDisplay)); 
ylabel('|derivada central| (mV/s)'); xlabel('tempo (s)');
fclose(fid);
return;