基于Matlab随机时间序列预测模型数据+源代码

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

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（t2）                      %数据样值点数
randn（‘state‘0）;
x = cos（2*pi*t*200） + randn（1N）;  % 200Hz cosine plus noise

%计算N个取样数据的取样数据自相关函数
rxx = zeros（1N）;   %保存取样数据自相关函数的变量
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（1PMAX+1）; %保存方程系数的中间变量
atemp2 = zeros（1PMAX+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）准则作为AR模型阶次的判断准则      
        a = atemp2（1:（P））;
        P = P - 1;
        deviation = veriance（P-1）;
        break;  
	end
	if （criteria == 2） & （AIC（P-1）         %利用AIC准则作为AR模型阶次的判断准则
        a = atemp2（1:（P））;
        P = P - 1;
        deviation = veriance（P-1）;
        break;
	end
	if（criteria == 3） & （BIC（P-1）        %利用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

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3735  2004-06-21 18:40  ARMODEL.m

     文件       1116  2005-11-13 09:24  zuoye.m

     文件     219136  2005-11-23 12:02  随机序列作业.doc

----------- ---------  ---------- -----  ----

               223987                    3