% 以下循环的目的是找出FIR滤波器合适的阶数
% 判据是当阶数增加而均方误差没有明显下降时，则认为阶数足够 while abs(e0-e1)>1e-6 % e1和e0不够接近则循环 N = N+1; e0 = e1;
Rxs = Rss(M:(M+N-1));
Rxx = Rww(M:(M+N-1))+Rss(M:(M+N-1)); R_xx = zeros(N); for j = 1:N for n = 1:N
R_xx(j,n) = Rxx(abs(j-n)+1); end end
h = inv(R_xx)*Rxs'; e1 = Rss(M)-h'*Rxs'; end
N % 显示N的最终值 e = e1;
% 主程序 clear; clc;
M = input('信号的长度 M = '); n = 1:M;
s = exp(-0.002*n).*sin(pi*n/50); % 仿真信号，可以自己生成，任意形式 % load ecgdata; % 实际心电信号 % s = x(1:M)';
w = 0.4*randn(1,M); % 白噪声，系数代表噪声相对强度 x = s+w; % 仿真信号
Rss = xcorr(s,s); %估计信号自相关函数
Rww = xcorr(w,w); % 估计噪声自相关函数 [h,e] = WH(Rss,Rww,M);
ss = filter(h,1,x); %用维纳滤波器滤波 figure; subplot(2,2,1); plot(n,s); title('信号'); subplot(2,2,2); plot(n,w); title('噪声'); subplot(2,2,3); plot(n,x); title('观测值'); subplot(2,2,4); plot(n,ss); title('信号估计'); figure; plot(n,ss-s); title('估计误差');
实验四 Yule-Walker方程
（一）实验目的
学习求解Yule-Walker方程，建立随机信号的AR模型。
（二）实验原理
随机信号可以看作是由当前激励白噪声w(n)以及若干次以往信号x(n-k)的线性组合产生，即所谓自回归模型（AR模型）
x?n??w?n???akx?n?k?k?1pH?z??11??akz?kk?1p
模型参数满足Yule-Walker方程
Rxx?m????akRxx?m?k?m?0
k?1p矩阵形式
R??1??R?0??R?1?R?0???????R?p?1?R?p?2??R?1?p???a1??R?1???a??R?2???R?2?p??2??????? ??????????????a???R?0???Rp?p????求解Yule-Walker方程，就可以得到AR模型系数
当模型阶次较大时，直接用矩阵运算求解的计算量大，不利于实时运算。利用系数矩阵的特性，人们提出了如L-D算法等快速算法。
（三）实验内容和步骤
编写求解Yule-Walker方程的程序，并对实际生理信号（例如脑电）建立AR模型。
对同一数据，使用matlab信号处理工具箱自带函数aryule计算相同阶数AR模型系数，检验程序是否正确。
用伪随机序列（白噪声）驱动AR模型，观察输出是否与真实脑电信号相似，对比真实信号与仿真信号的功率谱。
（四）思考题
对EEG建模后，该模型产生的信号是否能反映EEG信号的特征？
（五）实验报告要求
简述实验原理及目的；
按要求编程求解Yule-Walker方程，并对脑电信号建立AR模型； 使用白噪声驱动生成仿真脑电； 简要回答思考题。
附：参考程序
% 实验四 Yule-Walker方程 clear; clc; load eegdata;
x = x(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('仿真数据');