newpdprocess_backup.m

神经网络去噪 matlab程序 神经网络去除随机脉冲干扰

%newpdprocess
%局放的随机脉冲干扰去除程序
%调用方式：pddataout=newpdprocess(pddatain,flag_fir)
%pddatain为要处理的数据
%pddataout为处理后输出的数据
%flag_fir为fir滤波标志位，为1表示pddatain已经经过滤波，否则则没
function pddataout=newpdprocess(pddatain)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                  FIR滤波[500k 10m]                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

data_fir=pdfilter(pddatain);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           提取脉冲                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

bias=0.002;%脉冲的提前阈值为0.002v。
[pulse,pulse_place]=pddistill(pddatain);
pulse_character=chara(pulse);%计算脉冲的特征值
load pddata.mat chara_max;
pulse_character=charaguiyi(pulse_character,chara_max);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           脉冲判断                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load pddata.mat pdnet;%读取神经网络文件
output_sim=netsim(pdnet,pulse_character);
jam_place=pulse_place(:,find(output_sim<=0.5));%找出随机脉冲干扰位置

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           去除干扰                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pddataout=data_fir;
[r,q]=size(jam_place);
for i=1:q
    pddataout(jam_place(1,i):jam_place(2,i))=0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           绘图                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure;
subplot(2,1,1);
plot(pddatain);
title('Data before processing');
subplot(2,1,2);
plot(pddataout);
title('Data (net) after processing');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        FIR滤波子程序                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dataout=pdfilter(datain)
W1=500*1e3;%下截止频率
W2=10*1e6;%上截止频率
fs=25*1e6;%采样频率
Wn=[W1/(fs/2) W2/(fs/2)];%带通[500k 10m]
pdfir=fir1(150,Wn,'bandpass');
dataout=filter(pdfir,1,datain);%进行（500k 10m)的fir滤波


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        脉冲提取子程序                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%根据包络线来提取脉冲波形
%脉冲提取的原则：1、包络线边界的单调性开始改变
%                2、原信号数据点不低于阈值
%function [pulse,pulse_place]=pddistill(handledata,bias,cutoffband)
function [pulse,pulse_place]=pddistill(handledata,bias,cutoffband)
%pulse为提取出来的脉冲数组
%pulse_place为脉冲数组的位置坐标
%handledata为输入的原始数据（未经过[500k、10m]的fir滤波）
%bias为提取脉冲的阈值
%cutoffband为的低通滤波的上截止频率
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   带通滤波（[500k 10m]）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fs=25*1e6;%采样频率
w1=500*1e3;
w2=10*1e6;
w=[w1 w2]/(fs/2);
b=fir1(150,w,'bandpass');%滤波器设置：150阶，[500k 10m]，带通

data_bandpass_fir=filter(b,1,handledata);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   包络线提取：1、取模，2、低通滤波、3、幅值还原
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(nargin==1)
    bias=0.002;
    cutoffband=1e6;
end

data_abs=abs(handledata);%取模

w=cutoffband/(fs/2);
b=fir1(150,w,'low');
data_profile=filter(b,1,data_abs);%低通滤波

data_profile=data_profile*max(data_bandpass_fir)/max(data_profile);%幅值还原，得到包络线数据

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   脉冲提取
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%实现方法： 1、找出所有大于阈值的点
%           2、找出第一个大于阈值的点
%           3、定位第一个脉冲的起始位置（改变单调性）
%           4、定位第一个脉冲的结束位置（改变单调性）
%           5、找出第一个脉冲后的第一个大于阈值的点
%           6、定位第二脉冲
%           7、类推
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dataoverbias=wthresh(data_profile,'h',bias); 
pointoverbias=find(dataoverbias~=0);
pointnumber=length(pointoverbias);%找出所有大于阈值的点

if(sum(dataoverbias>0)==0)%所有点都小于阈值
    error('No pulse over bias found!');
end

dec=[data_profile 0]-[0 data_profile];
mul=[dec 0].*[0 dec];
turnpoint=find(mul<0)-1;%找出所有的转折点

if(turnpoint(1)==1)
    turnpoint(1)=[];
end
if(turnpoint(length(turnpoint))==length(data_profile))
    turnpoint(length(turnpoint))=[];
end%起始点与最后一点均为无效点

if(length(turnpoint>=1))
    for i=1:length(turnpoint)
        validturn=turnpoint(i);
        if(dec(validturn)>=0)
            turnpoint(i)=0;
        end
    end
end

turnpoint(find(turnpoint==0))=[];%找出所有有效转折点

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       方法一
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
j=1;
for i=1:max((length(turnpoint)-1),1)
    currentpulse=data_bandpass_fir(turnpoint(i):turnpoint(i+1));
    if(max(abs(currentpulse))>=bias)
        pulse(1:length(currentpulse),j)=currentpulse';
        pulse_place(:,j)=[turnpoint(i);turnpoint(i+1)];
        j=j+1;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       特征值计算子程序                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%特征提取：
%波形参素（正负峰值的幅值及位置、正负次峰值的衰减率、过阈值率（检验震荡性）、脉宽、前沿占脉宽比、波形因子：峰值/脉宽）、
%能量参数（均值，能量和，平均能量）、特征谱峰、
%频率参数（傅立叶变换后的频谱峰值和谱峰位置）、
%自相关模型参数、
%AR模型参数。

%Exam:data_chara=chara(data)
%
function charaout=chara(datain)
[r,q]=size(datain);
for i=1:q
    a=datain(:,i);
    if(max(abs(a))~=max(a))
        a=-a;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%  频率因素  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    freq_a=fft(a);
    [freq_max(i),freq_place(i)]=max(abs(freq_a));
    freq_angle(i)=angle(freq_a(freq_place(i)));%频谱峰值和所在位置
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%  自相关因素  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    cov_a(i)=cov(a);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%  波形因素  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    sizea=max(find(a~=0));
    [max_value(i),max_point(i)]=max(a);%正峰值及位置
    [min_value(i),min_point(i)]=min(a);%负峰值及位置
    max_valuevstime(i)=max_value(i)/sizea;%正波形因子
    min_valuevstime(i)=min_value(i)/sizea;%负波形因子
    pulse_width(i)=sizea;%脉冲宽度
    ascend_occupy(i)=max_point(i)/sizea;%上升沿占脉宽比
    a1=[0;a];
    a2=[a;0];
    pass_zero(i)=sum((a1.*a2)<0)/sizea;%过零率
    dec=a2-a1;
    mul=[dec;0].*[0;dec];
    surge(i)=sum(mul<0)/sizea;%振荡率

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%  能量因素  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    sizea=max(find(a~=0));
    aa=a(1:sizea);
    mean_value(i)=mean(abs(aa));%均值
    energy=aa.^2;
    mean_energy(i)=mean(energy);%平均能量
    sum_ener(i)=sum(energy);%能量和
end

charaout=[      max_value;    max_point;      min_value;      min_point;      max_valuevstime;      min_valuevstime; 
%                 正峰值         位置          负峰值           位置              正波形因子          负波形因子
                pulse_width;     ascend_occupy;     pass_zero;      freq_max;     freq_place;      freq_angle;     mean_value;
%                  脉宽              上升比           过零率          频峰           位置              频角           均值
                mean_energy;        sum_ener;       cov_a;           surge];
%                 平均能量          能量和          自相关因子      振荡率    


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                    特征值归一化子程序                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%对提取的特征值进行归一化处理
function charaout=charaguiyi(charain,data_max)
[r,q]=size(charain);
charaout=charain;
for i=1:q
    charaout(:,i)=charaout(:,i)./data_max;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                      神经网络仿真子程序                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%replace of a=sim(net,p)
%net:net structure input
%p:neural network input
%a:neural network output
function a=netsim(net,p)
if(nargin<2)
    error('wrong number of input');
end
[dim num]=size(p);%dim:dimention--输入向量的特征维数;num:输入向量的个数
b1=[];
b2=[];
b3=[];
B1=[];
B2=[];
B3=[];%bias initialize
[numLayer ans]=size(net.IW);%numLayer:the number of the layers
IW1=net.IW{1};%bias and weight matrix between input matrix and the first layer
b1=net.b{1};
IW2=net.LW{2};%bias and weight matrix between the first layer and the second
b2=net.b{2};
IW3=net.LW{numLayer*1+3};%bias and weight matrix between the second layer and the third
b3=net.b{3};
for i=1:num
    B1=[B1 b1];
    B2=[B2 b2];
    B3=[B3 b3];
end
a1=tansig(IW1*p+B1);
a2=tansig(IW2*a1+B2);
a=tansig(IW3*a2+B3);