armodel.m

功率谱估计的应用范围很广

function ARMODEL()
%现代数字信号处理作业
%采用MATLAB仿真实现AR模型，进行谱估计
%AR模型的理论公式： x(n) + a1*x(n-1) + a2*x(n-2) + …… +ap*x(n-p) = w(n)
%待估计数据样本：
%目前输入随机过程为叠加了正态分布噪声的正弦波

Fs = 1000;  
t = 0:1/Fs:15;
N = size(t,2)                      %数据样值点数
randn('state',0);
x = cos(2*pi*t*200) + randn(1,N);  % 200Hz cosine plus noise

%计算N个取样数据的取样数据自相关函数
rxx = zeros(1,N);   %保存取样数据自相关函数的变量
for m = 0:N-1
    sum = 0;
    for n= 1:N-m
        temp1 = x(n)*x(m+n);
        sum = sum + temp1;
    end
    rxx(m+1) = sum/N;
end

%采用Levison-Durbin算法求解AR模型的Yule-Walker模型
%需要确定AR模型理论公式中的参数：白噪声w(n)的方差、方程系数a1……ap（这里包括了模型的阶次）
PMAX = 100;                %设定AR模型最高阶次
atemp1 = zeros(1,PMAX+1); %保存方程系数的中间变量
atemp2 = zeros(1,PMAX+1); %保存方程系数的中间变量
deviationtemp1 = zeros; %保存白噪声w(n)方差的中间变量
deviationtemp2 = zeros; %保存白噪声w(n)方差的中间变量

%AR(1)模型：x(n) + a1*x(n-1) = w(n)
%其Yule-Walker方程： R(0)*1 + R(1)*a1 = deviation1;
%                    R(1)*1 + R(0)*a1 = 0;
%求解方程确定a1、deviation1
atemp1(1) = 1;
atemp1(2) = -rxx(2)/rxx(1);
atemp2 = atemp1;
deviationtemp1 = (    rxx(1)*rxx(1)   -   rxx(2)*rxx(2)   )/rxx(1);
deviationtemp2 = deviationtemp1;

%利用Levison-Durbin迭代算法计算AR模型参数
%根据FPE准则、AIC准则和BIC准则确定AR模型的阶次
%atemp1、deviation1保存第k次的运算结果
%atemp2、deviation2保存第k+1次的运算结果
FPE(1) = deviationtemp1*(N+2)/N;
AIC(1) = log(deviationtemp1) + 2/N;
BIC(1) = log(deviationtemp1) + log(N)/N;
veriance(1) = deviationtemp1;

criteria = 3 

for P = 2:PMAX
    sum1 = 0;
    sum2 = 0;
    for i = 2:(P+1)
        sum1 = atemp1(i)*rxx(i) + sum1;
    end
    
    for i = 1:(P+1)
        sum2 = atemp1(i)*rxx(P+2-i) + sum2;
    end
    
    deviationtemp1 = rxx(1) + sum1;
    dk = sum2;
    ref(P) = dk/deviationtemp1;
    deviationtemp2 = ( 1 - ref(P)*ref(P) )*deviationtemp1;
    
    for i = 2:(P+1)
    atemp2(i) = atemp1(i) - ref(P)*atemp1(P+2-i);
    end
    %计算AR(P)模型参数
    atemp1 = atemp2;
    veriance(P) = deviationtemp2;
    
    %利用阶次判断准则
    FPE(P) = veriance(P)*(1+P/N)/(1-P/N);
    AIC(P) = log(veriance(P)) + 2*P/N;
    BIC(P) = log(veriance(P)) + P*log(N)/N;
    
	if (criteria == 1) & (FPE(P-1)<FPE(P))
        %利用最终预测误差（FPE）准则作为AR模型阶次的判断准则      
        a = atemp2(1:(P));
        P = P - 1;
        deviation = veriance(P-1);
        break;  
	end
	if (criteria == 2) & (AIC(P-1)<AIC(P))
         %利用AIC准则作为AR模型阶次的判断准则
        a = atemp2(1:(P));
        P = P - 1;
        deviation = veriance(P-1);
        break;
	end
	if(criteria == 3) & (BIC(P-1)<BIC(P))
        %利用BIC准则作为AR模型阶次的判断准则
         a = atemp2(1:(P));
         P = P - 1;
         deviation = veriance(P-1);
         break; 
    end
end

%计算功率谱估计值
a = atemp2(1:(P+1));
deviation = veriance(P);
freqsample = 10;
NF = Fs/freqsample;
freq = [freqsample:freqsample:Fs];
delta_w = 2*pi/NF;
theta = [delta_w:delta_w:2*pi]; 
for k = 1:NF
    sum = 0;
    for i= 2:P+1
        theta1 = -1*theta(k)*(i-1);
        sum = a(i)*(    cos(theta1)+j*sin(theta1)   ) + sum;
    end

    Sxx(k) = deviation /(   (1+sum)*conj(1+sum) );
end
Sxx = 10*log10(Sxx/Fs);

P
veriance
deviation

%画出功率谱图
plot(freq,Sxx,'r-');
grid on
title('Yule-Walker 功率谱估计')
xlabel('频率 f/Hz')
ylabel('功率谱密度PSD dB/Hz')
legend('功率谱密度估计曲线')

%利用MATLAB内部函数考察AR模型的频率响应
figure
freqz(1,a)  % AR filter frequency response
title('AR模型频率响应')