flat_events.m

matlab下KL变换,用于去除随机噪声

clc; clear; close all  
  dt = 2./1000;
  tmax = 0.8;
  h = [0:10:10*(80-1)];
  tau = [0.1,0.2,0.6]; 

  p = [0,0.0004,-0.0008];
  amp = [0.3,0.3,0.3];
  f0 = 40;
  snr = 1.5;     %信噪比
  L = 5;
 
 nt = floor(tmax/dt)+1;
 nfft = 4*(2^nextpow2(nt));
 n_events = length(tau);
 nh = length(h);
 wavelet = ricker(f0,dt); 
 nw = length(wavelet);
 W = fft(wavelet,nfft);
 D = zeros(nfft,nh);
 i = sqrt(-1);

 delay = dt*(floor(nw/2)+1);   %半个子波延续时间
 
 for ifreq=1:nfft/2+1
  w = 2.*pi*(ifreq-1)/nfft/dt;
    for k=1:n_events
        Shift = exp(-i*w*(tau(k)+h*p(k)-delay));
        D(ifreq,:) = D(ifreq,:) + amp(k)* W(ifreq)*Shift;
    end
 end

% w-domain symmetries

 for ifreq=2:nfft/2
  D(nfft+2-ifreq,:) = conj(D(ifreq,:));
 end 

 d = ifft(D,[],1);
 d = real(d(1:nt,:));
 t = (0:1:nt-1)*dt; 
 dmax  = max(max(abs(d)));
 if 1
    op = hamming(L);
    Noise = conv2(randn(size(d)),op,'same');
    Noise=randn(size(d));
    Noisemax = max(max(abs(Noise)));
 
    d = d + Noise*(dmax/Noisemax)/snr; %由此式计算snr,就是信噪比
 end

 wigb(d,1,h,t)
 
 
  dip1=2;

[nt,nh]=size(d)
out1=zeros(nt,nh);
 
    
for k=1:nh

    for j=1:nt
         qq=j+dip1*(k-1);
         if qq>0
           if qq<nt+1 
             out1(j,k)=d(qq,k);
           end
         end 
     end
 end
%end
figure
wigb(out1,1,h,t)

 dip2=-4;

[nt,nh]=size(d)
out2=zeros(nt,nh);
 


%for i=1:Len
    
for k=1:nh

    for j=1:nt
         qq=j+dip2*(k-1);
         if qq>0
           if qq<nt+1 
             out2(j,k)=d(qq,k);
           end
         end 
     end
 end
%end
figure
wigb(out2,1,h,t)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%kl滤波
d2=kl(d,1);
out12=kl(out1,1);
out22=kl(out2,1);

dip3=-2;

[nt,nh]=size(out12)
out12_o=zeros(nt,nh);
 
    
for k=1:nh

    for j=1:nt
         qq=j+dip3*(k-1);
         if qq>0
           if qq<nt+1 
             out12_o(j,k)=out12(qq,k);
           end
         end 
     end
 end
%end
figure
wigb(out12_o,1,h,t)

 dip4=4;

[nt,nh]=size(out22)
out22_o=zeros(nt,nh);
 


%for i=1:Len
    
for k=1:nh

    for j=1:nt
         qq=j+dip4*(k-1);
         if qq>0
           if qq<nt+1 
             out22_o(j,k)=out22(qq,k);
           end
         end 
     end
 end
%end
figure
wigb(out22_o,1,h,t)

figure
wigb(d2+out12_o+out22_o,1,h,t)

figure
wigb([d,d2+out12_o+out22_o],1,0:10:(2*nh-1)*10,t)