ar.m

基于ar方法的频谱估计算法

clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 产生信号 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

FS=100; %设采样频率为100，则根据F1/FS=0.2，F2/FS=0.3 OR 0.25 ，可以得到F1，F2便于计算；
N=155; %数据长度 改变数据长度会导致分辨率的变化；
n=0:1/FS:N/FS; 
F1=20; %第一个SIN信号的频率；
F2=25; %第二个SIN信号的频率，取25或者30；
X=sin(2*pi*F1*n)+sin(2*pi*F2*n)+(1/20)*randn(size(n)) ; %产生信号，由信噪比为10dB推出噪声功率；信号长度从X（1）到X（N+1）；

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 产生自相关函数数组 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for m=1:N+1 %初始化R（m），R（m）用来存放自相关函数；由于R（0）在MATLAB里无效，所以从1开始到N+1；
R(m)=0;
end
for m=1:N+1 %做自相关函数的无偏估计；
S=0; 
for n=1:N+2-m
H=X(n)*X(n+m-1);
S=S+H;
end
R(m)=S/N;
end %估计完毕 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Levinson算法主体部分 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k=1:N+1 %初始化后面算法中要用到的数组，其中a(k,k)用来存放AR模型参数，FPE(k)是最终预测误差准则函数，
a(k,k)=0; % U2(k)是AR模型中的另一参数即误差功率；
FPE(k)=0;
U2(k+1)=0;
end
U2(1)=R(1); 
a(1,1)=-R(2)/R(1); %由Levinson算法，计算一阶模型参数a11；
U2(2)=(1-(abs(a(1,1)))^2)*U2(1); %由Levinson算法，计算一阶模型参数 误差功率；
FPE(1)=U2(2)*(N+2)/N; %计算一阶模型的最终预测误差准则函数；
S=0;
for k=2:N %由Levinson梯归公式计算K阶模型参数和FPE函数；
for l=1:k-1
M=a(k-1,l)*R(k-l+1);
S=S+M;
end
a(k,k)=-(R(k+1)+S)/U2(k);
for i=1:k-1
a(k,i)=a(k-1,i)+a(k,k)*a(k-1,k-i);
end
U2(k+1)=(1-(abs(a(k,k)))^2)*U2(k);
FPE(k)=U2(k+1)*(N+k+1)/(N-k+1);
S=0;
end %梯归计算完毕；
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 确定阶数P %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
min=FPE(1); %求出使得FPE函数取最小值的阶数P；
for k=2:N
if FPE(k)<min
min=FPE(k);
p=k;
end
end 
%p=2 %为了调整效果可以在这里自行指定阶数；
disp('输出模型参数a');
for k=1:p
disp(a(p,k));
end
disp('输出误差功率');
disp(U2(p+1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% AR模型参数确定后计算出功率谱 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Z=0;
W=0:0.01:pi; %功率谱以2*pi为周期，又信号为实信号，只需输出0到PI即可；
for k=1:p 
Z=Z+(a(p,k).*exp(-j*k*W));
end
PXX=U2(p+1)./((abs(1+Z)).^2); %得到功率谱函数；
F=W*FS/(pi*2); %将角频率坐标换算成HZ坐标，便于观察；
plot(F,abs(PXX))
disp('输出阶数P')
disp(p);
grid