%檔名稱  :  ECG_Plot
%實現功能  :  讀取MIT-BIH-DB檔案，讀取訊號，對訊號加噪聲，並利用濾波器去噪，並
%             實現影象輸出。
%參考資料  :  rddata.m    Author-Robert Tratnig 
%作者資訊  :  171848-張冰
%             [email protected]
%             
%修訂時間  :  2018年3月27日17點03分
%呼叫格式  :  無
%引數釋義  :  無

clc;clear all;

PATH= '';
HEADERFILE= '100.hea';
ATRFILE= '100.atr';
DATAFILE='100.dat';

%=========================對100.dat檔案資料的處理=========================%
%.hea檔案儲存了ECG的基本資訊
% 通過函式 fullfile 獲得標頭檔案的完整路徑
signalh= fullfile(PATH, HEADERFILE);
% 開啟標頭檔案，其識別符號為 fid1 ，屬性為'r'--“只讀”
fid1=fopen(signalh,'r');
% 讀取標頭檔案的第一行資料，字串格式
z= fgetl(fid1);
% 按照格式 '%*s %d %d %d' 轉換資料並存入矩陣 A 中
A= sscanf(z, '%*s %d %d %d',[1,3]);

nosig= A(1);    % 訊號通道數目
sfreq=A(2);     % 資料取樣頻率
SAMPLES2READ = 10*sfreq;    %取十秒資料

for k=1:nosig           % 讀取每個通道訊號的資料資訊
    z= fgetl(fid1);
    A= sscanf(z, '%*s %d %d %d %d %d',[1,5]);
    dformat(k)= A(1);           % 訊號格式; 這裡只允許為 212 格式
    gain(k)= A(2);              % 每 mV 包含的整數個數
    bitres(k)= A(3);            % 取樣精度（位解析度）
    zerovalue(k)= A(4);         % ECG 訊號零點相應的整數值
    firstvalue(k)= A(5);        % 訊號的第一個整數值 (用於偏差測試)
end;
fclose(fid1);
clear A;

%=========================對100.dat檔案資料的處理=========================%
%.dat檔案的資料格式讀取為每行三個位元組,即三個八位的二進位制數字,其內容含義為
%      0000 0000  ||             0000 0000              ||  0000 0000
%sign1(L)低八位資訊||左四位sign2(R)資訊,右四位sign1(L)資訊||sign2(R)低八位資訊
%將第二位元組的資訊處理後,後四位移至第一位元組最左位即得到完整的sign1
%將第二位元組的資訊處理後,前四位移至第一位元組最左位即得到完整的sign2.
signald = fullfile(PATH , DATAFILE);
fid2 = fopen(signald,'r');
A= fread(fid2, [3, SAMPLES2READ], 'uint8')';
fclose(fid2);

%對第二位元組做左位移運算，位移距離-4
%得到第二位元組左四位，即sign2的高四位，包括符號位,右高資訊
M_R_H = bitshift(A(:,2), -4);
%對第二位元組和1111做與運算，
%保留第二位元組右四位，即sign1的低四位，包括符號位,左高資訊
M_L_H = bitand(A(:,2), 15);
%對第二位元組和1000做與運算，
%保留第二位元組右邊第四位，獲取sign2符號位,並向左位移九位，與整體sign1進行運算
PRL=bitshift(bitand(A(:,2),8),9); 
%對第二位元組和10000000做與運算，
%保留第二位元組右邊第四位，獲取sign1符號位，並向左位移5位，與整體sign2進行運算
PRR=bitshift(bitand(A(:,2),128),5);

%M矩陣為sign1，2的儲存矩陣，儲存100.dat處理後資料
%將sign1(L)高位移至sign1低位前(A(:,1))
%將sign2(R)高位移至sign2低位前(A(:,3))
%最後將訊號符號位資訊去掉
M( : , 1)= bitshift(M_L_H,8)+ A(:,1)-PRL;
M( : , 2)= bitshift(M_R_H,8)+ A(:,3)-PRR;

%將sign的數值與零點做減法得到正負值
%再將得到的具有正負性的值與每mV的整數值相除，即得到電壓多少mV
M( : , 1)= (M( : , 1)- zerovalue(1))/gain(1);
M( : , 2)= (M( : , 2)- zerovalue(2))/gain(2);
%將我們設定的取樣個數除以頻率即得到這段樣品的測定時間。
TIME =(0:(SAMPLES2READ-1))/sfreq;
%釋放變數
clear A M_R_H M_L_H PRR PRL;

%===============================給訊號加白噪聲=============================%
%SNR 信噪比。訊號功率與噪聲功率的比值，不過一般取對數。即SNR=10*lg（A/B）。
%因為將倍數關係轉換為指數關係，所以設定分貝為單位。這裡加2%噪聲
SNR=10*log(100/2);
%利用matlab中的awgn函式給訊號加噪聲
Mwgn = awgn(M,SNR);

%==============================設計濾波器並濾波============================%
%三角濾波器
M_filter1 = [0,1,2,3,2,1,0];
%梯形濾波器
M_filter2 = [0,1,1,1,1,1,0];
%利用filter函式，將sign1進行與三角濾波器濾波並作平均，得到新的訊號值
Mflt(:,1) = filter(M_filter1,1,M(:,1) )/sum(M_filter1);
%利用conv2函式，將sign2進行與梯形濾波器濾波器做卷積，作平均，得到新的訊號值
%M_filter2需要旋轉180度，即做一次轉置
%效果與filter函式相同，這裡只是試一下不同的方法。
Mflt(:,2) = conv2(Mwgn(:,2) , M_filter2','same')/sum(M_filter2);

%=================================計算心率================================%
%觀察R值在0.5以上，並且在R值前後函式單調，所以在0.5範圍以上尋找極值確定R點位置
[R_V,R_L]=findpeaks(M(:,1),'minpeakheight',0.5);
%根據位置和取樣頻率來計算取樣區間內的平均心率
H_Rate = 60*(length(R_L)-1)/((R_L(length(R_L))-R_L(1))/sfreq);
%算出取樣的時間區間
R_Time = R_L(length(R_L))/sfreq;
%列印結果
disp('**********************診斷資訊**********************')
disp('**                                                **')
disp(['**','  該位測試物件在',num2str(R_Time,3),'秒內，的平均心率是：  '...
    num2str(H_Rate,3),'下  ','**']);
disp('**                                                **')
disp('**********************祝您健康**********************')

%===================================畫圖==================================%
%畫2個訊號在原圖/加噪聲（wgn）/濾波後（flt）的影象-----依次從上到下
%左列位訊號1，右列為訊號2
figure('NumberTitle', 'off', 'Name', 'ECG作圖');
clf, box on, hold on;
%因為影象較大較多，這裡讓畫圖面板最大化
set(gcf,'Position',get(0,'ScreenSize'))
%分配3*2的畫圖區間
s(1) = subplot(3,2,1);
s(2) = subplot(3,2,2);
s(3) = subplot(3,2,3);
s(4) = subplot(3,2,4);
s(5) = subplot(3,2,5);
s(6) = subplot(3,2,6);
string=['ECG signal ',DATAFILE];
%原圖
plot(s(1),TIME, M(:,1),'r');title(s(1),[string,' MLII']);
xlabel(s(1),'Time / s'); ylabel(s(1),'Voltage / mV');
plot(s(2),TIME, M(:,2),'b');title(s(2),[string,' V5']);
xlabel(s(2),'Time / s'); ylabel(s(2),'Voltage / mV');
%加噪聲圖（wgn）
plot(s(3),TIME, Mwgn(:,1),'r');title(s(3),[string,' MLII wgn']);
xlabel(s(3),'Time / s'); ylabel(s(3),'Voltage / mV');
plot(s(4),TIME, Mwgn(:,2),'b');title(s(4),[string,' V5 wgn']);
xlabel(s(4),'Time / s'); ylabel(s(4),'Voltage / mV');
%濾波後圖（flt）
%直接濾波原資料的圖
plot(s(5),TIME, Mflt(:,1),'r');title(s(5),[string,' MLII flt']);
xlabel(s(5),'Time / s'); ylabel(s(5),'Voltage / mV');
%濾波加噪聲後的圖
plot(s(6),TIME, Mflt(:,2),'b');title(s(6),[string,' V5 wgn flt']);
xlabel(s(6),'Time / s'); ylabel(s(6),'Voltage / mV');
%畫出兩個濾波器的形狀圖
figure('NumberTitle', 'off', 'Name', '濾波器');
clf, box on, hold on;
f(1) = subplot(2,1,1);
f(2) = subplot(2,1,2);
plot(f(1),M_filter1,'r');title(f(1),'三角濾波器');
plot(f(2),M_filter2,'b');title(f(2),'梯形濾波器');