信号处理—频域分析（幅值谱、相位谱）

本文介绍信号频域分析中幅值谱和相位谱的相关代码和分析过程！

这个内容相关博主可能也分享过，但是着急想用找不到一些参考代码，感觉不方便，所以还是想梳理一下这块的内容。

该内容参考了一些资料：

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

2、matlab官网的一些例子：

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

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

前言

说到信号处理，绕不开的话题就是信号的频域分析，频域分析的基础是傅里叶变换（FFT），它认为信号可以由正弦信号和余弦信号组成，而不论正弦信号还是余弦信号，它的两个重要参数分别是频率和初始相位。

在信号的频域分析中，幅值谱和相位谱是平稳信号分析两个手段，它能直观地反映信号的频域特性。

具体案例

这里有四个信号y1、y11、y2、y22，他们的表达式如下。

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

图1

图中，从上到下依次对应信号y1、y11、y2、y22，从图中能发现，在噪声的干扰下，y11和y22的波形发生了变形。

下面利用绘制上述这几个信号频域的幅值谱和相位谱如下。

y1信号的幅值谱（左）和相位谱（右）：

图2

从y1的幅值谱中能找到幅值为0.8的直流分量、频率为50Hz幅值为0.7的余弦分量、频率为300Hz幅值为0.8的余弦分量。

从y1相位谱中能找到频率为50Hz的余弦分量的初始相位为-0.25pi，频率为300Hz的余弦分量的初始相位为0.2pi，直流分量的初始相位为0（好奇请往下看）。（注意相位谱的单位为pi）

这验证了幅值谱和相位谱代码的准确性。

y11信号的幅值谱（左）和相位谱（右）

图3

y11是在y1的基础上增加了噪声，能明显从其幅值谱中看到幅值较低的随机噪声分量，这些噪声的初始相位是杂乱无章的。

y2信号的幅值谱（左）和相位谱（右）

图4

从y2的幅值谱中能找到幅值为0.8的直流分量、频率为120Hz幅值为0.7的余弦分量、频率为340Hz幅值为0.8的余弦分量。

从y2相位谱中能找到频率为50Hz的余弦分量的初始相位为-5pi/6，频率为340Hz的余弦分量的初始相位为-0.2pi，直流分量的初始相位为0（好奇请往下看）。（注意相位谱的单位为pi）

这验证了幅值谱和相位谱代码的准确性。

y22信号的幅值谱（左）和相位谱（右）

图5

y22是在y2的基础上增加了噪声，能明显从其幅值谱中看到幅值较低的随机噪声分量，这些噪声的初始相位是杂乱无章的。

综上所述，y1和y2是为了验证幅值谱和相位谱代码的准确性。

关于相位谱需要注意，默认信号是余弦信号求其初始相位，所以0.8直流分量=0.8cos（0），所以其初始相位是0，y2中0.7sin(240pi*t-pi/3)=0.7cos(240pi*t-5pi/6),所以图中展示的值为-5/6=-0.83333。

具体代码

代码主要分为两个：

1、main.m（主函数）

2、frequ_am_phase.m（幅值谱和相位谱计算函数）

说明：

frequ_am_phase.m函数

调用形式：

[freq,P1,Theta]=frequ_am_phase(y,fs,tol)

输入：

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

     fs:采样频率，标量

     tol:相位阈值参数，标量（可缺省，下面讲解）

输出

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

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

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

main.m（主函数）

%% 信号处理————频域分析（幅值谱和相位谱）
%% 作者：冷漠
%% 时间：2024年6月02日
%% 关注公众号 ："故障诊断与寿命预测工具箱",每天进步一点点
clc
clear all
close all
%
fs=1000;                                             %采样频率
L=1000;                                              %信号长度
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
y2=0.8+0.7*sin(2*pi*120*t-pi/3)+0.8*cos(2*pi*340*t-pi/5);   %信号y2
y22=y2+0.5*randn(size(y2));                                 %加噪信号y22

y=[y1;y11;y2;y22]';       %所有的信号，y的列为信号数量维度，y的行为每个信号的索引维度
figure;
subplot(411);plot(t,y1,'b');xlabel('时间(s)');ylabel('幅值');
subplot(412);plot(t,y11,'b');xlabel('时间(s)');ylabel('幅值');
subplot(413);plot(t,y2,'b');xlabel('时间(s)');ylabel('幅值');
subplot(414);plot(t,y22,'b');xlabel('时间(s)');ylabel('幅值');

[freq,P1,Theta]=frequ_am_phase(y,fs);   %幅值谱和相位谱     相位谱的阈值参数为1e-6
% [freq,P1,Theta]=frequ_am_phase(y,fs,0.1);   %幅值谱和相位谱 相位谱的阈值参数为0.1
% [freq,P1,Theta]=frequ_am_phase(y,fs,0);   %幅值谱和相位谱   相位谱的阈值参数为0
%画图
for i=1:size(P1,2)
    figure
    stem(freq,P1(:,i),'k');    %作图
    xlabel('频率(Hz)');ylabel('幅值');
    figure
    stem(freq,Theta(:,i),'k');    %作图
    xlabel('频率(Hz)');ylabel('相位 (\pi)');
end

frequ_am_phase.m（幅值谱和相位谱计算函数）
 

function [freq,P1,Theta]=frequ_am_phase(y,fs,tol)
%% 作者：冷漠
%% 时间：2024年6月01日
%% 关注公众号 ："故障诊断与寿命预测工具箱",每天进步一点点
%% 绘制信号频域的幅值谱和相位谱
%% 参数解释： 
%     y: 表示输入信号，它可以为一个矩阵，行X列,具体为单个信号的采样索引X信号数
%        比如y的大小为8192X12，表示一个有12个信号的数据矩阵，每个信号长度为8192
%        注意，如果仅有一个信号，则y应该是一个列向量
%        同时，y的行数尽量为偶数，奇数的话会引起程序索引的警告
%     fs:表示采样频率
%     tol:相位阈值参数
%     freq:表示幅值谱的横轴
%     P1:表示幅值谱的纵轴
%     Theta:表示相位谱的纵轴

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,:);
end

细节说明

1、输入信号的长度为偶数，奇数会引起frequ_am_phase函数中索引警告。

2、相位阈值参数tol默认选择为1e-6，它直接影响相位谱的质量。下图为主函数中执行：

[freq,P1,Theta]=frequ_am_phase(y,fs,0);   

    将[freq,P1,Theta]=frequ_am_phase(y,fs);注释掉，结果如下：

y1的幅值谱和相位谱

图6

y2的幅值谱和相位谱

图7

对比缺省参数下相位阈值参数为1e-6，即图2和图4,当相位阈值参数为0时，纯净的信号y1和y2的相位谱是杂乱无章的。这是因为matlab进行FFT变换后产生了计算误差。

借用（书籍：MATLAB数字信号处理85个实用按比例——入门到进阶  宋知用 编著）解释，取y1的FFT变换后的前5个分量观察其实部和虚部，结果如下：

其中，y1（1,1）为直流分量0.8对应的FFT变换结果（800/L，L为信号长度1000）。

初始相位求解的思路是把上述某个值，比如8e2+0i放到一个复数坐标系中，计算该点与X轴（实数轴）正向的夹角，夹角范围为[-pi,pi]，所以8e2+0i在实数正半轴上，所以初始相位是0，y1(51,1)为50Hz处的FFT结果为2.4749e2-2.4749e2i，它落在了第四象限45度分界线上，所以它的初始相位是-pi/4。

在上述这些y1（2,1）=-1.1437e-15-1.15e-13i，

y1（3,1）=-3.4681e-13+1.9273e-13i，

y1（4,1）=5.9391e-13-3.9776e-14i，

y1（5,1）=-8.1909e-13+5.586e-13i，

它们的数量级实部和虚部在1e-14至1e-13量级，计算其模达到了1e-12，几乎为0，然而在matlab中采用反三角函数计算初始相位时实部和虚部量级会抵消，导致最终初始相位分布在-pi到pi之间，从而造成了相位混乱。针对这种现象可以增加一个阈值参数来避免，如本代码中的tol。

本文采用的tol默认为1e-6,通过判断幅值谱上幅值小于1e-6的位置，将其置为0，从而达到避免相位混乱的现象。

为了能方便地获取频域幅值谱中意向谱峰的初始相位，tol可人为设定，如以这个案例为例，在主函数中执行：

[freq,P1,Theta]=frequ_am_phase(y,fs,0.1);   

设置相位阈值参数tol为0.1，那么幅值谱中幅值低于0.1的位置，其FFT变换都会置为0，这基本能够滤掉人为添加的噪声，因为我添加的正态分布的噪声其幅值从幅值谱中能发现，幅值低于0.1，那么针对含噪信号y11和y22也能获得干净的相位谱，结果如下：

tol为0.1时，y11的幅值谱和相位谱：

tol为0.1时，y22的幅值谱和相位谱：

总结

上述是一个测试信号频域幅值谱和相位谱分析的例子，封装了幅值谱和相位谱的函数，方便大家使用，其中一些细节也没有办法面面俱到，使用时参考注意细节。

幅值谱采用了缩放因子（信号的长度），并且将双边频谱转换为单边频谱，因为负频率没有物理含义，同时正频率部分乘2，这样能保证幅值谱和表达式中幅值一一对应。

相位谱设置了相位阈值参数，它对于提取幅值谱中意向谱峰的初始相位是有用。

上述这个例子主要参考了matlab的官方案例子（https://ww2.mathworks.cn/help/signal/ug/practical-introduction-to-frequency-domain-analysis_zh_CN.html#d126e17836）和书籍（MATLAB数字信号处理85个实用按比例——入门到进阶  宋知用 编著）

相关资料

附件

55qd

1、上述源码

     ①代码文件：

       maintf.m(主函数)；

      frequ_am_phase.m（幅值谱和相位谱计算函数）

2、相关参考

    ①https://ww2.mathworks.cn/help/signal/ug/practical-introduction-to-frequency-domain-analysis_zh_CN.html#d126e17836（Mathwork的一个频域分析例子）

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

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