基于BT方法的经典功率谱估计仿真

这是一份基于BT方法的经典功率谱估计仿真作业
这里附上实验代码：

// 基于BT方法的经典功率谱估计仿真
%步骤一：产生随机过程信号
A1 = sqrt(2000);
A2 = sqrt(2000);
A3 = sqrt(2.*10^2.7);
%unifrnd(a,b) 生成1个服从U(a,b)均匀分布的随机数
Q1 = unifrnd(0,2.*pi);
Q2 = unifrnd(0,2.*pi);
Q3 = unifrnd(0,2.*pi);
n = 0:1:255;
S1 = A1.*cos(0.2.*pi.*n+Q1);
S2 = A2.*cos(0.5.*pi.*n+Q2);
S3 = A3.*cos(0.54.*pi.*n+Q3);
%这里S1+S2+S3是原始图像，添加了均值是0，方差是1的高斯噪声。这里需要对高斯噪声进行归一化（方法一）
%Vn = imnoise(S1+S2+S3,'gaussian',0,1);%引入高斯噪声
%Un = S1+S2+S3+Vn;
Nnoise = 256;
% 对高斯白噪声，其方差和功率（单位为W）是一样的。（方法二）
% 因此，对方差，要做的只是将w变换成dbw，即dbw=10log10（w）
Vn = wgn(Nnoise,1,10*log10(1));%256*1
Un = S1+S2+S3+Vn';%1*256
%随机信号生成定义域[1~256]的256点信号
%步骤二：自相关函数
N = 256;
Rm = xcorr(Un,255,'biased');
%未截短后自相关函数生成定义域[1~511]的511点信号
Rm_small = Rm(192:1:320);
%保留中间，截取两端，截短后自相关函数生成定义域[1~129]的129点信号
%步骤三：计算信号的功率谱
NF=1024;%BT法中FFT的点数
BT=fftshift(fft(Rm_small,NF));%BT法计算功率谱
w=-(NF-1)/(2*NF):1/NF:(NF-1)/(2*NF);%保证0~1之间有1024个数，1023个间隔
subplot(121)
BT = 10*log10(abs(BT)/max(abs(BT)));%归一化功率谱
plot(w,BT,'b');
axis([-0.5 0.5 -40 0]);
xlabel('{\it\omega}/2{\pi}');
ylabel('归一化功率谱/dB');
line([0.1 0.1],[-40 0],'linestyle','--', 'Color','r', 'LineWidth',1);
line([0.25 0.25],[-40 0],'linestyle','--', 'Color','r', 'LineWidth',1 );
line([0.27 0.27],[-40 0],'linestyle','--', 'Color','r', 'LineWidth',1);

做的比较粗糙，对一些函数只是使用，对其内部的参数还不太了解其中的差异。
实验结果：