Matlab高斯白噪声通过互相关函数求功率谱密度

参考文章：百度文库文章
但是该文章有很多部分，不知是matlab更新了还是咋地，有错误，踩了很多坑，仅做参考。

产生高斯白噪声的方法在“信号处理”专栏中有写，其时域和频域图如下：

用matlab中求自（互）相关的xcorr函数，参考帮助文档

[r,lags]=xcorr(z,'biased');

其中z是上述高斯白噪声，r是自相关函数，lags是时间偏移量（索引），尤其注意’biased’参数，这是调试了半天才发现的问题。
帮助文档里描述如下：

相当于算出来的r除以了样本数，这是离散自相关函数计算步骤。
如果不加这个参数，自相关函数值会变得很大。
算出自相关函数之后，根据维纳辛钦定理，信号的功率谱密度和信号的自相关函数是一个傅里叶变换对。使用fft函数求自相关函数的傅里叶变换：

pdv=fft(r);
pdv=abs(fftshift(pdv))./length(pdv);

fftshift是频谱矫正函数，将fft后的结果以中心为分界，两边分别做镜像，得出来的结果再除以样本数，即双边频谱图。
而这个频谱图就是原信号的功率谱密度。

plot((0:length(pdv)-1)*L/length(pdv)-L/2,pdv)

注意这里横坐标是频率，需要做映射处理，L为样本数量

高斯白噪声其实看不出来啥，下面用一个没加噪声的纯信号来展示：

y=12*cos((2*pi)*100.*t)+15*cos((2*pi)*150.*t)+18*cos((2*pi)*210.*t)+10;

上面的功率谱是由双边频谱的平方（除直流）计算出来的，下面的是用维纳辛钦定理计算出来的，可以看到功率谱几乎相同。
信号的自相关函数如下：

编者水平有限，很多知识点都是似懂非懂，如有错误欢迎指出！
附代码：

fs=1000;%采样频率hz
T_N=1.5;%总时间s
t=1/fs:1/fs:T_N;%时间向量
L=T_N*fs;%样本长度
y=12*cos((2*pi)*100.*t)+15*cos((2*pi)*150.*t)+18*cos((2*pi)*210.*t)+10;%信号
subplot(2,1,1);
plot(t,y);
xlabel("时间/s")
ylabel("幅度/v")
title("时域")
fft_y=fft(y);%快速傅里叶变换
P = abs(fft_y/L);%取幅频特性，除以L
P = P(1:L/2+1);%截取前半段
P(2:end-1)=2*P(2:end-1);%单侧频谱非直流分量记得乘以2
f = fs*(0:(L/2))/L;%频率，最多到一半（奈奎斯特采样定理）
subplot(2,1,2);
plot(f,P);
xlabel("频率/Hz")
ylabel("幅度/v")
title("单边频谱")

figure(4)
[r,lags]=xcorr(y,'biased');%得到自相关函数的幅度和偏移量
subplot(2,1,1)
plot(lags,r)
xlabel("时间偏移/s")
ylabel("相关程度")
title("（自相关函数）")

fft_y=fftshift(fft_y);%频谱矫正
powerpu=abs(fft_y./L).^2;%双边功率谱
% subplot(2,1,1)
% plot((1:size(powerpu,2))-751,powerpu)
% xlabel("频率/Hz")
% ylabel("功率/W")
% title("功率谱")

pdv=fft(r,size(r,2));%对自相关函数快速傅里叶变换
pdv=abs(fftshift(pdv));%频谱矫正，让正半轴部分和负半轴部分的图像分别关于各自的中心对称，得到双边谱
subplot(2,1,2)
plot((1:length(pdv))*L/length(pdv)-L/2,pdv./length(powerpu))
xlabel("频率/Hz")
ylabel("功率谱W/Hz")
title("（功率谱）")

% z1=0.1*randn(1,201);%产生方差N(0,0.12)高斯白噪声
% [r1,lags]=xcorr(z1,'unbiased');%自相关函数的估计
% plot(lags,r1);
% f1=fft(r1);
% f2=fftshift(f1);%频谱校正
% l1=((0:length(f2)-1)*200/length(f2)-100);%功率谱密度x轴
% y4=abs(f2);
% figure(2)
% plot(l1,y4);