本文介绍信号频域分析中幅值谱和相位谱的相关代码和分析过程！
这个内容相关博主可能也分享过，但是着急想用找不到一些参考代码，感觉不方便，所以还是想梳理一下这块的内容。
该内容参考了一些资料：
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)');
endfrequ_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个实用按比例——入门到进阶 宋知用 编著）
更多内容，请关注公众号“故障诊断与寿命预测工具箱”，每天进步一点点。
扫一扫
1108
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
