ex.rar_ex

%带通子波 clc clear % zf=[]; % zt=[]; % zt0=[]; % a=5; % b=5*2^(1/2); % d=100; % c=d/2^(1/2); % f1=[0:9999]*0.1; % zf=(f1<a|f1>d)*0+(f1>=a&f1<=b).*sin(pi*(f1-a)/(2*(b-a)))+(f1>b&f1<c)*1+(f1>=c&f1<=d).*cos(pi*(f1-c)/(2*(d-c))); % t1=[-5000:4999]*0.001; % zt=ifft(zf); % zt0=ifftshift(real(zt)); % zt0=zt0/max(zt0); % subplot(122),plot(f1,zf),title('5-100Hz带通频谱'),xlabel('f/Hz'),ylabel('\itAMP'),legend('带通频谱'),axis([0 150 0 1.2]) % subplot(121),plot(t1,zt0),title('5-100Hz带通子波'),xlabel('t/s'),legend('带通子波'),axis([-0.1 0.1 -0.4 1.2]) % % % %求主旁瓣能量比 % % P=sum(abs(zt0(4994:5000)))/sum(abs(zt0(4988:4994))) % % Q=sum(abs(zt0(4994:5000)))/sum(abs(zt0(4988:5000))) % % % % % %B_L小波作积分 % yt=[]; % yf=[]; % yt0=[]; % yf0=[]; % yt1=[]; % yt2=[]; % yf1=[]; % % % m=10; % n=2; % yf(1)=0; % for f=[1:9999]*0.1; % yf1=@(g)16*(sin(pi*0.4*f./g)).^4.*((1+2*(sin(pi*0.4*f./g)).^2)./(1-(2/3)*(sin(pi*0.4*f./g)).^2)./(3-8*(sin(pi*0.4*f./g)).^2+8*(sin(pi*0.4*f./g)).^4)).^(1/2)./(2*pi*0.8*f./g).^2./(g+m); % yf(n)=quad(yf1,0.1,80)/(log(80+m)-log(0.1+m)); % n=n+1; % end % f=[0:9999]*0.1; % yt=ifft(yf); % t=[-5000:4999]*0.001; % yf1=abs(yf); % yf0=yf1/max(yf1); % yt1=ifftshift(real(yt)); % yt0=yt1/max(yt1); % % % % % % % % % %对积分后的B_L子波作带通滤波 % xf=[]; % xf0=[]; % xt=[]; % xt0=[]; % xf=yf.*zf; % xf0=abs(xf); % xf0=xf0/max(xf0); % xt=ifft(xf); % xt0=ifftshift(real(xt)); % xt0=xt0/max(xt0); % % % % subplot(231),plot(f,yf0,'g',f,abs(zf),'b',f,xf0,'r'),title('0-80'),xlabel('f/Hz'),ylabel('AMP'),axis([0 150 0 1.2]),hold on % subplot(234),plot(t,xt0,'r'),title('0-80'),xlabel({'t/s';'整形后的带通'}),axis([-0.1 0.1 -0.6 1.2]),hold on % % % % % subplot(122),plot(f,yf0,'g',f,abs(zf),'b',f,xf0,'r'),title('经宽带B_L子波整形后的带通频谱'), % % xlabel('f/Hz'),ylabel('AMP'),text(90,0.8,'5-100Hz带通'),text(120,0.15,'宽带B_L'),text(50,0.2,'整形后的带通'),axis([0 150 0 1.2]),hold on % % subplot(121),plot(t,xt0,'r'),title('经宽带B_L子波整形后的带通剖面'),legend('整形后的带通'),xlabel('t/s'),axis([-0.1 0.1 -0.6 1.2]),hold on % % B_L小波 f=[0:9999]*0.1; t=[-5000:4999]*0.001; g0=1; hf=16*(sin(pi*0.4*f./g0)).^4.*((1+2*(sin(pi*0.4*f./g0)).^2)./(1-(2/3)*(sin(pi*0.4*f./g0)).^2)./(3-8*(sin(pi*0.4*f./g0)).^2+8*(sin(pi*0.4*f./g0)).^4)).^(1/2)./(2*pi*0.8*(f+eps)./g0).^2; ht=ifft(hf); ht0=ifftshift(ht)*(-1)*10^4/6; %ht0=ht0/max(ht0); %hf0=hf/max(hf); subplot(121),plot(t,ht0),xlabel('t/s'),title('Battle-Lemarie子波'),hold on subplot(122),plot(f,hf),axis([0 20 0 1]),xlabel('f/Hz'),ylabel('AMP'),title('频谱图'),hold on % % % % % % 宽带Ricker子波 % % p0=10; % % q0=80; % % hf=(exp(-(f/q0).^2)-exp(-(f/p0).^2))/pi^(1/2); % % ht=(q0*exp(-(pi*q0*t).^2)-p0*exp(-(pi*p0*t).^2))/(q0-p0); % % hf0=hf/max(hf); % % ht0=ht/max(ht); % % % % % % sf=hf.*zf; % % sf0=abs(sf); % % sf0=sf0/max(sf0); % % st=ifft(sf); % % st0=ifftshift(real(st)); % % st0=st0/max(st0); % % % % % %Ricker子波 % % g0=50; % % ht=(1-(pi*g0*t1).^2).*exp(-1*(pi*g0*t1).^2); % % hf=fft(ht)/10; % % % % ht0=ht/max(ht); % % % % % % % % % subplot(221),plot(f,yf0,'k',f,abs(zf),'g',f,hf0,'b'),title({'0.8f/g替代f并除以（g+10）以30-80Hz积分';'红线为滤波前的子波频谱'}), % % xlabel('f/Hz'),legend('宽带B_L小波','5-100Hz带通'),axis([0 150 0 1]),hold on % % subplot(222),plot(f,xf0,'r',f,sf0,'b'),title('经带通滤波后的子波频谱'),legend('宽带B_L带通频谱'),xlabel('f/Hz'),axis([0 150 0 1]),hold on % % subplot(224),plot(t,xt0,'r',t,st0,'b'),title('经带通滤波后的子波时域'),legend('宽带带通B_L'),xlabel('t/s'),axis([-0.1 0.1 -0.6 1]),hold on % % subplot(223),plot(t,yt0,'r',t,ht0,'b'),title('经带通滤波前的子波时域'),legend('宽带B_L'),xlabel('t/s'),axis([-0.1 0.1 -0.6 1]),hold on % % % % figure,subplot(121),plot(f,yf0,'r'),title('0.8f/g替代f并除以（g+10）以10-80Hz积分'), % % xlabel('f/Hz'),legend('宽带B_L频谱'),axis([0 150 0 1]),hold on % % subplot(122),plot(t,yt0,'r'),title('宽带B_L子波时域'),legend('宽带B_L剖面'), % % xlabel('t/s'),axis([-0.1 0.1 -0.6 1]),hold on % % % % % % % % figure,subplot(121),plot(f,zf,'r'),title('5-100Hz带通'), % % xlabel('f/Hz'),legend('带通频谱'),axis([0 150 0 1]),hold on % % subplot(122),plot(t,zt0,'r'),title('带通剖面'),legend('带通剖面'), % % xlabel('t/s'),axis([-0.1 0.1 -0.6 1]),hold on % % % % % % figure,subplot(221),plot(f,yf0,'g',f,abs(zf),'b',f,xf0,'r'),title({'宽带BL对带通整形';'红线为整形后的带通频谱'}), % % xlabel('f/Hz'),text(100,0.8,'5-100Hz带通'),text(110,0.3,'宽带B_L'),text(50,0.2,'整形后的带通'),axis([0 150 0 1]),hold on % % subplot(222),plot(t,zt0,'b'),title('5-100Hz带通剖面'),legend('带通剖面'),xlabel('f/Hz'),axis([-0.1 0.1 -0.6 1]),hold on % % subplot(224),plot(t,xt0,'r'),title('整形后的带通'),legend('整形后的带通'),xlabel('t/s'),axis([-0.1 0.1 -0.6 1]),hold on % % subplot(223),plot(t,yt0,'g'),title('宽带B_L'),legend('宽带B_L'),xlabel('t/s'),axis([-0.1 0.1 -0.6 1]),hold on % % % % % figure,subplot(121),plot(f,yf0,'g',f,abs(zf),'b',f,xf0,'r'),title({'0.1-80Hz宽带BL对带通整形';'红线为整形后的带通频谱'}), % % xlabel('f/Hz'),text(100,0.8,'5-100Hz带通'),text(110,0.3,'宽带B_L'),text(50,0.2,'整形后的带通'),axis([0 150 0 1]),hold on % % subplot(122),plot(t,zt0,'b',t,xt0,'r'),title('整形后的带通'),legend('整形后的带通'),xlabel('t/s'),axis([-0.1 0.1 -0.6 1]),hold on % % % % % % % % % % subplot(233),plot(f,yf0,'g',f,abs(zf),'b',f,xf0,'r'),title('1/(g+10)5-80Hz'), % % xlabel('f/Hz'),axis([0 150 0 1]),hold on % % subplot(236),plot(t,xt0,'r'),title('整形后的带通'),legend('5-80Hz'),xlabel('t/s'),axis([-0.1 0.1 -0.6 1]),hold on % % % %求剖面各种参数 % yt2=yt0(1:5000); % ytmin=min(yt0) % N1=find(yt2==min(yt2)) % dt2=0.002*(5000-N1) % N2=find(yt2<=0.00001); % dt1=0.002*(5000-max(N2)) % maxN2=max(N2) % yt4=yt2(1:N1); % N3=find(yt4<=0.001); % dt3=0.002*(5000-length(N3)) % % % %求频谱各种参数 % yf0=yf0(1:10000); % a=[]; % b=[]; % c=[]; % a=yf0.*yf0; % a0=sum(a(1:end)); % b=a.*f(1:10000); % b0=sum(b(1:end)); % f0=b0/a0 % f1=f(1:10000)-f0*ones(1,10000); % c=f1.^2.*a; % c0=sum(c(1:end)); % df=(c0/a0)^(1/2) % % % ztmin=min(zt0) % % % %求主旁瓣能量比 % % p=sum(abs(yt0(5000-max(N2):5000)))/sum(abs(yt0(4988:5000-max(N2)))) % % q=sum(abs(yt0(5000-max(N2):5000)))/sum(abs(yt0(4988:5000))) % % % % % % % % % subplot(121),plot(f,abs(xf),'r',f,abs(hf),'b'),title('原Ricker子波50Hz和频域10-80Hz积分并经5-80Hz带通滤波后频谱'),legend('积分并经5-80Hz带通滤波后频谱','原Ricker子波50Hz'), % % xlabel('f/Hz'),axis([0 150 0 1]) % % subplot(122),plot(t,xt0,'r',t,ht0,'b'),title('原Ricker子波50Hz和频域10-80Hz积分并经5-80Hz带通滤波后的子波'),legend('积分并经滤波后的子波','原Ricker子波50Hz'), % % xlabel('t/s'),axis([-0.1 0.1 -0.6 1])