clear all
close all
clc
x=1:50;
xc=length(x);
x0=25;%%球心位置
h=10;
for i=1:xc
d1(i)=10*h/((x(i)-x0)^2+h^2)^(3/2);%%G*m=10
end
a=d1*500+randn(1,xc)*1;%%randn(a,b)返回正态分布随机数a×b的矩阵%%加了随机噪声原始异常扩大500倍
% load durial.txt
% a=durial(:,4)';
xc=length(a);
% a(20:30)=;
%%%%%kuobian
c2=2^nextpow2(xc);%%取比xc大的最邻近的2^n值
D=c2-xc;
zc=(a(3)- a(1))/3;
yc=(a(xc-3)-a(xc))/3;
ks=(D+1)/2;
for k=1:ks
kz(k)=a(1)-k*zc;
ky(k)=a(xc)-k*yc;
end
kzb=fliplr(kz);
za=[kzb a ky];
xc=length(za);%%扩编
%%%%%%
x1=linspace(0,pi,c2/2);
x2=linspace(-pi,pi,c2);
dfy=fft(za);
pw2=log(abs(dfy));%%abs绝对值%%功率谱
df=fftshift(dfy);
pw=log(abs(df));%%功率谱
k0=0.15*pi;%%截止频率k0
n2=8;%%滤波器阶数为8
for i=1:c2/2
k(i)=1/(1+(x1(i)/k0)^n2);%%巴特沃斯滤波器
end
k1=fliplr(k);
k=[k1 k];
dff=df.*k;%%.*矩阵元素对应相乘
dff=ifftshift(dff);
pw3=log(abs(dff));
vdf=real(ifft(dff));%%取实
vdf1=vdf(ks:xc-ks-1);
subplot 321
plot(a);
title('球体重力异常加噪声主剖面曲线')
subplot 322
plot(za);
title('扩边后的异常曲线')
subplot 323
plot(pw2);
title('扩边后异常曲线的功率谱')
subplot 324
x=x2;
[AX,H1,H2] = plotyy(x,pw,x,k,'plot');
set(AX(1),'XColor','k','YColor','b');
set(AX(2),'XColor','k','YColor','r');
title('移位后的功率谱和巴特沃斯型滤波器')
subplot 325
plot(pw3);
title('滤波反移位后的功率谱')
subplot 326
plot(d1*500);
title('未加噪声前异常曲线和频率域滤波后异常曲线对比')
hold on
plot(vdf1,'r')
