MATLAB ECG 数字信号处理 医学信号处理

Question

实验要求编写求解Yule-Walker方程的程序，并对实际生理信号（例如心电、脑电）建立AR模型。
对同一数据，使用Matlab信号处理工具箱自带函数aryule计算相同阶数AR模型系数，检验程序是否正确。
用伪随机序列（白噪声）驱动AR模型，观察输出是否与真实心电、脑电信号相似，对比真实信号与仿真信号的功率谱。
但是我多次导入了心电信号都没有出正确的图
我找的心电信号部分如下，是在https://archive.physionet.org/cgi-bin/atm/ATM里找的
想知道哪里错了

clear; clc;
load ecgdata;
x = ecgdata (1:1024); % 长度可以任意选择，但信号越长计算量越大
% load eegdata;
% x = eegdata (1:1024); % 长度可以任意选择，但信号越长计算量越大

p = 12; % 尝试改变模型阶数，观察效果
Rxx = xcorr(x,'biased');
Rtemp = zeros(1,p);
Rl = zeros(p,1);
for k = 1:length(Rtemp)
    Rtemp(k) = Rxx(length(x)-1+k);
    Rl(k) = Rxx(length(x)+k);
end
Rs = toeplitz(Rtemp); % 生成自相关系数矩阵(Toeplitz型)
A = -inv(Rs)*Rl; % AR模型系数估计
Sw = [Rtemp(1),Rl']*[1;A]; % 白噪声方差估计

% 采用malab自带函数估计模型系数
[a,E] = aryule(x,p); % a--系数，E--预测误差，k--反射系数
da = a(2:end)-A' % 自编程序求解是否正确？

w = randn(size(x));
x2 = filter(1,a,w); % 仿真数据
figure;
subplot(3,1,1);
plot(x);
title('真实数据');
subplot(3,1,2);
plot(x2);
title('仿真数据');

Answer 1

可以借鉴下

% Filtering
for i = 1:N,
    Pbar(:,:,i) = Pminus;
    Xbar(:,i) = Xminus;

    if(i<order+1)
        H = [-data(i-1:-1:1)' zeros(1,order-i+1)];
    else
        H = -data(i-1:-1:i-order)';
    end

    Yminus = H * Xminus + Vmean;
    inov = data(i)-Yminus;

    K = Pminus * H'/(H * Pminus * H' + alpha*R);
    Pplus = ( (eye(order) - K * H) * Pminus * (eye(order) - K * H)' + K * R * K' )/alpha;
    Xplus = Xminus + K*(inov);                                % A posteriori state estimate

    Xminus = A * Xplus + Wmean;                              % State update
    Pminus = A * Pplus * A' + Q;

    ARCoefs(:,i) = Xplus;
%     P(i) = max(diag(Pplus));
    Phat(:,:,i) = Pplus;
    Xhat(:,i) = Xplus;
end

% Smoothing
PSmoothed = zeros(size(Phat));
X = zeros(size(Xhat));
PSmoothed(:,:,N) = Phat(:,:,N);
X(:,N) = Xhat(:,N);
for k = N-1:-1:1,
    S = Phat(:,:,k) * A' /Pbar(:,:,k+1);
    X(:,k) = Xhat(:,k) + S * (X(:,k+1) - Xbar(:,k+1));
    PSmoothed(:,:,k) = Phat(:,:,k) - S * (Pbar(:,:,k+1) - PSmoothed(:,:,k+1)) * S';
end

ARCoefs2 = X;

 

noise0 =  NoiseGenerator('white',SignalPower,SNR,N,1);
noise1 =  NoiseGenerator('colored',SignalPower,SNR,N,fs,beta,10);
noise2 =  NoiseGenerator('MA',SignalPower,SNR,N,fs,0);
noise3 =  NoiseGenerator('EM',SignalPower,SNR,N,fs,100);
noise4 =  NoiseGenerator('BW',SignalPower,SNR,N,fs,100);
noise5 =  NoiseGenerator('mixture',SignalPower,SNR,N,fs,[w_bw,w_em,w_ma],1000);

x0 = data + noise0;
x1 = data + noise1;
x2 = data + noise2;
x3 = data + noise3;
x4 = data + noise4;
x5 = data + noise5;

% Plot results
figure;
hold on
plot(t,noise0,'g');
plot(t,x0,'r');
plot(t,data,'b');
grid;
xlabel('time (sec.)');
legend('Noise','Noisy ECG','Original ECG');
title('White Gaussian Noise');

figure;
hold on
plot(t,noise1,'g');
plot(t,x1,'r');
plot(t,data,'b');
grid;
xlabel('time (sec.)');
legend('Noise','Noisy ECG','Original ECG');
title('Colored Gaussian Noise');

figure;
hold on
plot(t,noise2,'g');
plot(t,x2,'r');
plot(t,data,'b');
grid;
xlabel('time (sec.)');
legend('Noise','Noisy ECG','Original ECG');
title('Real Muscle Artifact Noise');

figure;
hold on
plot(t,noise3,'g');
plot(t,x3,'r');
plot(t,data,'b');
grid;
xlabel('time (sec.)');
legend('Noise','Noisy ECG','Original ECG');
title('Real Electrode Movement Noise');

figure;
hold on
plot(t,noise4,'g');
plot(t,x4,'r');
plot(t,data,'b');
grid;
xlabel('time (sec.)');
legend('Noise','Noisy ECG','Original ECG');
title('Real Baseline Wander Noise');

figure;
hold on
plot(t,noise5,'g');
plot(t,x5,'r');
plot(t,data,'b');
grid;
xlabel('time (sec.)');
legend('Noise','Noisy ECG','Original ECG');
title('Mixture or BW, EM, and MA Noise');