autocorr_signal_noise = xcorr(signal_noise); %求线性自相关
fft_autocorr_signal_noise = fft(autocorr_signal_noise); %求自相关的fft
abs_fft_autocorr_signal_noise = abs(fft_autocorr_signal_noise)/(2*points-1); %求幅值
max_abs_fft_autocorr_signal_noise = max(abs_fft_autocorr_signal_noise);
abs_fft_autocorr_signal_noise = abs_fft_autocorr_signal_noise/max_abs_fft_autocorr_signal_noise;
log10_abs_fft_autocorr_signal_noise = 10*log10(abs_fft_autocorr_signal_noise);
Delta_f1 = (0:2*points-2)*fs/(2*points-1);
figure;
plot(Delta_f1(1:points)/1000,log10_abs_fft_autocorr_signal_noise(1:points));
xlabel('KHz');
ylabel('归一化功率谱P(k)/dB');
string = ['间接法（自相关法BT）信号点数=',num2str(signal_points),',采样点数=',num2str(points),...
',信号频率=',num2str(f0/1000),'KHz,采样率=',num2str(fs/1000),'KHz,信噪比SNR=',...
num2str(SNR),'dB'];
title(string);
Bartlett法的Matlab程序
首先将截断的信号分成等长的不重叠的若干段，然后计算每一段的离散傅里叶变换，并取模值的平方，然后对应点求和。所有段计算并求和完成后，归一化处理，并取以10为底的对数然后乘以10，绘制归一化后的功率谱图。
M = 128; %每段的点数
L = points/M; %段数
Bartlett_Spectrum = zeros(1,points);
for i=1:L
temp_signal = zeros(1,points); %清零
start_point = (i-1)*M;
temp_signal(1:M) = signal_noise((start_point+1):(start_point+M)); %截取第i段信号
fft_temp_signal = fft(temp_signal);
squar_fft_temp_signal = fft_temp_signal.*conj(fft_temp_signal); %求当前段的功率谱
Bartlett_Spectrum = Bartlett_Spectrum + squar_fft_temp_signal; %累加
end
max_Bartlett_Spectrum = max(Bartlett_Spectrum);
Bartlett_Spectrum = Bartlett_Spectrum/max_Bartlett_Spectrum; %归一化
log10_Bartlett_Spectrum = 10*log10(Bartlett_Spectrum); %求对数
Delta_f2 = (0:points-1)*fs/points;
figure;
plot(Delta_f2(1:points/2)/1000,log10_Bartlett_Spectrum(1:points/2));
xlabel('KHz');
ylabel('归一化功率谱P(k)/dB');
string = ['Bartlett法,信号点数=',num2str(signal_points),',采样点数=',num2str(points),...
',信号频率=',num2str(f0/1000),'KHz,采样率=',num2str(fs/1000),'KHz,信噪比SNR=',...
num2str(SNR),'dB,段数L=',num2str(L),',每段点数M=',num2str(M)];
title(string);
Welch法的Matlab程序
首先将截断的信号分成等长的有重叠的若干段，然后计算每一段的离散傅里叶变换，并取模值的平方，然后对应点求和。所有段计算并求和完成后，归一化处理，并取以10为底的对数然后乘以10，绘制归一化后的功率谱图。
M = 128; %每段的点数
Overlap = M-M/16; %重叠点数
L = (points - Overlap)/(M-Overlap); %段数
Welch_Spectrum = zeros(1,points);
for i=1:L
temp_signal = zeros(1,points); %清零
start_point = (i-1)*(M-Overlap);
temp_signal(1:M) = signal_noise((start_point+1):(start_point+M)); %截取第i段信号
fft_temp_signal = fft(temp_signal);
squar_fft_temp_signal = fft_temp_signal.*conj(fft_temp_signal); %求当前段的功率谱
Welch_Spectrum = Welch_Spectrum + squar_fft_temp_signal; %累加
end
max_Welch_Spectrum = max(Welch_Spectrum);
Welch_Spectrum = Welch_Spectrum/max_Welch_Spectrum; %归一化
log10_Welch_Spectrum = 10*log10(Welch_Spectrum); %求对数
Delta_f3 = (0:points-1)*fs/points;
figure;
plot(Delta_f3(1:points/2)/1000,log10_Welch_Spectrum(1:points/2));
xlabel('KHz');
ylabel('归一化功率谱P(k)/dB');
string = ['Welch法,信号点数=',num2str(signal_points),',采样点数=',num2str(points),...
',信号频率=',num2str(f0/1000),'KHz,采样率=',num2str(fs/1000),'KHz,信噪比SNR=',...
num2str(SNR),'dB,段数L=',num2str(L),',每段点数M=',num2str(M),',重叠点数=',...
num2str(Overlap)];
title(string);
Matlab运行结果及分析
设置相应参数运行上述程序得到如下图：
图1信号的时域图
图2周期图法得到的功率谱图 图3 BT法得到的功率谱(M=128)
图4 Batlett法得到的功率谱（4段） 图5
Welch法得到的功率谱(128点)
图6 BT法得到的功率谱(M=64) 图7 BT法得到的功率谱(M=32)
图8 Batlett法得到的功率谱（8段） 图9 Welch法得到的功率谱(64点)
图1为生成的仿真信号，一共512个点。
图2为用周期图法对数据直接求出功率谱。由于主瓣宽度B = 2/512 = 0.003906 < (150-140)/750 = 0.013所以f0与f1能够分开。
图3是利用BT法求出的功率谱，M=128,依然能将f0与f1分开；对于与图2其主瓣更宽了。当M=64时，如图6 f0与f1不能完全分开，只是在波形的顶部看出是两个频率的分量。当M=32,如图7 f0与f1就已经无法分开了。
图4是利用Batlett法求出的功率谱，共分4段，每段128个点，可以看到f0与f1能够被分开，相对于图2其更加平滑，但是主瓣的宽度比图2的要宽。当分为8段，每段64个点时，如图8其变的更加平滑但f0与f1已经不能完全分开了。
图5是利用Welch方法求出的平均周期图，每段128个点，重合120个点，谱变得更加的平滑，效果与图4基本一致。若每段64个点，重合60个点，如图9
f0与f1已不能分开，效果与图8基本一致。