如下的一段程序,其中pr是spectrum.welch计算出来的功率谱密度
Fs=2000;
a=sin(2*pi*50*(1:1000)/Fs);
h = spectrum.welch;
Hpsd = psd(h,a(500:1000),'Fs',Fs,'NFFT',512);
pr=10*log10(Hpsd.data);
以下是波叠加法,功率谱密度生成时间序列的程序
NFFT=512;
t=(1:1000)/Fs;
num=length(t);
NN=length(pr);
deta=1000/(NFFT/2);
omega=deta*(0:NFFT/2);
theta=2*pi*rand(1,NN);
for j=1:num
for i=1:NN
ff(i)=sqrt(2*pr(i)*deta)*cos(2*pi*omega(i)*t(j)-theta(i));
end
f(j)=sum(ff);
end