do-ar-fun.m

分别使用库函数和自编函数实现AR模型谱估计

%使用周期图法以及改进的周期法实现谱估计
%使用库函数和自编函数实现AR模型谱估计
%
close all;
clear all;
clc;
colordef black;
N=128;
n=1:128;
a1=sqrt(20);
a2=sqrt(2);
w=randn(128,1);
x=a1*sin(2*pi*0.2*n)+a2*sin(2*pi*0.213*n)+w';%数据源序列
figure(1);
avec=arrcest(x,4,1,2,'unbiased');%使用高阶分析工具箱中的函数分析
freqz(1,avec,1024,1);%显示
title('库函数谱估计ls方法');
grid on;
%fft
figure(2);
subplot(211);
yfft=fft(x,N);
plot(n/N,10*log(abs(yfft)));
xlabel('f');
ylabel('dB');
title('周期图法');
grid on;
axis([0,1,-20,70]);
subplot(212);
w=hanning(128)';
plot(n/N,10*log(abs(fft(w.*x))),'g');
xlabel('f');
ylabel('dB');
title('加hanning窗后的周期图法');
grid on;
%----------------------------------------------------------------------------
%xcorr
xxcorr=xcorr(x,'unbiased');
figure(3);
n1=1:2*N-1;
plot(n1/(2*N),10*log(abs(fft(xxcorr))));
xlabel('f');
ylabel('dB');
title('间接谱估计法');
grid on;

%ls
pp=6;%4
qq=20;%qq>=pp
%matrix
for i=1:pp
    for j=1:pp
        recorr(i,j)=xxcorr(qq+i+j-1-pp);%生成自相关序列
    end
end
disp('自相关序列的维数=');
recorrlegnth=size(recorr)
%r
for i=i:pp
    r(i)=xxcorr(qq+i);
end
disp('r的维数=');
rlength=size(r)
r=-r';
yyy=(inv(recorr'*recorr)*recorr')*r;
disp('求解序列的维数');
yylength=length(yyy)
for i=1:yylength
    yyy1(i)=yyy(yylength+1-i);%序列反转
end
figure(4);
freqz(1,yyy1,1024,1);
title('自定义函数实现谱估计');
grid on;