BT法和周期图法估计自相关函数

clc,clear
%产生零均值、方差为1的复高斯白噪声序列
N=256;
noise=(randn(1,N)+1i*randn(1,N))/sqrt(2);
%产生三个复正弦信号
f1=0.15;
f2=0.17;    %复正弦信号归一化频率
f3=0.26;
SNR1=30;
SNR2=30;    %复正弦信号信噪比
SNR3=27;
A1=10^(SNR1/20);
A2=10^(SNR2/20);     %复正弦信号幅度
A3=10^(SNR3/20);
signal1=A1*exp(1i*2*pi*f1*(0:N-1));
signal2=A2*exp(1i*2*pi*f2*(0:N-1));     %复正弦信号
signal3=A3*exp(1i*2*pi*f3*(0:N-1));
%产生观察样本
un=signal1+signal2+signal3+noise;
%周期图法
NF=1024; %周期图法中FFT的点数
Spr=fftshift((1/NF)*abs(fft(un,NF)).^2); %计算信号的周期图
%BT法
M=64;%自相关函数的单边长度
r=xcorr(un,M,'biased');%计算自相关函数
NF=1024;%BT法中FFT的点数
BT=fftshift(fft(r,NF));%BT法计算功率谱
%作图
w=-(NF-1)/(2*NF):1/NF:(NF-1)/(2*NF);
figure(1)