信号特征提取

In machine learning, pattern recognition and in image processing, feature extraction starts from an initial set of measured data and builds derived values (features) intended to be informative and non-redundant, facilitating the subsequent learning and generalization steps, and in some cases leading to better human interpretations. Feature extraction is related to dimensionality reduction.

When the input data to an algorithm is too large to be processed and it is suspected to be redundant (e.g. the same measurement in both feet and meters, or the repetitiveness of images presented as pixels), then it can be transformed into a reduced set of features (also named a feature vector). Determining a subset of the initial features is called feature selection.^[1] The selected features are expected to contain the relevant information from the input data, so that the desired task can be performed by using this reduced representation instead of the complete initial data.

I. 原始信号

x=x1;
fs=15360;
N=length(x);
t=(0:N-1)/fs;

figure
plot(t,x);
title('sig1波形图');
xlabel('时间 ');
ylabel('振幅 ');

II. 时域信号

N=length(x);
p1=mean(x); %均值 1/N*sum(x)
p2=sqrt(sum((x-p1).^2)/N); %均方根值
p3=(sum(sqrt(abs((x-p1))))/N).^2; %方根幅值
p4=mean(abs((x-p1))); %绝对平均值
p5=sum((x-p1).^3)/N; %偏斜度
p6=sum((x-p1).^4)/N; %峭度
p7=sum((x-p1).^2)/N; %方差
p8=max(abs((x-p1)));%峰值 最大值
p9=min((x-p1));%最小值
p10=p8-p9;%峰峰值

% val返回有量纲指标
val=[p1; p2; p3; p4; p5; p6; p7; p8; p9; p10];

%以上都是有量纲统计量，以下是无量纲统计量
f1=p2./p4; %波形指标 f1=p2/(sum(x)/N);
f2=p8./p2; %峰值指标 E[MAX(X)]=P8? 峰值/有效值
f3=p8./p4; %脉冲指标
f4=p8./p3; %裕度指标
f5=p5./(p2.^3); %偏斜度指标
f6=p6./(p2.^4); %峭度指标,峭度/标准差四次方
% f6=kurtosis(x);%峭度指标

% factor返回无量纲指标
factor=[f1; f2; f3; f4; f5; f6];

valfactor=[val;factor];

III. 频域信号

fre_line_num=max(size(y));
p1 = mean(y);
% 均值频率 特征1,反映频域振动能量的大小；
p2 = sum((y-p1).^2)/fre_line_num;
% 标准差 特征2，表示频谱的分散或者集中程度；
p3 = sum((y-p1).^3)/(fre_line_num*sqrt(p2^3));
% 特征3，表示频谱的分散或者集中程度；
p4 = sum((y-p1).^4)/(fre_line_num*p2^2);
% 特征4，表示频谱的分散或者集中程度；
meanf = sum(f.*y)/sum(y);
p5 = meanf;
% 频率中心 特征5，反映主频带位置的变化；
sigma = sqrt(sum((f-meanf).^2.*y)/fre_line_num);
p6 = sigma;
% 特征6，表示频谱的分散或者集中程度；
p7 = sqrt(sum(f.^2.*y)/sum(y));
% 均方根频率 特征7，反映主频带位置的变化；
p8 = sqrt(sum(f.^4.*y)/sum(f.^2.*y));
% 特征8，反映主频带位置的变化；
p9 = sum(f.^2.*y)/sqrt(sum(y)*sum(f.^4.*y));
% 特征9，反映主频带位置的变化；
p10 = sigma/meanf;
% 特征10，表示频谱的分散或者集中程度；
p11 = sum((f-meanf).^3.*y)/(sigma.^3*fre_line_num);
% 特征11，表示频谱的分散或者集中程度；
p12 = sum((f-meanf).^4.*y)/(sigma.^4*fre_line_num);
% 特征12，表示频谱的分散或者集中程度；
p13 = sum(sqrt(abs(f-meanf)).*y)/(sqrt(sigma)*fre_line_num);
% 特征13，表示频谱的分散或者集中程度；
fac=[p1;p2;p3;p4;p5;p6;p7;p8;p9;p10;p11;p12;p13];

IV. 3维频谱-FFT

load V200.txt;
load V250.txt;
load V300.txt;
load V400.txt;

sample_point=6144;%采样点数为6144点
filenum=4;
yumatrix=[V200,V250,V300,V400];
rotate_speed=[406 531 719 906];

%求频谱
for k=1:1:filenum
    temp_value(:,k)=fft(yumatrix(:,k))/sample_point;
    for l=1:1:sample_point/2
        frequency_value(l,k)=2*abs(temp_value(l,k));
    end
end

for k=1:1:filenum
    sample_frequency(k)=3200;%计算采样频率
    f(k)=sample_frequency(k)/sample_point;%画图的频率
    beipin(k)=rotate_speed(k)/60;
    b(k)=f(k)/beipin(k);
    for l=1:1:sample_point/2
        z(l,k)=rotate_speed(k);
    end
end
for m=1:1:filenum
    for n=1:1:200
        final_f(n,m)=n*b(m);
        final_value(n,m)=frequency_value(n,m);
        final_z(n,m)=z(n,m);
    end
end

figure
plot3(final_f,final_z,final_value);
xlabel('倍频');
ylabel('转速（rpm)');
zlabel('幅值（um)');
grid on;

V. hilbert包络

% hilbert包络
y_hilbert_1=hilbert(x_input);
y_hilbert_2=abs(y_hilbert_1); % 包络信号

VI. FFT频谱/解调

N=length(x_input);
fs=15360;
f=(0:N/2-1)*fs/N; % f=(0:N-1)*fs/N等价，为什么？

% 功率谱,解调频谱图
p=abs(fft(x_input-mean(x_input)),N)*2/N;
% p=abs(fft(x_input),N)*2/N;
% p=yfft.*conj(yfft)/N;

figure
plot(f,p(l:N/2),'k-');
xlabel('频率 f/Hz');
ylabel('幅值 P/W'); % (m/s2)

VII. 小波系数分解与重构

m=6;
[c_11,l_11]=wavedec(y_11,m,'db6');

% 对第六层低频系数进行重构——
a_11=wrcoef('a',c_11,l_11,'db6',6);

% 对分解的1-m层的高频系数（细节信号）进行重构
    for i=1:m
        d_11(m+1-i)=wrcoef('d',c_11,l_11,'db6',m+1-i);

        subplot(m+2,1,i+2)
        plot(d_11(m+1-i))
        ylabel(['d',num2str(m+1-i)])
    end

% 分解
wpt=wpdec(xdatac,3,'db1','shannon');
% 重构
s330=wprcoef(wpt,[3,0]);
s331=wprcoef(wpt,[3,1]);
s332=wprcoef(wpt,[3,2]);
s333=wprcoef(wpt,[3,3]);
s334=wprcoef(wpt,[3,4]);
s335=wprcoef(wpt,[3,5]);
s336=wprcoef(wpt,[3,6]);
s337=wprcoef(wpt,[3,7]);

% 小波包分解，重构轴承振动信号，Hilbert包络，FFT进行频谱分析
%采样频率
fs=12000;
%计算其自相关序列
xdatac=xcorr(xdata);
[t,d]=wpdec(xdatac,3,'db5');% Wavelet packet decomposition 1-D

%求每组系数>>>>>阈值=均值+标准方差
thr0=mean([3,0])+std([3,0]);
thr1=mean([3,1])+std([3,1]);
thr2=mean([3,2])+std([3,2]);
thr3=mean([3,3])+std([3,3]);
thr4=mean([3,4])+std([3,4]);
thr5=mean([3,5])+std([3,5]);
thr6=mean([3,6])+std([3,6]);
thr7=mean([3,7])+std([3,7]);
%THR=[thr0,thr1,thr2,thr3,thr4,thr5,thr6,thr7];

%消噪
xd0=wpdencmp([3,0],'s',3,'db5','threshold',thr0,0);
xd1=wpdencmp([3,1],'s',3,'db5','threshold',thr1,0);
xd2=wpdencmp([3,2],'s',3,'db5','threshold',thr2,0);
xd3=wpdencmp([3,3],'s',3,'db5','threshold',thr3,0);
xd4=wpdencmp([3,4],'s',3,'db5','threshold',thr4,0);
xd5=wpdencmp([3,5],'s',3,'db5','threshold',thr5,0);
xd6=wpdencmp([3,6],'s',3,'db5','threshold',thr6,0);
xd7=wpdencmp([3,7],'s',3,'db5','threshold',thr7,0);
%XD=wpdencmp([t,d],'s',3,'db5','threshold',THR,0);%用一个公式代替

%计算每组系数的能量
p0=sum(xd0^2);
p1=sum(xd1^2);
p2=sum(xd2^2);
p3=sum(xd3^2);
p4=sum(xd4^2);
p5=sum(xd5^2);
p6=sum(xd6^2);
p7=sum(xd7^2);
%找到最大的能量所对应的序数
pmax=max(p0,p1,p2,p3,p4,p5,p6,p7);
if pmax==p0 xdmax=xd0;
elseif pmax==p1 xdmax=xd1;
elseif pmax==p2 xdmax=xd2;
elseif pmax==p3 xdmax=xd3;
elseif pmax==p4 xdmax=xd4;
elseif pmax==p5 xdmax=xd5;
elseif pmax==p6 xdmax=xd6;
else xdmax=xd7;
end
%计算其自相关序列
rce=xcorr(xdmax);
%hilbert包络&FFT频谱分析

---