自己编写算法的功率谱密度的三种matlab实现方法(2)

nfft=1024

nfft=10000

相关法：

clear;

Fs=1000; %采样频率

n=0:Fs;%产生含有噪声的序列 nfft=1024;

xn=sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+randn(size(n)); cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数 CXk=fft(cxn,nfft); %求出功率谱密度 Pxx=abs(CXk);

index=0:round(nfft/2-1); f=index/nfft;

plot_Pxx=10*log10(Pxx(index+1)); plot(f,plot_Pxx); xlabel('频率'); ylabel('功率/DB'); grid on;

nfft=256

nfft=1024

Burg法：

clear

Fs=1000 %设置关键变量，可通过调节这些变量观察不同效果 f1=0.2; f2=0.213;

nfft=128; %取样点数

p=50; %阶数p应该选择在N/3

f=0:1/1000:0.5; n=1:Fs;

xn=sin(2*pi*f1*n)+sqrt(2)*sin(2*pi*f2*n)+randn(size(n)); figure;

plot(n,xn); title('burg时域'); xn= xn(:); N=length(xn); ef = xn; eb = xn; a = 1; for l=1:p

efp = ef(2:end);%m-1阶前向预测误差 ebp = eb(1:end-1);%m-1阶后向预测误差

num = -2.*ebp'*efp;%1km分子多项式

den = efp'*efp+ebp'*ebp;%1km的分母多项式

k(l) = num ./ den;%计算反射系数

% 更新前向和后向预测误差

ef = efp + k(l)*ebp;%各阶前向预测误差 eb = ebp + k(l)*efp;%各阶后向预测误差

% 计算模型参数

a=[a;0] + k(l)*[0;conj(flipud(a))];%AR模型参数a end

a1=a(2:p+1);

for i=1:length(f) %循环递推 sum=0; for k=1:p

sum=sum+a1(k)*exp(-m*2*pi*f(i)*k); end

Pbrg(i)=delta/(abs(1+sum))^2;

Pbrg_f(i)=10*log10(Pbrg(i));%求出功率谱 end figure

plot(f,Pbrg_f); title('burg频域');

nfft=128