功率谱密度的计算

功率谱密度指的是在每个频率上信号具有的功率大小，即 带宽内信号的平均功率，单位为 或者 。

①实数值信号输入

表示信号的功率谱密度。当信号为实数时， 具有对称性。若采样点数为 ，则 关于 对称，为保持总功率不变，可将双边谱中除了第一个点外的其他点的功率谱密度乘2变为单边功率谱密度。

【例】在MATLAB中使用fft()函数求功率谱密度。除此之外，还可以使用序列的自相关函数求功率谱密度（维纳-辛钦定理），在此先不作解释。

fs=1000; %采样频率为1kHz
fin=200; %输入信号为500Hz
N=1000; %采样点数为1000
n=0:1/fs:1-1/fs; %采样点数为1000个
In=sin(2*pi*200*n); %输入信号
In_PSD=(abs(fft(In)).^2)/(fs*N); %求功率谱密度
plot(n*fs,10*log10(In_PSD));
grid on
title('双边功率谱密度')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')


双边功率谱密度

clc;clear;
fs=1000; %采样频率为1kHz
fin=200; %输入信号为500Hz
N=1000; %采样点数为1000
n=0:1/fs:1-1/fs; %采样点数为1000个
In=sin(2*pi*200*n); %输入信号
In_PSD=(abs(fft(In)).^2)/(fs*N); %求功率谱密度
In_PSD(2:1+N/2)=2*In_PSD(2:1+N/2);
plot(n(1:1+N/2)*fs,10*log10(In_PSD(1:1+N/2)));
grid on
title('单边功率谱密度')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')


单边功率谱密度

②加窗减轻频谱泄露

当截断的信号不是整数个周期时，会导致频谱泄露。

频谱泄露 ，就是信号频谱中各谱线之间相互影响，使测量结果偏离实际值，同时在谱线两侧其他频率点上出现一些幅值较小的假谱，导致频谱泄露的原因是 采样频率 和信号频率的不同步，造成周期采样信号的相位在始端和终端不连续。

hann窗的幅值修正系数为2，功率修正系数为1.633。加窗都是通过使截断的边界变得平缓，来减少高频分量以减轻频谱泄露的。

【例】加窗之后可明显看到频谱泄露被减轻了。

clc;clear;
fs=500; %采样频率
fin=50; %输入信号频率
N=128; %采样点数为N
n=(1:N)/fs; %采样点数为N个
In=sin(2*pi*fin*n)+0.3*sin(2*pi*(fin+3)*n); %输入信号
In_PSD=(abs(fft(In)).^2)/(fs*N); %求不加窗功率谱密度
In_PSD(2:1+N/2)=2*In_PSD(2:1+N/2);
subplot(2,1,1);
plot(n(1:1+N/2)*fs*fs/N,(In_PSD(1:1+N/2)));
grid on
title('功率谱密度不加窗')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (V^2/Hz)')
window=hann(N); %汉宁窗
k=1.63; %修正系数 使总功率不变
In_PSD_wind=(abs(fft(1.63*In.*window')).^2)/(fs*N);%求加窗功率谱密度
In_PSD_wind=k*In_PSD_wind;
In_PSD_wind(2:1+N/2)=2*In_PSD_wind(2:1+N/2);
subplot(2,1,2);
plot(n(1:1+N/2)*fs*fs/N,(In_PSD_wind(1:1+N/2)));
grid on
title('功率谱密度加窗')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (V^2/Hz)')


加窗前后的功率谱密度

③归一化的输入频率

输入频率与采样频率的比值叫做归一化频率，也叫做数字频率。

频谱范围变为0-1。

【例】

clc;clear;
fs=500; %采样频率
fin=50; %输入信号频率
f=fin/fs; %归一化输入频率
N=128; %采样点数为N
n=(1:N); %采样点数为N个
In=sin(2*pi*f*n); %输入信号
In_PSD=(abs(fft(In)).^2)/(N*2*pi); %求不加窗功率谱密度
plot([1/N:1/N:1],10*log10(In_PSD(1:N)));
grid on
title('输入频率归一化')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')


具有归一化输入频率的功率谱密度

功率谱的计算

不同于功率谱密度，可从它们的单位中看出不同。功率谱的单位中不包含“ ”,就是功率的单位，它表示在某个频率间隔内的总功率。例如采样率为 ,采样点数为 的离散信号，它的功率谱密度用 表示，那么它的功率谱 为

.

因为在频域中，各条谱线之间的距离为 ，可将这一间隔内的功率谱密度看作相等，计算功率谱。

【例】同功率谱密度，也可以求出双边/单边功率谱，加窗的功率谱，反正都是先求功率谱。

fs=1000; %采样频率为1kHz
fin=200; %输入信号为500Hz
N=1000; %采样点数为1000
n=0:1/fs:1-1/fs; %采样点数为1000个
In=sin(2*pi*200*n); %输入信号
In_PSD=(abs(fft(In)).^2)/(fs*N); %求功率谱密度
In_PS=In_PSD*fs/N;
plot(n*fs,10*log10(In_PS));
grid on
title('双边功率谱')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB)')


双边功率谱

频谱和能量谱的计算

频谱

就是离散傅里叶变换,也有单边/双边频谱的说法。频谱中还包含了相位成分，功率谱密度平方之后将相位丢掉了，不能分析相位。

,

能量谱

帕斯瓦尔定理，信号的总能量既可以按照每单位时间内的能量在整个时间内的积分计算出来，也可以按照每单位频率内的能量在整个频率范围内的积分而得到。

即时域的能量和频域的能量相等。DFT形式的帕塞瓦尔定理为

.

其中 就是能量谱密度。

【例】分别计算方差为1，均值为0的高斯白噪声的能量，结果相等。

N = 2^11;
noise=1e0*randn(1,2048)';
energy_f=sum((abs(fft((noise(1:N))))).^2)/N^2; %频域计算噪声的能量
energy_t=sum((abs((noise(1:N)))).^2)/N; %时域计算噪声的能量

书山有路勤为径，学海无涯苦作舟。