f011x2_preprocess_data.m

对所采集到的非线性信号进行数据预处理,消除趋势项

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%                                清理，初始化                            %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% Syntax for clearing
clear            % Remove items from workspace, freeing up system memory
clc              % Clear Command Window
clf              % Clear current figure window
close all hidden % removal draws only those lines that are not obscured
                 % by other objects in the field of view. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%                 读入数据：方法1：用语句 load                            %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% load an matrix (matrix name and file name are the same )
% load  f_111.txt;             % a data file for test
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%                 读入数据：方法2：用语句 fopen                           %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Read formatted data from file
% f_i1=fopen('abs_bisp_h.txt')
% abs_bisp_hh=fscanf(f_i1,'%f',[128 128]);
% fclose(f_i1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%                 读入数据：方法3：用自定义函数 FLRead                     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Read a row data to a vectory
% f_name='f_011.txt';        % f_011.txt 删去第一行中文解释的文件
 f_name='f_111.txt';
% f_name='f_211.txt';
% f_name='f_311.txt';
% f_name='f_411.txt';
% f_name='f_511.txt';

begin=5000;      % 开始行
numline=5002;    % 总行数
Nocolume=2;      % 所选列
x=FLRead(f_name,begin,numline,Nocolume);
                 % 输入参数
                 % s代表文件名，为一字符串 例如: s='1.txt'
                 % begin代表开始行 例如: begin=10
                 % numline代表返回的行数 例如: numline=1024
                 % Nocolume代表返回的列，即返回哪一列 例如: Nocolume=3,
                 % 表示把第三列的数返回给Result
%%%%%%%%%%%% Write formatted data to an ASCII file
x=x/1000;                        % 将原数据缩小
% f_o1=fopen('f011x2.txt','w');     % f_011x1.txt 存原文件原的第1列数据(位移)
 f_o1=fopen('f111x2.txt','w');
% f_o1=fopen('f211x2.txt','w');
% f_o1=fopen('f311x2.txt','w');
% f_o1=fopen('f411x2.txt','w');
% f_o1=fopen('f511x2.txt','w');
for k=1:length(x) 
  fprintf(f_o1,'   %8.3f\n',x(k));
end
fclose(f_o1);                       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%                               消除趋势项                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
delta_width=150;      % 分三段，每段宽度
T=1;
%start_x=120;           % f011x2: 140-620----lilear
start_x=140;         % f111x2: 140-620----lilear
%start_x=150;         % f211x2: 140-620----lilear
%start_x=130;         % f311x2: 130-610----lilear
%start_x=120;         % f411x2: 120-600----lilear
%start_x=130;         % f511x2: 130-610----lilear
area_1=0;
area_3=0;
for i=1:delta_width
    area_1=area_1+x(start_x+i).*T;
    area_3=area_3+x(start_x+2*delta_width+i).*T;
end;    
data_ratio=0.5*(area_3-area_1)/delta_width/delta_width;
for i=1:3*delta_width
    x_trend(i)=data_ratio*i*T;                                      % 趋势项
    x_del_trend(i)=x(start_x+i)-x_trend(i);                  % 消除趋势项信号
end;            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                    中数法: 消除确定性成分                                %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%% 中数法：滤去低频的确定性信号，并获得零均值有色噪声
   ok_data=true;
   k=1;
   while ok_data
      k=k+1;
      x_median(k)=(x(k-1)+2*x(k)+x(k+1))/4;                    % x_median 中数
      x_del_median(k)=x(k)-x_median(k);                      % 消除中数后的数据
      if k==length(x)-1
         ok_data=false
      end;
   end;                   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%                               小波滤波                                       %
% x:原始数据；x0:提取的较高频率的信号；y:低频信号==》xc:滤去低频的时间序列:xc=x-y  % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    scales=3;        % a=2--4
    x_wavelet_high=cwt(x,scales,'db5');                      % 小波提取的高频成分
    
    scales=96;       % a=128
    x_wavelet_low=cwt(x,scales,'db5');                       % 小波提取的低频成分  
    
    for k=1:length(x_wavelet_low);
      x_wavelet_del_low(k)=x(k)-x_wavelet_low(k); % 原数据滤去小波的低频成分后信号 
    end;  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%                        五点滑动平均法平滑处理                                  % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%555%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   moving=23;    %平滑次数不能太大,取[1：3],moving 越大,滤去越多低频部分
   n=length(x);
   a=x;
   for  k=1:moving
      b(1)=(3*a(1)+2*a(2)+a(3)-a(4))/5;
      b(2)=(4*a(1)+3*a(2)+2*a(3)+a(4))/10; 
      for j=3:n-2
         b(j)=(a(j-2)+a(j-1)+a(j)+a(j+1)+a(j+2))/5;
      end
      b(n-1)=(a(n-3)+2*a(n-2)+3*a(n-1)+4*a(n))/10;
      b(n)=(-a(n-3)+a(n-2)+2*a(n-1)+3*a(n))/5;
      a=b;
   end
   x_moving=a;                                         % 五点滑动平均法获得平滑数据	
   x_del_moving=x-x_moving';                      % 滤去五点滑动平均法获得的平滑数据
% x_del_moving=x-x_moving;  % 有时采用（注意 行 或 列 向量）

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%                                绘制图                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%------------------------------------------------------------------%
%                          消除趋势图                               % 
%------------------------------------------------------------------%
    figure(1)  
    %-----------------
    subplot(3,1,1)
    k=1:3*delta_width;
    plot(k,x(start_x+k)); 
    title('原始数据');
    %-----------------
    subplot(3,1,2)
    k=1:3*delta_width;
    plot(k,x(start_x+k),'b',k,x_trend(k),'r--'); 
    title('原始数据与趋势成分比较图');
    %-----------------
    subplot(3,1,3)
    k=1:3*delta_width;
    plot(k,x_del_trend(k));grid
    ylabel(' a=128')
    title('原始数据滤去趋势成分成分后信号');
%------------------------------------------------------------------%
%                          中数法图                                 % 
%------------------------------------------------------------------%
   figure(2)    
   % x=x(2:length(y)+1);
   %----------------绘制原始数据 x-t图----------------------
   subplot(3,1,1);
   k=1:1000;
   plot(k,x(k));        %原始数据：绘制原始数据 x-t图
   title('原始数据');
   %----------------原始数据与确定性成分比较图 --------------
   subplot(3,1,2);
   k=50:150;
   plot(k,x(k),'b',k,x_median(k),'r--'); 
   title('原始数据与中数法获得低频的确定性成分比较图');
   %---------------- 绘制有色噪声---------------------------
   subplot(3,1,3);
   k=1:500;
   plot(k,x_del_median(k));
   title('用中数提取的较高频的信号');% 可以用于系统分析的有色噪声数据 
%------------------------------------------------------------------%
%                    小波法（WAVELET）图                            % 
%------------------------------------------------------------------%
   figure(3)   
   %----------------绘制原始数据 x-t图----------------------
   subplot(5,1,1);
   k=1:1000;
   plot(k,x(k));     %原始数据：绘制原始数据 x-t图
   title('原始数据');
   %-----------------
   subplot(5,1,2)
   k=1:1000;
   plot(k,x_wavelet_high(k));grid
   ylabel(' a=2');
   title('原始数据与小波法获得高频的成分比较图');
   %-----------------
   subplot(5,1,3)
   k=1:1000;
   plot(k,x_wavelet_low(k));grid
   ylabel(' a=96')
   title('原始数据与小波法获得低频的成分比较图');
   %-----------------
   subplot(5,1,4)
   k=1:500;
   plot(k,x_wavelet_del_low(k));grid
   title('原始数据滤去小波法获得的低频的成分后数据');
   %-----------------
   subplot(5,1,5)
   c=cwt(x,[10,30,60,90,120,150],'db5','plot');
   title('小波法图像');
%-------------------------------------------------------------------%
%                    五点滑动平均法平滑图                            % 
%-------------------------------------------------------------------%
    figure(4)  
    %-----------------
    subplot(3,1,1)
    k=1:1000;
    plot(k,x(k)); 
    title('原始数据');
    %-----------------
    subplot(3,1,2)
    k=1:1000;
    plot(k,x(k),'b',k,x_moving(k),'r--'); 
    title('原始数据与五点滑动平均法获得低频成分比较图');
    %-----------------
    subplot(3,1,3)
    k=1:500;
    plot(k,x_del_moving(k));grid
    ylabel('x_del_moving')
    title('原始数据滤去五点滑动平均法的低频成分后信号');  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               计算方差                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   c = xcorr(x);
   the_sigma=sqrt(abs(c));
   % fprintf(1,'the_sigma  % 8.5f\n',the_sigma);
   
   for i=1:length(x)
      if abs(x(i))>=3.*the_sigma;
          x(i)=0;
      end; % if
   end;  % for   
   xt=x';
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                    处理后数据存盘                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 趋势项
save('f111x2_del_trend.txt','x_del_trend','-ASCII');               % 滤去趋势项后的信号
%% 中数
save('f111x2_del_median.txt','x_del_median','-ASCII');                 % 滤去中数后的信号
%% 小波
save('f111x2_wavelet_low.txt','x_wavelet_low','-ASCII');             % 小波提取的低频信号
save('f111x2_wavelet_high.txt','x_wavelet_high','-ASCII');      % 小波提取的较高频率的信号
save('f111x2_wavelet_del_low.txt','x_wavelet_del_low','-ASCII');  % 小波法滤去低频后的信号
%% 五点滑动平均法
save('f111x2_del_moving.txt','x_del_moving','-ASCII');    % 五点滑动平均法滤去低频后的信号