近年來，心臟病已成為危害人類健康最常見的疾病。為了有效預防心臟疾病的發生，往往需要更加準確地采集與診斷心電信號，以便于更好地反映心臟情況。心電信號作為人體生理信號，對于識別心臟異常和心臟疾病具有重要的參考價值。心電信號，是由心臟的電活動引起心臟周圍的導電組織和體液反映到身體表面上而產生的電壓變化。正常心電信號主要由P波、QRS復合波、T波與U波組成。由于心電信號頻率與振幅較低，在采集心電信號時往往伴隨著不可避免的噪聲，因而降噪成為心電信號處理的核心。心電信號中存在多種低頻和高頻噪聲，常見的噪聲有加性高斯白噪聲(AdditiveWhiteGaussianNoise，AWGN)、基線漂移(BaselineWander，BW)、肌電干擾(MuscleArtefacts，MA)、運動偽影(ElectrodeMotion，EM)、工頻干擾(PowerLineInterference，PLI)。

基線漂移(BW):BW是一種低頻噪聲信號，主要是由呼吸、身體運動、電極接觸不良與皮膚電極阻抗所引起的。BW噪聲會使心電信號中的ST段與其他低頻成分失真。

工頻干擾(PLI):PLI主要是由供電設備產生的，它會扭曲低振幅心電信號波(P波或T波)的振幅與形狀。

肌電干擾(MA):MA作為最主要的心電噪聲信號之一，是由肌肉的收縮與顫動或身體突然運動而產生的電活動。肌電信號導致心電信號的局部波失真。

運動偽影(EM):由測量心電信號的電極與皮膚之間的阻抗隨著相對位移而發生改變所引起的，與P波、T波的頻譜幾乎完全重合。

ECG處理主要是從采集到有噪聲的信號中通過技術提取純信號，以便于更精確地分析與診斷心臟問題。因此，希望設計一種能夠完全消除噪聲而不影響信號的技術。然而事實上，這是不可能的。因為噪聲的頻率與信號的頻率在一個階段非常接近以至于重疊，會無法濾除噪聲或會使信號失真。因此，盡最大可能去除噪聲而不影響信號是目前最主要的技術。近年來，眾多學者采用不同學科不同領域方法來研究如何降噪，取得了較滿意效果。

鑒于此，提出一種基于集成經驗模態分解的心電信號降噪和基于希爾伯特變換的R峰檢測，運行環境為MATLAB 2018，主代碼如下：

%% Data Loading
ecg=load ('ecg1.mat');          % loading the signal 
ecg=struct2cell(ecg);
ecg=cell2mat(ecg);
ecg = (ecg - 1024)/200;     % you have to remove "base" and "gain"
ecg1=ecg(1,:);              %Noisy ecg
ecg2=ecg(2,:);              %filtered ecg
Fs =500;                    % sampling frequecy
t =linspace(0,length(ecg1)/Fs,length(ecg1)); %time vector%% ECG signal denoising
imf=eemd(ecg1,.2,70); %Apply the EEMD to the noisy signal .2->ratio of the standard deviation 70->ensemble number
imfs=imf';             %transpose the imf's matrix
reconstruction=imfs(4,:)+imfs(5,:)+imfs(6,:);  %We consider that these 3 imf's possess the important information%4 order Butterworth filter bandpass .05-230Hz. 
fclowpass=230; % Low pass cut-off frequency 230Hz
fchighpass=.05; % Low pass cut-off frequency .05Hz
filterorder=4;  %filter order[b,a]=butter(filterorder,[filterorder*fchighpass/Fs,2*fclowpass/Fs]);
filtered_ECG=filter(b,a,reconstruction);%% Emphasizing R peaks of the ECG
%Getting the maxima and minima of the ECG signal, to emphasize the R peaks
decg=(1/Fs)*(diff(filtered_ECG));  %derivative of the ecg
hecg=hilbert(decg); %hilbert transform of the derivative. 
envelope=abs(hecg);  %It returns the envelope of the ecg%% R peaks detection 
maximum=(max(envelope));
Threshold=.6*(maximum); 
[pks,locs] = findpeaks(envelope,'MinPeakHeight',Threshold);
time=(1/Fs)*length(ecg1);
timefactor=60/time;
cardiacFreq=round(timefactor*length(pks));
%% Plots
figure (1)
plot(t,ecg1); xlabel('time (s)'); ylabel('mV'); title('Raw ECG');
figure(2)
subplot(7,1,1);
plot(t,imfs(1,:)); xlabel('time (s)'); ylabel('mV'); title('Original ECG');
subplot(7,1,2);
plot(t,imfs(2,:)); xlabel('time (s)'); ylabel('mV'); title('1st IMF');
subplot(7,1,3);
plot(t,imfs(3,:)); xlabel('time (s)'); ylabel('mV'); title('2nd IMF');
subplot(7,1,4);
plot(t,imfs(4,:)); xlabel('time (s)'); ylabel('mV'); title('3rd IMF');
subplot(7,1,5);
plot(t,imfs(5,:)); xlabel('time (s)'); ylabel('mV'); title('4th IMF');
subplot(7,1,6);
plot(t,imfs(6,:)); xlabel('time (s)'); ylabel('mV'); title('5th IMF');
subplot(7,1,7);
plot(t,imfs(7,:)); xlabel('time (s)'); ylabel('mV'); title('6th IMF');
figure (3)
subplot(7,1,1);
plot(t,imfs(8,:)); xlabel('time (s)'); ylabel('mV'); title('7th IMF');
subplot(7,1,2);
plot(t,imfs(9,:)); xlabel('time (s)'); ylabel('mV'); title('8th IMF');
subplot(7,1,3);
plot(t,imfs(10,:)); xlabel('time (s)'); ylabel('mV'); title('9th IMF');
subplot(7,1,4);
plot(t,imfs(11,:)); xlabel('time (s)'); ylabel('mV'); title('10th IMF');
subplot(7,1,5);
plot(t,imfs(12,:)); xlabel('time (s)'); ylabel('mV'); title('11th IMF');
subplot(7,1,6);
plot(t,imfs(13,:)); xlabel('time (s)'); ylabel('mV'); title('12th IMF');
subplot(7,1,7);
plot(t,imfs(14,:)); xlabel('time (s)'); ylabel('mV'); title('13th IMF');
figure (4)
plot(t,reconstruction); xlabel('time (s)'); ylabel('mV'); title('IMFs reconstruction');
figure(5)
plot(t,filtered_ECG); xlabel('time (s)'); ylabel('mV'); title('Filtered ECG (IMF+bandPass)');
figure(6)
plot(t(1:9999),decg); xlabel('time (s)'); ylabel('d(ECG)/dt'); title('ECG derivative');
figure(7)
plot(t(1:9999),hecg); xlabel('time (s)'); ylabel('H(d(ECG)/dt)'); title('Hilbert Transform of derivate');
figure(8)
plot(t(1:9999),envelope); xlabel('time (s)'); ylabel('B(d(ECG)/dt)'); title('Envelope');
figure (9)
plot(t(1:9999),envelope,locs,pks,'or');
％完整代碼：mbd.pub/o/bread/mbd-ZJablZ5x
legend('ECG','R peaks','Location','NorthWest');

工學博士，擔任《Mechanical System and Signal Processing》《中國電機工程學報》《控制與決策》等期刊審稿專家，擅長領域：現代信號處理，機器學習，深度學習，數字孿生，時間序列分析，設備缺陷檢測、設備異常檢測、設備智能故障診斷與健康管理PHM等。