时间序列分析的一些模型Matlab源码

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模型频率响应')