AR模型的系统函数H(z)可以表示为：
矩阵形式
可以看出通过矩阵求逆解Yule-Walker方程的运算量很大，通过观察可以看出Yule_Walker方程的自相关矩阵有一系列良好的性质，可以通过递推的方法求p个参数 ，从而避免矩阵求逆的运算，减少运算量。
令p=1，可以得到1阶AR模型Yule-Walker方程：
clearclc% N = 128;% p = 64;N = input('请选择采样点数N:');p = input('请选择AR模型阶数p（建议N/3<p<N/2):');n = 0:N-1;dt = 1;t = n*dt;xn = 2^0.5*sin(2*pi*0.1*t+pi/3) + 10*sin(2*pi*0.13*t+pi/4); %+ randn(1,N);R = xcorr(xn,'biased');% 求xn的自相关函数Rx = zeros(1,p+1);% 预分配内存% 取对称序列的后半部分for i = 1:p+1 Rx(i) = R(i+N-1);endRxRmatl = zeros(p,p);% 生成p个方程对应的矩阵式for m = 1:p for n = 1:p Rmatl(m,n) = Rx(max(m,n)-min(m,n)+1); endendRmatlRmatr = zeros(p,1);for m = 1:p Rmatr(m,1) = -Rx(m+1);end% 输出参数aiai = (inv(Rmatl)*Rmatr)';% (Rmatl\Rmatr)'% 根据自相关函数和ai求解系统增益G = Rx(1);for i = 1:p G = G+ai(i)*Rx(i);endfprintf('系统增益G=%f\n',G^0.5);[H,w] = freqz(G^0.5,[1,ai],N);% 在2*pi范围内N个等分点求系统函数% Px = (G^0.5)*(abs(H)).^2;Hf = abs(H);f = w/(2*pi);plot(w/(2*pi),Hf);[magsor,fsor] = findpeaks(Hf,f,'SortStr','descend');% 降序排列for i = 1:2 fprintf('预测频率f%d=%f\n',i,fsor(i));endclearclc% N = 8;% p = 4;N = input('请选择采样点数N:');p = input('请选择AR模型阶数p（建议N/3<p<N/2):');n = 0:N-1;dt = 1;t = n*dt;xn = 2^0.5*sin(2*pi*0.1*t+pi/3) + 10*sin(2*pi*0.3*t+pi/4);% + randn(1,N);R = xcorr(xn,'biased');% 求xn的自相关函数Rx = zeros(1,p+1);% 预分配内存% 取对称序列的后半部分for i = 1:p+1 Rx(i) = R(i+N-1);end% 求出1阶AR模型的系数和预测误差功率a(1,1) = -Rx(2)/Rx(1);rou(1) = Rx(1)*(1-(a(1,1))^2);for m = 2:p kmsum = 0; for i = 1:m-1 kmsum = kmsum+a(m-1,i)*Rx(m-i+1); end k(m) = -(Rx(m+1)+kmsum)/rou(m-1); a(m,m) = k(m); for i = 1:m-1 a(m,i) = a(m-1,i)+k(m)*a(m-1,m-i); end; rou(m) = rou(m-1)*(1-(k(m))^2);end% rou% plot(rou)xlabel('p');ylabel('rou');G = (rou(p))^2;[H,w] = freqz(G,[1,a(p,1:p)],N);% 在2*pi范围内N个等分点求系统函数% Px = (G^0.5)*(abs(H)).^2;Hf = abs(H);f = w/(2*pi);plot(w/(2*pi),Hf);[magsor,fsor] = findpeaks(Hf,f,'SortStr','descend');% 降序排列for i = 1:2 fprintf('预测频率f%d=%f\n',i,fsor(i));endimport numpy as npimport scipy.signalimport matplotlib.pyplot as pltimport time# 定义了寻找峰值的函数，返回值是元素在数组中的位置def findpeaks(data): length = len(data) weizhi = [] for i in range(1,length-1): if data[i] > data[i-1] and data[i] > data[i+1]: weizhi.append(i) return weizhi# 定义了自相关函数def xcorr(data): length = len(data) Rx = [] datareverse = data[::-1] # 将数组data从后向前取值，每次步进值为1 Rx = scipy.signal.convolve(data,datareverse)/length # 自相关函数 return RxN = input('请选择采样点数N:') # input输入的N是str类型p = input('请选择AR模型阶数p（建议N/3<p<N/2):') # 输入阶数N = int(N) # 把N转换成int类形p = int(p)dt = 1t = np.arange(0,N,dt)# 数组是从0开始的，t=[0 1 2 ... N-1]# 信号频率f1 = 0.1f2 = 0.13x = np.sqrt(2)*np.sin(2*np.pi*f1*t+(np.pi)/3) + 10*np.sin(2*np.pi*f2*t+(np.pi)/4)# 输入信号#print(x)Rx1 = xcorr(x)# 求自相关Rx = Rx1[N-1:N+p]#print('Rx=',Rx)#start = time.clock()#start = time.perf_counter()# 生成自相关系数矩阵Rmatl = np.ones([p,p])for i in range(0,p): for j in range(0,p): Rmatl[i][j] = Rx[max(i,j)-min(i,j)]#print('Rmatl=',Rmatl)Rmatr = np.ones([p,1])for i in range(0,p): Rmatr[i][0] = -Rx[i+1]#print('Rmatr=',Rmatr)# 求系数aiai = np.ones([p,1])Rmatl_inv = np.linalg.inv(Rmatl) # 对左面自相关系数矩阵求逆ai = np.dot(Rmatl_inv,Rmatr) # 得到各个ai的值aiT = ai.T # 求ai的转置#print(aiT)#end = time.clock()#end = time.perf_counter()#print('运行时间=',end-start)# 求系统增益GG1 = Rx[0]for i in range(1,p): G1 = G1 + aiT[0][i]*Rx[i]G = np.sqrt(G1)#print('系统增益G=',G)a = np.insert(aiT,0,[1])# 系统函数分母的系数w,h = scipy.signal.freqz(G,a,worN = N)Hf = abs(h)f = w/(2*np.pi)plt.plot(f,Hf)plt.show()import numpy as npimport scipy.signalimport matplotlib.pyplot as pltimport time# 定义了寻找峰值的函数，返回值是元素在数组中的位置def findpeaks(data): length = len(data) weizhi = [] for i in range(1,length-1): if data[i] > data[i-1] and data[i] > data[i+1]: weizhi.append(i) return weizhi# 定义了自相关函数def xcorr(data): length = len(data) Rx = [] datareverse = data[::-1] # 将数组data从后向前取值，每次步进值为1 Rx = scipy.signal.convolve(data,datareverse)/length # 自相关函数 return RxN = input('请选择采样点数N:') # input输入的N是str类型p = input('请选择AR模型阶数p（建议N/3<p<N/2):') # 输入阶数N = int(N) # 把N转换成int类形p = int(p)dt = 1t = np.arange(0,N,dt)# 数组是从0开始的，t=[0 1 2 ... N-1]# 信号频率f1 = 0.1f2 = 0.3# 输入信号x = np.sqrt(2)*np.sin(2*np.pi*f1*t+(np.pi)/3) + 10*np.sin(2*np.pi*f2*t+(np.pi)/4)#print(x)Rx1 = xcorr(x)# 求自相关Rx = Rx1[N-1:N+p]#print('Rx=',Rx)#start = time.perf_counter()#设定初值a = np.zeros([p,p])a[0][0] = -Rx[1]/Rx[0]rou = np.zeros(p)rou[0] = Rx[0]*(1-np.square(a[0][0]))k = np.zeros(p)#p阶AR模型参数Levinson_Durbin递推算法for m in range(1,p): kmsum = 0 for i in range(0,m): kmsum = kmsum+a[m-1][i]*Rx[m-i] k[m] = -(Rx[m+1]+kmsum)/rou[m-1] a[m][m] = k[m] for i in range(0,m): a[m][i] = a[m-1][i]+k[m]*a[m-1][m-i-1] rou[m] = rou[m-1]*(1-np.square(k[m]))#end = time.perf_counter()#print('运行时间=',end-start)G = np.sqrt(rou[p-1])a = np.insert(a[p-1],0,[1])# 系统函数分母的系数w,h = scipy.signal.freqz(G,a,worN = N)Hf = abs(h)f = w/(2*np.pi)plt.plot(f,Hf)plt.show()联系客服
