如何优雅地进行频谱分析—— 一行代码实现绘制MATLAB频谱、功率谱图_matlab频谱分析代码

之前的文章里讲了关于信号频谱、能量谱的相关理论和MATLAB编程实现方法：

Mr.看海：信号频域分析方法的理解（频谱、能量谱、功率谱、倒频谱、小波分析）

Mr.看海：频域特征值提取的MATLAB代码实现（频谱、功率谱、倒频谱）

不过对于刚刚开始进行这方面研究，或者对MATLAB不甚熟悉，抑或者不想花太多时间用在编程这件事上的同学，目前缺少一种足够便捷、高效、可靠的代码。

所以笔者想做这样一件事：编写一个函数，同学们只需要导入自己要分析的数据，并按照自己的需求设置部分选项，调用一行代码就可以实现绘制频谱图/功率谱图。

今天他来了。

一、频谱分析函数的基本用法

首先我们生成一段测试信号，它是由50Hz和120Hz的两段正弦信号叠加形成的：

%% 1.生成仿真信号
Fs = 1000;            % 采样频率          
T = 1/Fs;             % 采样间隔    
L = 1500;             % 信号长度
t = (0:L-1)*T;        % 时间轴
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
figure('color','w')
plot(t,S,'k');xlabel('时间');ylabel('幅值')

测试信号图像为：

我将这个频谱分析函数命名为pSpecAnalysis

现在调用一下这个函数：

pSpecAnalysis(S,Fs,[]);

此时就可以画出上边信号的频谱图：

频谱的主频、幅值都可以完美对上

所以在调用pSpecAnalysis时只需要输入时域信号和采样频率，足够简洁了吧。

如果想要分析你自己的信号，只需要替换函数的入口参数就可以。

二、频谱分析函数的扩展用法之——去趋势

但是如果仅仅实现画频谱图，还不足以专门写一篇文章。

细心地同学可能注意到了pSpecAnalysis函数还有第三个入口参数，在上边调用的时候给了空值。

这是我留下的结构体选项options参数，可以实现各种扩展功能。

第一个就是去趋势。

先说什么是去趋势。

不知道有没有同学在做频谱分析的时候，发现频谱图在频率为0时有一个很高的尖峰，就像这样：

这是由于信号中包含直流分量导致的，为了更好地提取各频段频谱特征，就需要从数据中减去均值或最佳拟合线，也就是去趋势。

如果想要对数据进行去趋势再画频谱，就要用到options结构体中的Detrend字段，调用函数的时候这样写：

options.Detrend = 1;
%    options.Detrend ：去趋势，设置options.Detrend = 1时，将在做fft前去除data趋势，
%                      当画出的频谱图在0频率时出现峰值，可将该选项置1，可以有效抑制0频率的峰值，默认情况不去趋势
[P1,f1] = pSpecAnalysis(S,Fs,options);

此时再画出来的频谱图就是这样的了：

二、频谱分析函数的扩展用法之——画图形式

这里对画图的形式也设置了可选项。

除了可以画上边的单张频谱图，还可以选择将原始信号和频谱图画在一起，就像这样：

此时的函数调用写法是（设置plotType字段）：

options.plotType = 2;
%    options.plotType ：画图种类选择，当options.plotType = 1时，只绘制信号频谱图（默认）
%                                    当options.plotType = 2时，绘制信号时域图和频谱图在同一张图上（上下两张）
[P1,f1] = pSpecAnalysis(S,Fs,options);

三、频谱分析函数的扩展用法之——绘制功率谱

功率谱也是很常用的，如果要画功率谱，像这样调用：

options.analysisType = 'psd';
%    options.analysisType ：频谱分析类型，当options.analysisType = 'fft'时，进行快速傅里叶变化FFT
%                                        当options.analysisType = 'psd'时，计算功率谱
[P1,f1] = pSpecAnalysis(S,Fs,options);

画出来的图是这样：

另外画功率谱图时还可以将纵坐标改为分贝（dB）表示，此时这样调用：

%    options.analysisType ：频谱分析类型，当options.analysisType = 'fft'时，进行快速傅里叶变化FFT
%                                        当options.analysisType = 'psd'时，计算功率谱
options.analysisType = 'psd';
 
%    options.psdUnit      ：功率谱单位，当options.psdUnit = 1时，单位为W/Hz(默认)
%                                      当options.psdUnit = 2时，单位为分贝dB
options.psdUnit = 1;
 
[P1,f1] = pSpecAnalysis(S,Fs,options);

画出的图像这样：

当然，在画功率谱图的时候，去趋势和画图样式也是同时可以设置的。

四、函数介绍

函数的主要功能介绍完了，整体看一下函数介绍：

function [p,f] = pSpecAnalysis(data,fs,options)
% 对数据进行fft变换
% 输入：
% data为待分析信号，需要是一维实数
% fs为采样频率，即每秒钟采样点数
% options为可选选项，使用结构体的方式导入，安装自身需求选择是否设置。具体包括：
%    options.Detrend ：去趋势，设置options.Detrend = 1时，将在做fft前去除data趋势，
%                      当画出的频谱图在0频率时出现峰值，可将该选项置1，可以有效抑制0频率的峰值，默认情况不去趋势
%    options.analysisType ：频谱分析类型，当options.analysisType = 'fft'时，进行快速傅里叶变化FFT
%                                        当options.analysisType = 'psd'时，计算功率谱
%    options.psdUnit      ：功率谱单位，当options.psdUnit = 1时，单位为W/Hz(默认)
%                                      当options.psdUnit = 2时，单位为分贝dB
%    options.plotType ：画图种类选择，当options.plotType = 1时，只绘制信号频谱图（默认）
%                                    当options.plotType = 2时，绘制信号时域图和频谱图在同一张图上（上下两张）

函数预留了扩展接口，如果同学们还想要加入什么功能可以留言告诉我。

上边的测试代码和pSpecAnalysis函数，可以在下边链接获取：

频谱/功率谱快速画图程序 | 工具箱文档

编程不易，感谢支持~

最后

后边可能会不定期地将一些常用方法封装一下，作为专栏的番外了。

类EMD分解方法介绍完之后，将开始重点讲时间序列分析和深度学习。

敬请期待！