matlab谱估计

实验目的：

一、.理解参数模型谱估计的基本思路和方法优缺点。

二、掌握利用MATLAB编程进行AR模型谱估计的基本方法、计算步骤和函数使用。

1.产生一个多点频信号，利用上述参数谱估计方法进行功率谱估计；对比经典谱估计与参数谱估计几种方法的处理结果，并进行分析。

程序代码及编程思路：

clear all;

N=128;

p=40;

NFFT=2048;

Fs=2;

n=0:N-1;

randn(‘state’,0);%保持随机数不变

x=cos(0.3pin)+cos(0.4pin)+cos(0.44pin)+randn(size(n));

%产生3个点频的信号

[P,f]=periodogram(x,[],NFFT,2);%周期图法

[Py,fy]=pyulear(x,p,NFFT,2);%AR模型

[Px,F]=pwelch(x,hamming(65),overlap,NFFT);%Welch法

%使用汉明窗，叠合32点，FFT点数为2048

subplot(311);

plot(f,10*log§);

grid on;

title(‘周期图法’);

axis([0 1 -60 60]);

subplot(312);

plot(fy,10*log(Py));

grid on;

title([‘Yule-Walker方法，阶次为’,num2str§]);

axis([0 1 -60 60])

subplot(313)

plot(F,Px);

title([‘Welch法，汉明窗，叠合点数为’,num2str(overlap)])

%%以下函数是通过逐个计算参数进行的

r0=xcorr(x’,p,‘biased’);%计算自相关

r=r0(p+1:2p+1)’;%自相关的个数2p+1个，取后一半

a(1,1)=-r(2)/r(1);%教材14.2.14a公式

sigma(1)=r(1)-(abs(r(2))^2)/r(1);%教材14.2.14b公式

%备注：因为在matlab中，向量、矩阵起始点都是从1开始的

%所以教材上的sigma(0)即为本代码中的sigma(1)

%下方代码中反射系数k也是同理

for m=2:p%教材14.2.15公式

k(m)=-(r(m+1)+sum(a(m-1,1:m-1).*r(m👎2)))/sigma(m-1);

a(m,m)=k(m);

for i=1:m-1

a(m,i)=a(m-1,i)+k(m)*conj(a(m-1,m-i));

end

sigma(m)=sigma(m-1)*(1-abs(k(m))^2);

end

par=sigma§./(abs(fft([1,a(p,:)],NFFT)).^2);%计算AR模型谱估计

figure;

ff =
-0.5:1/NFFT:0.5-1/NFFT;%归一化

plot(ff,10*log(fftshift(par)));

title(‘使用详细步骤逐参数求解方式的AR模型图’)