AR过程功率谱估计MATLAB实现代码带备注，没吊用函数，自己写的

clc
close all
clear all


 rng（1）;%设置种子 产生固定白噪声（详见帮助文档）

N=1024;%产生高斯白噪声方法1
pee=2;%功率
e=wgn（1Npee）;
pw=sum（abs（e）.^2）/N;   % 验证功率
mean=sum（e）/N;         %验证均值
d=var（e）;                %验证方差

% %产生任意均值和方差白噪声方法 2
% % rand函数产生的是幅值为[0，1]均匀分布的白噪声，均值为1/2，方差为1/12
% % 使用下面的公式可以产生指定均值和方差的均匀分布的白噪声
% v=2;%方差
% m=0;%均值
% white_noise =  （rand（1N） - 0.5 ） * sqrt（ 12 * v ） + m;%给定均值与方差的白噪声
% mea=mean（white_noise）;
% va=var（white_noise）;
% Q=sum（abs（white_noise）.^2）/N;

[ba]=cheby2（1080[380 400]/750‘bandpass‘）;%滤波器
[HH2W22]=freqz（ba）;
x=filter（bae）;%AR建模产生样本信号x（n）
figure（1）
plot（W22/pi10*log10（abs（HH2）））;
title（‘窄带滤波器幅频响应‘）;
xlabel（‘角频率 （x \pi）‘）;
ylabel（‘幅度（db）‘）;


figure（2）
subplot（121）
plot（1:Ne）;
title（‘方差为2的1024点高斯白噪声‘）
xlabel（‘个数‘）;
ylabel（‘幅度（db）‘）;
subplot（122）
plot（1:Nx）;
title（‘AR建模产生信号x（n）‘）
xlabel（‘个数‘）;
ylabel（‘幅度（db）‘）;

x1=fft（x）;%直接法估计
x11=abs（x1）.^2;
x111=10*log10（（x11）/N）;
N1=length（x111）;


pow=0;
for ii=1:length（x）
    pow=pow+（abs（x（ii）））.^2;
end
pow111=pow/length（x）;
pow1111=10*log10（pow111）;%信号功率（时域）

%AR模型估计功率谱
%点数为1024，假设pe为300
r=zeros（1N）;
M=500;%M>>pe  M 的选择参见书130
pe=300;