信号处理—频域分析（功率谱、功率谱密度）

前言

书接上回，幅值谱和相位谱是信号频域分析的一个重要手段。在这其中，幅值谱和相位谱主要针对满足傅里叶变换条件的信号，比如能量信号。而工程和实际中接触的信号往往都是随机信号，随机信号不满足FFT变换的条件（在宋老师的书中（见参考链接）提到随机信号不满足傅里叶变换绝对可积条件），不能直接进行FFT变换，顾一般不采用幅值谱和相位谱进行分析。在目前故障诊断领域，对相位谱的关注不多，所以功率谱可以看出对标幅值谱的存在。

关于频域分析中的一些频谱、能量谱、功率谱、功率谱密度、自功率谱、互功率谱、随机信号的定义等等，请跳转下面这条链接了解：

信号处理（频谱，能量谱，功率谱）--学习笔记

https://blog.csdn.net/weixin_45647721/article/details/128008836

友情提示：上面链接的内容无法保证一定正确，请根据自己理解，自行判断。不过这是冷漠看了这么多科普，说的比较清晰、比较准确的一个。

1.内容概括

本文介绍信号频域分析中功率谱和功率谱密度（PSD）的相关代码和分析过程！

请注意该功率谱和功率谱密度采用的是非参数法中的周期图法！周期图法！周期图法！！！

虽然本文提供了功率谱和功率谱密度的封装函数，但是推荐用Matlab的periodogram函数，本文只是重现了这个函数的结果，让大家更加透彻的理解。

除此之外，估计PSD的方法matlab提供了很多，如参数法、子空间方法等等。Matlab对应的函数具体如下：

有需要的小伙伴自行阅读哈！！！

该内容参考了一些资料：

1、书籍：MATLAB数字信号处理85个实用按比例——入门到进阶  宋知用 编著

2、书籍：从这里学NVH噪声、振动、模态分析的入门与进阶   谭祥军 编著

3、书籍：机械工程测试技术基础 第3版 熊诗波编著

4、matlab官网的一些例子和教程：

1）频域分析实践介绍

https://ww2.mathworks.cn/help/signal/ug/practical-introduction-to-frequency-domain-analysis_zh_CN.html#d126e17836

2）使用 FFT 获得功率频谱密度估计

https://ww2.mathworks.cn/help/signal/ug/power-spectral-density-estimates-using-fft.html

3）频谱估计

https://ww2.mathworks.cn/help/signal/nonparametric-spectral-estimation.html

5、了解功率谱密度和功率谱

https://www.bilibili.com/video/BV1ui421Z7XP/?spm_id_from=333.337.search-card.all.click&vd_source=50b116844fc6444ae380b4347979702f（一个涵盖了幅值谱、功率谱和功率谱密度的视频，强烈推荐！！！）

6、一些有价值的科普（公众号搜索），理论讲述更加全面

1）信号处理进阶之路，文章名称：能量域特征提取方法https://mp.weixin.qq.com/s/W-qsi0e1aWC3pz3c4dy8kA

2）看海的城堡，文章名称：信号频域分析方法的理解（频谱、能量谱、功率谱和倒频谱）：https://mp.weixin.qq.com/s/RpCoRAiZlo96wWYgW09lnw

3）声振之家，文章名称：几种经典功率谱估计方法的实现（Matlab）及其局限性https://mp.weixin.qq.com/s/xvYS-YnUtlXntRsKXAkkMQ

代码采用了Matlab 2024a进行运行，欢迎大家测试和提出问题！

2.理论铺垫

功率谱和功率谱密度是针对随机信号而言，那么首先要认清什么是随机信号？

借用熊老师关于信号分类的定义，如下图：

确定性信号：可以表示为一个确定的时间函数，因而可确定其任何时刻的量值，这种信号称为确定性信号。

周期信号：按照一定时间间隔周而复始重复出现，无始无终的信号。

非周期信号：不具有周期重复性的信号。

准周期信号：两种以上的周期信号合成的，但其组成分量间无法找到公共周期，因而无法按某一时间间隔周期而重复出现的信号。

瞬变非周期信号：除准周期信号之外的其他非周期信号，是一些或在一定时间区间内存在，或随着时间的增长而衰减至零的信号。

随机信号：一种不能准确预测其未来瞬时值，也无法用数学关系式来描述的信号。但是，它具有某些统计特征，可以用概率统计方法由其过去来估计其未来。如噪声信号、含噪信号，工程中没有绝对意义的无噪信号，所以都是随机信号，所以都用功率谱和功率谱密度分析（个人理解）。

随机信号扩展：

对随机信号按照历程所作的各次长时间观测记录成为样本函数。在同一试验条件下，全部样本函数的集合就是随机过程。随机过程有平稳过程和非平稳过程之分。平稳随机过程是指其统计特征参数不随时间而变化的随机过程，否则为非平稳随机过程。

在平稳随机过程中，若任一单个样本函数的时间平均统计特征等于该过程的集合平均统计特征，这样的平稳随机过程叫各态历经（遍历性）随机过程。工程中，很多随机信号分析时我们都前提认为所取信号具有各态历经性（大家默认，都不强调），有的虽然不满足严格的各态历经过程，但可以忽略那些不具有各态历经性分量，它们可能不是我们的重点研究对象。

在宋老师的书中提到：随机信号在时间上是无限的，因此是能量无限、功率有限的信号（我认为这意思就是随机信号就是功率信号）。而能量无限的信号不满足傅里叶变换绝对可积条件（狄利克雷条件的第三点，信号在一个周期内是绝对可积的。），因此随机信号的傅里叶变换是不存在的。但是随机信号的功率是有限的，采用功率谱可以从统计的角度来描述随机信号的频域特性。

这个我认为有必要说明一下，感觉很关键，不知道方法为啥来的，作用对象是啥，稀里糊涂的使用不是一名合格的工程师！

3.具体案例

周期图法绘制信号的功率谱，就是取一段有限长的信号，进行傅里叶变换然，然后取模的平方，再除以信号的长度。

这里有一个含噪信号y11,表达式如下：

这里，noise表示一个符合标准正态分布的随机噪声。采样频率设置为1000，信长度为2000，即信号的时间长度为2s。上述这个信号的时域波形如下图所示。

准确的来说，含噪信号是随机信号，不说它的幅值谱和相位谱了，但是如果认为该随机信号是各态历经，幅值谱其实和功率谱、功率谱密度提供的信息差不多，幅值谱也是进行了标准化，去除了信号长度对幅值谱的干扰，也是可以进行信号间信息比较的。但是，为了严谨，我就不放了，但是代码会画，大家自行取舍。

上图信号的功率谱和功率谱密度如下图所示：

无论是从上述功率谱和功率谱密度中，都能发现信号中存在两个突出的频率分量50Hz和300Hz。幅值数值的意义需要比较不同信号才有价值，我对数值没有直观感受，感觉不到。实际分析信号时，一般都只看功率谱密度。

注意：上述功率谱和功率谱密度画图时采用pow2db函数，将功率转换为了分贝数，关系式为：ydb=10log10(y)。周期图法中，功率谱和功率谱密度的区别在于，功率谱密度是功率谱再除以采样频率fs，可以理解为去除采用频率差异的影响。所以功率谱和功率谱密度的单位分别为dB和dB/Hz。

除此之外，值得注意的是该图中有两条曲线，一条蓝色的一条红色的，蓝色的为自己写的测试函数，红色的为matlab的periodogram函数，它们的结果是一模一样的，这也验证了测试函数的正确性。下图为计算的两个结果的偏差。具体如何计算两个结果的偏差大家可以看一下具体代码部分。

综上所述，该案例开展了采用周期图法含噪信号的功率谱和功率谱密度分析，测试了matlab函数和自写测试函数之间一致性（一致有前提请看05细节说明），验证了自写函数的正确性。值得说一句的是，这个自写函数不是我自己写的，是在一个老外的视频中看到的，链接如下：

https://www.bilibili.com/video/BV1ui421Z7XP/?spm_id_from=333.337.search-card.all.click&vd_source=50b116844fc6444ae380b4347979702f

4.具体代码

功率谱和功率谱密度，matlab提供了函数可以进行分析，我也在幅值谱和相位谱的基础上增加了功率谱和功率谱密度，并与matlab提供的周期图法periodogram函数结果进行比较。

我推荐matlab提供的周期图法periodogram函数，当然我进行验证的函数也提供，大家自行选择。

代码主要分为两个：

1、main1.m（主函数）

2、frequ_am_phase_power_psd.m（幅值谱、相位谱、功率谱和功率谱密度计算函数）

说明：

frequ_am_phase_power_psd.m函数

调用形式：

[freq,P1,Theta,P1_power,P1_psd]=frequ_am_phase_power_psd(y,fs,tol)

输入：

信号y，矩阵，行X列=单个信号的采样索引X信号数，比如信号的大小为8192X12，表示一个有12个信号的数据矩阵，每个信号长度为8192。注意，如果仅有一个信号，则y应该是一个列向量。同时，y的行数尽量为偶数，奇数的话会引起程序索引的警告。

 fs:采样频率，标量

 tol:相位阈值参数，标量（可缺省，默认10^-6）

输出

freq:表示幅值谱的横轴，向量，Hz

P1:表示幅值谱的纵轴，矩阵，单个信号的采样索引X信号数，类比信号y

Theta:表示相位谱的纵轴，矩阵，单个信号的采样索引X信号数，类比信号y

P1_power:表示功率谱的纵轴，矩阵，单个信号的采样索引X信号数，类比信号y

P1_psd:表示功率谱密度的纵轴，矩阵，单个信号的采样索引X信号数，类比信号y

main1.m（主函数）

%% 信号处理————频域分析（功率谱和功率谱密度）
%% 作者：冷漠
%% 时间：2024年6月16日
%% 关注公众号 ："故障诊断与寿命预测工具箱",每天进步一点点
clc
clear all
close all
%
fs=1000;                                             %采样频率
L=2000;                                              %信号长度
t=(0:L-1)/fs;                                          %时间序列
y1=0.8+0.7*cos(2*pi*50*t-pi/4)+0.8*cos(2*pi*300*t+pi/5);    %信号y1
y11=y1+0.5*randn(size(y1));                                 %加噪信号y11
​
​
y=y11';       %分析的信号，y的列为信号数量维度，y的行为每个信号的索引维度
figure('Position', [100, 100, 600, 100]);
plot(t,y11,'b');xlabel('时间(s)');ylabel('幅值');
​
[freq,P1,Theta,P1_power,P1_psd]=frequ_am_phase_power_psd(y,fs);   %幅值谱和相位谱     相位谱的阈值参数为1e-6
​
%画图
[P11_power P11_freq_power]=periodogram(y,rectwin(size(y,1)),size(y,1),fs,'power');
[P11_psd P11_freq_psd]=periodogram(y,rectwin(size(y,1)),size(y,1),fs,'psd');
​
​
%幅值谱
figure
stem(freq,P1(:,1),'k');    %作图
xlabel('频率(Hz)');ylabel('幅值');title('幅值谱');
%相位谱
figure
stem(freq,Theta(:,1),'k');    %作图
xlabel('频率(Hz)');ylabel('相位 (\pi)');title('相位谱');
%功率谱
figure
plot(freq,pow2db(P1_power(:,1)),'b');hold on;  %作图
plot(P11_freq_power,pow2db(P11_power(:,1)),'r');
xlabel('频率(Hz)');ylabel('功率(dB)');title('功率谱');
legend('测试函数','matlab函数');
fprintf('测试信号y的功率谱中测试函数与matlab函数之间的差为%f \n',sum(P1_power(:,1)-P11_power(:,1)));
%功率谱密度
figure
plot(freq,pow2db(P1_psd(:,1)),'b');hold on;  %作图
plot(P11_freq_psd,pow2db(P11_psd(:,1)),'r');
xlabel('频率(Hz)');ylabel('功率(dB/HZ)');title('功率谱密度');
legend('测试函数','matlab函数');
fprintf('测试信号y的功率谱密度中测试函数与matlab函数之间的差为%f \n',sum(P1_psd(:,1)-P11_psd(:,1)));

​function [freq,P1,Theta,P1_power,P1_psd]=frequ_am_phase_power_psd(y,fs,tol)（幅值谱、相位谱、功率谱和功率谱密度计算函数）

function [freq,P1,Theta,P1_power,P1_psd]=frequ_am_phase_power_psd(y,fs,tol)
%% 作者：冷漠
%% 时间：2024年6月01日
%% 关注公众号 ："故障诊断与寿命预测工具箱",每天进步一点点
%% 绘制信号频域的功率谱和功率谱密度
%% 参数解释： 
%     y: 表示输入信号，它可以为一个矩阵，行X列,具体为单个信号的采样索引X信号数
%        比如y的大小为8192X12，表示一个有12个信号的数据矩阵，每个信号长度为8192
%        注意，如果仅有一个信号，则y应该是一个列向量
%        同时，y的行数尽量为偶数，奇数的话会引起程序索引的警告
%     fs:表示采样频率
%     tol:相位阈值参数
%     freq:表示幅值谱的横轴
%     P1:表示幅值谱的纵轴
%     Theta:表示相位谱的纵轴
%     P1_power:表示功率谱的纵轴
%     P1_psd:表示功率谱密度的纵轴
​
if nargin==2
    tol=1e-6;  %计算误差的默认阈值
end
​
L=size(y,1);         % 信号长度
% Y=fft(y,2^nextpow2(L));          % FFT 快速傅里叶变换
Y=fft(y,L);          % FFT 快速傅里叶变换
freq=(0:L/2)*fs/L;   % 设置频率刻度  横轴Hz
%幅值谱
P2 = abs(Y/L);
P1 = P2(1:L/2+1,:);
P1(2:end-1,:) = 2*P1(2:end-1,:);  %纵轴 幅值
​
%相位谱
P2(2:end-1,:)=2*P2(2:end-1,:);
for i=1:size(Y,2)
    Y(P2(:,i)<tol,i) = 0;
    theta(:,i) = angle(Y(:,i))/pi;
end
Theta=theta(1:L/2+1,:);
​
%功率谱
P2_power=abs(fft(y,L)).^2/(L^2);   %双边功率谱
% P2_power=abs(fft(y,L)).^2/(L);   %双边功率谱
P1_power=P2_power(1:L/2+1,:);  %单边功率谱
P1_power(2:end-1,:)=2*P1_power(2:end-1,:);  %纵轴 功率
​
%功率谱密度
P2_psd=abs(fft(y,L)).^2/(L*fs);   %双边功率谱密度
P1_psd=P2_psd(1:L/2+1,:);  %单边功率谱密度
P1_psd(2:end-1,:)=2*P1_psd(2:end-1,:);  %纵轴 功率密度
​
end

5.细节说明

1、冷漠推荐大家采用matlab自带的函数（periodogram）进行功率谱和功率谱密度分析。

功率谱使用方法如下：

[P11_power P11_freq_power]=periodogram(y,rectwin(size(y,1)),size(y,1),fs,'power');

功率谱密度使用如下：

[P11_psd P11_freq_psd]=periodogram(y,rectwin(size(y,1)),size(y,1),fs,'psd');

在该方法中，一些参数设置如下：

输入：

y为所分析的信号，维度：采样点数X信号的数量，比如：y的大小为8192X12，表示一个有12个信号的数据矩阵，每个信号长度为8192。更多细节参考matlab关于该函数的官方文档。

rectwin(size(y,1))为窗函数，这个窗函数是矩形窗。

size(y,1)表示该信号的采样点数，即信号长度。

fs为采样频率

'power'表示功率谱，'psd'表示功率谱密度。

输出：

P11_power（P11_psd）：功率谱（功率谱密度）的纵轴。

P11_freq_power（P11_freq_psd）：功率谱（功率谱密度)的横轴。

2、自测函数与periodogram函数结果一致有三个前提（来自老外视频总结）：

1）采用的窗函数均为矩形窗函数。

2）截取的随机信号要是各态历经的平稳随机过程，这关系到功率谱的正确性。

3）进行离散信号的分析。

上述的第2)和3)在进行信号分析时就是成立的，只有当采用矩形窗函数时，自测函数和periodogram函数结果才一致。2)和3)不懂得推荐看老外视频。

3、周期图法求功率谱密度时，本文只进行了单次周期图法。换句话说，可以滑窗选择信号，进行平均。同时，改变窗函数，避免频谱泄露，如汉明窗等。这些内容，可参考（公众号：声振之家，文章名称：几种经典功率谱估计方法的实现（Matlab）及其局限性https://mp.weixin.qq.com/s/xvYS-YnUtlXntRsKXAkkMQ

4、功率谱中大家可以看一下自测函数我写得部分，这里有一个细节，就是我是取FFT变换后的模的平方，除以信号长度的平方。在其他科普中，大家经常看到的一般都是取FFT变换后的模的平方，除以信号的长度，这里的关键就是这个长度有没有平方。先说结论，如果没有平方和periodogram函数不一致，而且我认为应该有平方，理由如下：

大家可以计算一下余弦信号的功率，余弦信号的表达式为y=Acos（wt+a），它的功率为A^2/2，它的功率是一个余弦的平方积分，大家可以自行计算，我们在幅值谱中为了保证信号的幅值和表达式一致性，我们是取FFT变换后的模，除以信号的长度，即abs(FFT(y))/N=A/2（之所以除以2，有一个双边频谱向单边频谱的过渡，即双边频谱的幅值是单边的一半）,所以功率谱中A^2/4=（abs(FFT(y))/N）^2，这一平方就是N的平方，A^2/4向A^2/2过渡，就是再把双边频谱向单边频谱做个转换，所以老外说的对，这与matlab中的一致，佩服人家的严谨，有啥疑问欢迎大家评论。

6.总结

上述是一个随机信号的功率谱和功率谱密度的例子，采用的是周期图法，在原先封装了幅值谱和相位谱的函数基础上增加了功率谱和功率谱密度，方便大家使用，但是关于功率谱和功率谱密度还是推荐大家用matlab的periodogram函数。

这个例子参考了很多的例子，感谢巨人们的贡献。最后，尤其推荐大家看一下这个视频，能让你对幅值谱、功率谱和功率谱密度有个很好的理解。

https://www.bilibili.com/video/BV1ui421Z7XP/?spm_id_from=333.337.search-card.all.click&vd_source=50b116844fc6444ae380b4347979702f（感谢国内网友的搬砖！！！）

7.相关资料

附件

exl5

1、上述源码

     ①代码文件：

       main1.m(主函数)；

     frequ_am_phase_power_psd.m（幅值谱、相位谱、功率谱、功率谱密度测试函数）

2、相关参考

①书籍：MATLAB数字信号处理85个实用按比例——入门到进阶  宋知用 编著

②书籍：从这里学NVH噪声、振动、模态分析的入门与进阶   谭祥军 编著

③书籍：机械工程测试技术基础 第3版 熊诗波编著

④matlab官网的一些例子和教程：

a）频域分析实践介绍

https://ww2.mathworks.cn/help/signal/ug/practical-introduction-to-frequency-domain-analysis_zh_CN.html#d126e17836

b）使用 FFT 获得功率频谱密度估计

https://ww2.mathworks.cn/help/signal/ug/power-spectral-density-estimates-using-fft.html

c）频谱估计

https://ww2.mathworks.cn/help/signal/nonparametric-spectral-estimation.html

⑤了解功率谱密度和功率谱

https://www.bilibili.com/video/BV1ui421Z7XP/?spm_id_from=333.337.search-card.all.click&vd_source=50b116844fc6444ae380b4347979702f（一个涵盖了幅值谱、功率谱和功率谱密度的视频，强烈推荐！！！）

⑥、一些有价值的科普，理论说的更加全面

a）信号处理进阶之路，文章名称：能量域特征提取方法 https://mp.weixin.qq.com/s/W-qsi0e1aWC3pz3c4dy8kA

b）看海的城堡，文章名称：信号频域分析方法的理解（频谱、能量谱、功率谱和倒频谱）https://mp.weixin.qq.com/s/RpCoRAiZlo96wWYgW09lnw

c）声振之家，文章名称：几种经典功率谱估计方法的实现（Matlab）及其局限性https://mp.weixin.qq.com/s/xvYS-YnUtlXntRsKXAkkMQ

更多内容，请关注公众号“故障诊断与寿命预测工具箱”，每天进步一点点。