参考文章：百度文库文章
但是该文章有很多部分，不知是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;%采样频率hzT_N=1.5;%总时间st=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);%取幅频特性，除以LP = P(1:L/2+1);%截取前半段P(2:end-1)=2*P(2:end-1);%单侧频谱非直流分量记得乘以2f = 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);