摘要：clear all;load idealECG.mat;load nECG.mat;fs = 128;N = length(idealECG);t = (0:N-1)/fs;D_t = 1/fs;% ECG modelling parameterstheta_
clear all;load idealECG.mat;load nECG.mat;fs = 128;N = length(idealECG);t = (0:N-1)/fs;D_t = 1/fs;% ECG modelling parameterstheta_i = [2*pi/3, 11*pi/12, pi, 13*pi/12, 3*pi/2];% theta_i = [-1*pi/3, -1*pi/12, 0, pi/12, pi/2];a_i = [1.2, -5, 30, -7.5, 0.75];b_i = [0.25, 0.1, 0.1, 0.1, 0.4];omega = 2*pi/0.8;% Synthetic ECG signalTheta = zeros(size(nECG));z = zeros(size(nECG));for n = (2:N)Theta(n) = mod(Theta(n-1) + omega*D_t,2*pi);z_dot = -1.5e-3*sin(2*pi*0.2*n/fs);for i = (1:5)D_theta = Theta(n)-theta_i(i);z_dot = z_dot - a_i(i)*D_theta*exp(-0.5*D_theta^2/b_i(i)^2);endz(n) = z(n-1) + D_t*z_dot;endfigure(1);plot(t,z);xlim([0,8]);
a_i = [1.2, -5, 30, -7.5, 0.75]*20;% EKFTheta = zeros(size(nECG));z = zeros(size(nECG));R = eye(2)*1e-6;Q = 0.027;sigma = eye(2)*1;C = [0 1];for n = (2:N)%%%%%%%%% Prediction Theta(n) = mod(Theta(n-1) + omega*D_t,2*pi);z_dot = 0;for i = (1:5)D_theta = Theta(n)-theta_i(i);z_dot = z_dot - a_i(i)*D_theta*exp(-0.5*D_theta^2/b_i(i)^2);endz(n) = z(n-1) + D_t*z_dot;s = [Theta(n), z(n)]';%%%%%%%%% Jacobian of state transition G = eye(2);for i = (1:5)D_theta = Theta(n)-theta_i(i);G(2,1) = G(2,1) - D_t * a_i(i)*(1 - D_theta^2 / b_i(i)^2) * exp(-0.5*D_theta^2 / b_i(i)^2);end%%%%%%%%% Error covariance of prediction sigma = G*sigma*G' + R;%%%%%%%%% Kalman Gain K = sigma*C'/(C*sigma*C' + Q);%%%%%%%%% Correction using input signal s = s + K*(nECG(n) -C*s);Theta(n) = s(1);z(n) = s(2);sigma = (eye(2) - K*C)*sigma;end% MSE = mean((z-idealECG).^2)figure(2);subplot(3, 1, 1);plot(t,nECG);title('Noisy ECG Signal');xlabel('Samples');ylabel('Amplitude');xlim([5,10]);subplot(3, 1, 2);plot(t,z);title('Filtered ECG Signal (EKF)');xlabel('Samples');ylabel('Amplitude');xlim([5 10]);subplot(3,1,3);plot(t,idealECG);xlim([5,10]);title('ideal ECG');
% EKFTheta = zeros(size(nECG));z = zeros(size(nECG));R = eye(2)*1;Q = 0.027;sigma = eye(2)*1;C = [0 1];for n = (2:N)%%%%%%%%% Prediction Theta(n) = mod(Theta(n-1) + omega*D_t,2*pi);z_dot = 0;for i = (1:5)D_theta = Theta(n)-theta_i(i);z_dot = z_dot - a_i(i)*D_theta*exp(-0.5*D_theta^2/b_i(i)^2);endz(n) = z(n-1) + D_t*z_dot;s = [Theta(n), z(n)]';%%%%%%%%% Jacobian of state transition G = eye(2);for i = (1:5)D_theta = Theta(n)-theta_i(i);G(2,1) = G(2,1) - D_t * a_i(i)*(1 - D_theta^2 / b_i(i)^2) * exp(-0.5*D_theta^2 / b_i(i)^2);end%%%%%%%%% Error covariance of prediction sigma = G*sigma*G' + R;%%%%%%%%% Kalman Gain K = sigma*C'/(C*sigma*C' + Q);%%%%%%%%% Correction using input signal s = s + K*(nECG(n) -C*s);Theta(n) = s(1);z(n) = s(2);sigma = (eye(2) - K*C)*sigma;end% MSE = mean((z-idealECG).^2)figure(2);subplot(3, 1, 1);plot(t,nECG);title('Noisy ECG Signal');xlabel('Samples');ylabel('Amplitude');xlim([5,10]);subplot(3, 1, 2);plot(t,z);title('Filtered ECG Signal (EKF)');xlabel('Samples');ylabel('Amplitude');xlim([5 10]);subplot(3,1,3);plot(t,idealECG);xlim([5,10]);title('ideal ECG');
% EKFTheta = zeros(size(nECG));z = zeros(size(nECG));R = eye(2)*1e-6;Q = 0.027;sigma = eye(2)*1;MSE = 1;C = [0 1];for j = linspace(0.01,0.07,60)for k = linspace(0.01,0.07,60)R = [k, 0; 0, j];for n = (2:N)%%%%%%%%% Prediction Theta(n) = mod(Theta(n-1) + omega*D_t,2*pi);z_dot = 0;for i = (1:5)D_theta = Theta(n)-theta_i(i);z_dot = z_dot - a_i(i)*D_theta*exp(-0.5*D_theta^2/b_i(i)^2);endz(n) = z(n-1) + D_t*z_dot;s = [Theta(n), z(n)]';%%%%%%%%% Jacobian of state transition G = eye(2);for i = (1:5)D_theta = Theta(n)-theta_i(i);G(2,1) = G(2,1) - D_t * a_i(i)*(1 - D_theta^2 / b_i(i)^2) * exp(-0.5*D_theta^2 / b_i(i)^2);end%%%%%%%%% Error covariance of prediction sigma = G*sigma*G' + R;%%%%%%%%% Kalman Gain K = sigma*C'/(C*sigma*C' + Q);%%%%%%%%% Correction using input signal s = s + K*(nECG(n) -C*s);Theta(n) = s(1);z(n) = s(2);sigma = (eye(2) - K*C)*sigma;endMSE_jk = mean((z-idealECG).^2);if (MSE> MSE_jk)MSE = MSE_jk;k_m = k;j_m = j;z_m = z;endendend% MSE = mean((z-idealECG).^2)figure(2);subplot(3, 1, 1);plot(t,nECG);title('Noisy ECG Signal');xlabel('Samples');ylabel('Amplitude');xlim([5,10]);subplot(3, 1, 2);plot(t,z_m);title('Filtered ECG Signal (EKF)');xlabel('Samples');ylabel('Amplitude');xlim([5 10]);subplot(3,1,3);plot(t,idealECG);xlim([5,10]);title('ideal ECG');
知乎学术咨询：https://www.zhihu.com/consult/people/792359672131756032?isMe=1担任《Mechanical System and Signal Processing》《中国电机工程学报》等期刊审稿专家，擅长领域：信号滤波/降噪，机器学习/深度学习，时间序列预分析/预测，设备故障诊断/缺陷检测/异常检测。分割线分割线分割线分割线分割线分割线分割线分割线分割线分割线一维时间序列信号的稀疏度度量方法（MATLAB R2018A）算法运行环境为MATLAB R2018A，执行一维信号的稀疏度量方法，包括峰度(Kurt)、负熵(NE)、d -范数(DN)、2-范数与1-范数之比(L2/L1)、基尼指数(GI)、修正平滑指数(MSI)、基尼指数2 (GI2)、基尼指数3 (GI3)、广义基尼指数(GGI)、完全广义基尼指数等。算法可迁移至金融时间序列，地震信号，机械振动信号，语音信号，声信号，生理信号（EEG,EMG）等一维时间序列信号。完整代码可通过知乎付费咨询获得：https://www.zhihu.com/consult/people/792359672131756032 来源：我可不太甜
免责声明：本站系转载，并不代表本网赞同其观点和对其真实性负责。如涉及作品内容、版权和其它问题，请在30日内与本站联系，我们将在第一时间删除内容!
