内容发布更新时间 : 2025/7/18 17:27:02星期一 下面是文章的全部内容请认真阅读。
noise = randn(1,len);% 产生均值为0,方差(功率)为1,数据长度为3000的高斯白噪声
a0=[1 0.72 0.88];% 已知模型的输入系数 x=filter(1,a0,noise);%
产
生
的
信
号
模
型
是
x(n)=a0(1)*noise(n-2)+a0(2)*noise(n-1)+a0(3)*noise(n);
%% 利用L-D算法预测模型参数 m = 50;
[a,P] = Levinson_Durbin(x,m); aa = [ones(m,1),a];% 预测出来的系数 %预测误差功率图 figure
plot(1:m,P(1:m),'rh') ylim([0,6])
title('预测误差功率随迭代次数的变化'); xlabel('迭代次数') ylabel('预测误差功率') % 利用aryule验证 y = aryule(x,m); figure subplot(1,2,1) plot(1:51,aa(m,:),'r')
title('编写程序运行得出的系数') subplot(1,2,2) plot(1:51,y,'k')
title('aryule运行得出的系数')
6
4.2、LD算法子函数: %% Levinson-Durbin algorithm function [a,P]=Levinson_Durbin(x,m) % x为输入信号 % m为阶数 a = zeros(m,m); N = length(x); P = zeros(1,m + 1); for i = 1:N
P(1) = P(1) + x(i)^2; end
P(1) = P(1)/N;% 计算预测误差功率的初始值P_0
f = zeros(m + 1,N);% m阶前向预测误差 g = zeros(m + 1,N);% m阶后向预测误差 f(1,:) = x; g(1,:) = x; K = zeros(1,m); % 开始迭代: for i = 1:m numerator = 0; denominator = 0; for j = i + 1:N
numerator = numerator - f(i,j)*g(i,j - 1);
denominator = denominator + (f(i,j)^2 + g(i,j - 1)^2)/2; end
K(i) = numerator/denominator;% 反射系数 for j = 1:i - 1
7
a(i,j) = a(i - 1,j) + K(i)*a(i - 1,i - j); end a(i,i) = K(i); if i >= 2
P(i) = (1 - K(i)^2) * P(i - 1);% 更新预测误差功率 end for j = i + 1:N
f(i + 1,j) = f(i,j) + K(i)*g(i,