matlab randn函数_功率谱密度函数

最近做时间序列，知道了功率谱密度的概念，matlab里可以直接求，但是却不知道其真正涵义，所以在网上收集了一些资料，然后整理总结，理解的可能不对，部分是他人的原话摘抄。

引入：一个随机震动过程的特征可以用数学期望、方差和相关函数来描述。在工程技术问题中，广泛采用从频率域来描述一个随机振动过程特征的功率谱函数。功率谱密度函数能够反映随机振动的功率关于频率的分布密度。

一、频谱密度，

频谱密度：设一个能量信号为s(t),则它的频谱密度S（w）可以由付氏变换求得。

S(w)=F(s(t))

能量信号的频谱密度S（f）和功率信号C（jnw）（比如一个周期信号）的频谱主要区别有：

（1）S（f）是连续谱，而C（jnw）是离散谱；

（2）S（f）单位是幅度/频率，而C（jnw）单位是幅度；（这里都是指其频谱幅度）

（3）能量信号的能量有限，并连续的分布在频率轴上，每个频率点上的信号幅度是无穷小的，只有df上才有确定的非0振幅；

功率信号的功率有限，但能量无限，它在无限多的离散频率点上有确定的非0振幅

二、功率谱密度

功率谱：也称功率谱密度（PSD），单位是功率/Hz。针对功率有限信号的(能量有限信号用能量谱密度)，所表现的是单位频带内信号功率随频率的变换情况，即信号功率在频域的分布状况。

三、计算方法

1、周期图法：它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。

2、自相关法：根据维纳-辛钦定理，先估计相关函数，再经傅立叶变换得功率谱估计。功率谱与自相关函数是一个傅氏变换对。功率谱具有单位频率的平均功率量纲，所以标准叫法是功率谱密度。通过功率谱密度函数，可以看出随机信号的能量随着频率的分布情况。像白噪声就是平行于w轴，在w轴上方的一条直线。

四、MATLAB实现功率谱函数计算

1、旧版本中较为简单，直接采用psd函数即可

Fs = 1000; t = 0:1/Fs:.296;

x = cos(2*pi*t*200)+randn(size(t));

h = spectrum.welch; % 创建一个韦尔奇谱估计。

Hpsd = psd(h,x,'Fs',Fs); % 计算PSD 值

plot(Hpsd) % 绘制 PSD.图

2、新版本采用很多方法，常用的是pwelch(来自matlab官方说明文档)

[pxx,f] =pwelch(x,window,noverlap,f,fs)returns the two-sided Welch PSD estimates at the frequencies specified in the vector,f. The vector,f, must contain at least 2 elements. The frequencies infare in cycles per unit time. The sample rate,fs, is the number of samples per unit time. If the unit of time is seconds, thenfis in cycles/sec (Hz).

Window, specified as a row or column vector or an integer. Ifwindowis a vector,pwelchdividesxinto overlapping segments of length equal to the length ofwindow, and then multiplies each signal segment with the vector specified inwindow. Ifwindowis an integer,pwelchis divided into segments of length equal to the integer value, and a Hamming window of equal length is used. If the length ofxcannot be divided exactly into an integer number of segments withnoverlapnumber of overlapping samples,xis truncated accordingly. If you specifywindowas empty, the default Hamming window is used to obtain eight segments ofxwithnoverlapoverlapping samples.（窗口函数，窗函数是频谱分析中一个重要的部分，窗函数修正了由于信号的非周期性并减小了频谱中由于泄露而带来的测量不准确性。快速傅里叶变换假定了时间信号是周期无限的。但在分析时，我们往往只截取其中的一部分，因此需要加窗以减小泄露。窗函数可以加在时域，也可以加在频域上，但在时域上加窗更为普遍。截断效应带来了泄漏，窗函数是为了减小这个截断效应，其设计成一组加权系数。窗口函数包含了，1）矩形窗（Rectangle Window） 调用格式：w=boxcar(n)，根据长度 n 产生一个矩形窗 w。2）三角窗（Triangular Window） 调用格式：w=triang(n)，根据长度 n 产生一个三角窗 w。3）汉宁窗（Hanning Window） 调用格式：w=hanning(n)，根据长度 n 产生一个汉宁窗 w。4）海明窗（Hamming Window） 调用格式：w=hamming(n)，根据长度 n 产生一个海明窗 w。5）布拉克曼窗（Blackman Window） 调用格式：w=blackman(n)，根据长度 n 产生一个布拉克曼窗 w。6）恺撒窗（Kaiser Window） 调用格式：w=kaiser(n,beta)，根据长度 n 和影响窗函数旁瓣的β参数产生一个恺撒窗w。）

Noverlap, Number of overlapped samples, specified as a positive integer smaller than the length ofwindow. If you omitnoverlapor specifynoverlapas empty, a value is used to obtain 50% overlap between segments.

rng default
fs = 1000;
t = 0:1/fs:5-1/fs;
x = cos(2*pi*100*t) + randn(size(t));

[pxx,f] = pwelch(x,500,300,500,fs);
plot(f,10*log10(pxx))
xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')

[px,fx]=pwelch(x,Windows,noverlap,nfft,fs);

其中，
x --------- 时程信号；
nfft--------傅里叶变换（fft）点数；
Windows--窗函数，通常设为hanning(nfft)；
nooverlap-重叠的点数，通常为nfft/2；
px---------求得的信号x的功率谱密度；
fx----------对应px的频率序列。

需要说明的是，若中间有个变量不知如何设定（如noverlap），却仍然想设定该变量后面的值，

可采用中括号[]代替，[]则表示该变量并未设定具体值，而是采用函数默认值来取。

【】MATLAB功率谱函数psd与pwelch等价问题_艾米儿_新浪博客

【】[图文]第五章 功率谱密度函数 - 百度文库