用matlab计算时序数据的时域、频域特征，附代码_时间序列的频域特征

用matlab计算时序数据的时域、频域特征

主思路都是利用滑动窗对数据先进行分割，对每一段分别计算特征。
废话少说，上代码：

%对于输入的一个序列片段，计算以下时域统计特征，以 行 的形式返回特征值
function T=regular_statitic_features(data)
% data为输入的振动信号，rpm为对应的转速信号
%计算信号的时域特征
ma = max(data); 			%最大值
mi = min(data); 			%最小值	
me = mean(data); 			%平均值
pk = ma-mi;			%峰-峰值
av=mean(abs(data));     %整流平均值
r=rms(data);           %有效值
p=peak(data);          %峰值
va = var(data);		%方差
st = std(data);		%标准差
ku = kurtosis(data);   %峭度
sk = skewness(data);   %偏度
rm = rms(data);        %均方根
S = rm/av;			%波形因子
C = pk/rm;			%峰值因子
I = pk/av;			%脉冲因子
YD=marginfactor(data); %裕度因子

% 傅里叶变化获得频谱
[f,result_FFt]=transToFFT(data,20000);
% 计算频域特征
F1=fc(f,result_FFt);   %重心频率
F2=msf(f,result_FFt);  %均方频率
F3=vf(f,result_FFt);   %频率方差
F4=BandEnergy(f,result_FFt); % 计算350-700Hz频带能量
F5=RPSD(result_FFt); % 相对功率谱熵

T=[ma mi me pk av r p va st ku sk rm S C I YD F1 F2 F3 F4 F5];
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

以下是一些特征的具体计算过程：

%重心频率
function F1=fc(f,result_FFt)
    F1=sum(f.*result_FFt)/(sum(result_FFt));
end

%均方频率
function F2=msf(f,result_FFt)
    f=f.^2;
    F2=sum(f.*result_FFt)/(sum(result_FFt));
end

%频率方差
function F3=vf(f,result_FFt)
    f=(f-fc(f,result_FFt)).^2;
    F3=sum(f.*result_FFt)/(sum(result_FFt));
end

%频带能量
function F4=BandEnergy(f,result_FFt)
    startIndex=sum(f<=350);
    endIndex=sum(f<700);
    F4=sum(result_FFt(startIndex:endIndex));
end

%相对功率谱熵
function F5=RPSD(result_FFt)
    Fi=result_FFt.*result_FFt;
    Pi=Fi/sum(Fi);
    F5=-sum(Pi.*log2(Pi))/log2(length(result_FFt));
end

function pk=peak(a)
%求信号的峰值
pk=max(a)-min(a);
end

function crestfactor=crestfactor(a)
%求信号的峰值因子
crestfactor=peak(a)/rms(a);
end

function shapefactor=shapefactor(a)
%求信号的波形因子
shapefactor=rms(a)/mean(abs(a));
end

function impulsefactor=impulsefactor(a)
%求信号的脉冲因子
impulsefactor=peak(a)/mean(abs(a));
end

function marginfactor=marginfactor(a)
%求信号的裕度因子
marginfactor=peak(a)/mean(sqrt(abs(a)))^2;
end

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

其中用到了计算频谱的函数transToFFT，代码如下：

function [f,result_FFt]=transToFFT(data,fs)
%用于求解数据的快速傅里叶变换结果，以便于绘制频谱图
N=length(data);
%去均值
data=data-mean(data);
%求频率分辨率
df=fs/(N-1);
f=(0:N-1)*df;
Y=fft(data)/N*2;
result_FFt=abs(Y);
%取一半
result_FFt=result_FFt(1:ceil(N/2));
f=f(1:ceil(N/2));
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14

用的时候放到一个路径下就好啦，另外数据分段的代码如下：

%输入一段长序列，利用滑动窗的思路切割成数据片段，返回以 行 为片段的矩阵
function data_segments=enframe(data,win,inc)
%win:窗长；inc：移贞距离
data=data';   %根据数据是行或列决定
data_segments=[];
left=1;
right=win;
while right<=length(data)
    seg=data(left:right);
    data_segments=[data_segments;seg];
    left=left+inc;
    right=right+inc;
end
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14

之前还写了两个文章，一个用python，一个用的matlab，但都不是很详尽，且可能有些错误，因此重新写了一篇，可以结合着看。另外两篇里面有几个参考文献，大家可以理解记忆。

python代码：https://blog.csdn.net/weixin_44620044/article/details/107239960?spm=1001.2014.3001.5502

matlab代码：https://blog.csdn.net/weixin_44620044/article/details/105617749?spm=1001.2014.3001.5502
转载请注明出处，谢谢！