bt.m

BT法估计功率谱

clear;
Fs=1000;
N=256;
M=64;
t=0:1/Fs:(N-1)/Fs;
%信号的幅度
A1=sqrt(10^(30/10));
A2=A1;
A3=sqrt(10^(27/10));
%产生均值为0，方差为1的复高斯白噪声序列
phi=2*pi*rand;
Noise=randn(1,N)*exp(j*phi);
%三个复正弦信号
S=A1*exp(i*2*pi*t*150)+A2*exp(i*2*pi*t*170)+A2*exp(i*2*pi*t*260);
X=S+Noise; %产生的随机信号
temp1=abs(fft(X,2*N));
temp2=(temp1.*temp1)/N;
r0=ifft(temp2);
rm=[r0(N+1:2*N) r0(1:N)];  %自相关函数
k=1024;
w=-pi:2*pi/k:pi;
%BT
num=rm(N-M:N+M);
Sw_BT=0;
for kk=1:2*M+1
    Sw_BT=Sw_BT+num(kk)*exp(-i*w*(kk-M-1));
end
Sw_BT=10*log10(abs(Sw_BT)/max(abs(Sw_BT)));

plot(w/2*pi,Sw_BT);
title('BT法 M=64,N=256')