clc;clear
%% 导入数据
load(‘ct.mat’)
Fs=12800;
a=[];c=[];w=[];
%% 1时域特征提取
for i=2:8
y=ct(:,i);
a(1,i) = max(y); %最大值
a(2,i)= min(y); %最小值
a(3,i) = mean(y); %平均值
ma=a(1,i) ;mi=a(2,i);
a(4,i) = ma-mi; %峰-峰值
a(5,i) = mean(abs(y)); %绝对值的平均值(整流平均值)
a(6,i) = var(y); %方差
a(7,i) = std(y); %标准差
a(8,i) = kurtosis(y); %峭度
a(9,i)= skewness(y); %偏度
a(10,i) = rms(y); %均方根
av=a(5,i) ;rm=a(10,i) ;pk=a(4,i);
a(11,i) = rm/av; %波形因子
a(12,i) = pk/rm; %峰值因子
a(13,i) = pk/av; %脉冲因子
a(14,i) = mean(sqrt(abs(y))).^2;
xr=a(14,i);
a(15,i) = pk/xr; %裕度因子
%% 2频域特征提取
[ frequencystruct ]=frequencyDomainFeatures( y,Fs);
c(1,i)=frequencystruct.MF;%频率均值
c(2,i)=frequencystruct.FC;%重心频率
c(3,i)=frequencystruct.RMSF;%频率均方根
c(4,i)=frequencystruct.RVF;%频率标准差
%% 3小波特征提取
[ waveletstruct ] = waveletFeatures(y);
w(1,i)=waveletstruct.E1;
w(2,i)=waveletstruct.E2;
w(3,i)=waveletstruct.E3;
w(4,i)=waveletstruct.E4;
w(5,i)=waveletstruct.E5;
w(6,i)=waveletstruct.E6;
w(7,i)=waveletstruct.E7;
w(8,i)=waveletstruct.E8;
w(9,i)=waveletstruct.p1;
w(10,i)=waveletstruct.p2;
w(11,i)=waveletstruct.p3;
w(12,i)=waveletstruct.p4;
w(13,i)=waveletstruct.p5;
w(14,i)=waveletstruct.p6;
w(15,i)=waveletstruct.p7;
w(16,i)=waveletstruct.p8;
w(17,i)=waveletstruct.qyshang;
end
a=a’;c=c’;w=w’;

function [ frequencystruct ] = frequencyDomainFeatures( src,fs)
%计算频域统计特征
%对信号进行FFT变换********
FS=fs;
N=length(src);n=0:N-1;
freq=n*fs/N;
f=abs(fft(src,N)*2/N);

x=f(1:N/2); %纵坐标 频率幅值
freq=freq(1:N/2)’; %横坐标 频率值

% plot(freq,x);
% title(‘原始信号频域波形’);
% xlabel(‘频率/hz’);
% ylabel(‘幅值/v’);

%***计算频域特征值
frequencystruct.MF=mean(x); %平均频率
frequencystruct.FC=sum(freq.*x)/sum(x);%重心频率
frequencystruct.RMSF=sqrt(sum([freq.^2].*x)/sum(x));%频率均方根
frequencystruct.RVF=sqrt(sum([(freq-frequencystruct.FC).^2].*x)/sum(x));%频率标准差

%第一级边频带幅值和=上下最靠近啮合频率的边频带对应序号的频率幅值 这里随便写
%已知：行星轮个数、转频，可以找到啮合频率，然后得到边频带，即可计算

% frequencystruct.FSB=x(1000)+x(2000);%第一级边频带幅值和
% frequencystruct.FSI=frequencystruct.FSB/2;%边频带指数
% frequencystruct.FM0=(max(x)-min(x))/frequencystruct.FSB;%FM0
% frequencystruct.FSLF=frequencystruct.FSB/std(x,1);%边频带等级因子

end

function [ waveletstruct ] = waveletFeatures( src)
%计算小波域统计特征
%wpcoef：找小波包分解的系数，N层小波包分解后将信号分成2^{N个子带，没个子带的信号长度是原是信号的1/2}N倍，其实就是将原始信号化成2^N段，每段的长度是相等的且比原信号短
%wprcoef：是小波包分解系数重构，将分解的各个子带的信号拓展为原始信号那么长，
%两者的区别就是wprcoef是wpcoef的拓展（长度与原始信号一样），他们每个节点的能量是一样的

%%小波熵是一个统称，不是一种具体的熵算法。
%旗下包括小波时间熵、小波尺度熵(小波特征尺度熵)、小波时频熵、小波能量熵(小波能谱熵)
% 小波相关尺度熵(小波相关特征尺度熵)、小波奇异熵、小波方差谱熵、小波空间特征熵(小波空间特征尺度熵)

n=3;%3层小波包分解
T=wpdec(src,n,‘db6’); %3层小波包分解 得到8个子带
for i=1:2^n %%重构分量，使信号和原始信号一样长
E(i)=norm(wprcoef(T,[n,i-1]),2)*norm(wprcoef(T,[n,i-1]),2); %求第i个节点的二范数，乘以平方才是能量
end

E_total=sum(E);%总能量

%8个子带的小波能量比
waveletstruct.p1=E(1)/E_total;
waveletstruct.p2=E(2)/E_total;
waveletstruct.p3=E(3)/E_total;
waveletstruct.p4=E(4)/E_total;
waveletstruct.p5=E(5)/E_total;
waveletstruct.p6=E(6)/E_total;
waveletstruct.p7=E(7)/E_total;
waveletstruct.p8=E(8)/E_total;

%绘制条形图
% p=[waveletstruct.p1,waveletstruct.p2,waveletstruct.p3,waveletstruct.p4,…
% waveletstruct.p5,waveletstruct.p6,waveletstruct.p7,waveletstruct.p8];
% bar§;
% set(gca,‘XTickLabel’,{’(3,0)’,’(3,1)’,’(3,2)’,’(3,3)’ ,’(3,4)’,’(3,5)’,’(3,6)’,’(3,7)’});

%小波能量熵(小波能谱熵)
waveletstruct.energyE=-sum(E/E_total.*log(E/E_total));

%8个子带的小波熵(小波尺度熵) 在各尺度上求的是shannon熵=小波尺度熵(小波特征尺度熵)
waveletstruct.E1=wentropy(wpcoef(T,[3,0]),‘shannon’);%
waveletstruct.E2=wentropy(wpcoef(T,[3,1]),‘shannon’);%
waveletstruct.E3=wentropy(wpcoef(T,[3,2]),‘shannon’);%
waveletstruct.E4=wentropy(wpcoef(T,[3,3]),‘shannon’);%
waveletstruct.E5=wentropy(wpcoef(T,[3,4]),‘shannon’);%
waveletstruct.E6=wentropy(wpcoef(T,[3,5]),‘shannon’);%
waveletstruct.E7=wentropy(wpcoef(T,[3,6]),‘shannon’);%
waveletstruct.E8=wentropy(wpcoef(T,[3,7]),‘shannon’);%

%小波奇异谱熵
%%先求取奇异值，将奇异值构造向量即得到奇异值谱，然后计算熵值即可
%第一步：小波分解，得到各分量信号，对系数重构，得到完整系数
for i=1:2^n %%重构分量，使信号和原始信号一样长
q(:,i)=wprcoef(T,[n,i-1]); %求系数
end

%奇异值分解
s=svd(q);
S=sum(s);
dim=length(s);
for i=1:dim
p(i)=s(i)/S;
end
waveletstruct.qyshang=-sum(p.*log§);%小波奇异熵

end