FK滤波效果分析matlab

clearclc;
%%地震数据的读取
fid=fopen（‘F:\Documents\sample‘‘rb‘‘ieee-be‘）;
head1=fread（fid3200‘*int8‘）;%Ascii区
head2=fread（fid100‘*int32‘）;%二进制区

    Dao=double（zeros（1257））;
for i=1:1:57
    daohead（:i）=fread（fid60‘*int32‘）;%道头
    d（1:1500i）=fread（fid1500‘*float‘‘ieee-be‘）;%道数据
end
figure（1）;
wigb（d（::）1）;
set（gca‘ydir‘‘reverse‘）;
set（gca‘xaxislocation‘‘top‘）;
xlabel（‘道数‘）;
ylabel（‘时间（ms）‘）;
set（gca‘ylim‘[01500]）;
set（gca‘xlim‘[054]）;
%F-K域分析
M=2048;N=64;  %采样点数   
y=fft2（d（::）MN）;  %二维傅里叶变换
yy=fftshift（y）;  %正周期变换到主周期
f=0:M;k=0:N;
dt=0.001;  %1ms的时间采样间隔
f=f/M.*（1/dt）; 
dx=5;  %5m的道间距
k=k/N.*（1/dx）; 
figure（2）
subplot（121）;
imagesc（kfabs（y））;  %取绝对值进行成像
set（gca‘YDir‘‘normal‘）;
xlabel（‘k‘）;
ylabel（‘f（Hz）‘）;
subplot（122）;
imagesc（kfabs（yy））;  %取绝对值进行成像
set（gca‘YDir‘‘normal‘）;
xlabel（‘k‘）;
ylabel（‘f（Hz）‘）;




%滤波器的设计及滤波处理
filter=ones（MN）;
for i=1:N
    for j=1:M
        if j>2*（i-N/2）+M/2&j>（-2）*（i-N/2）+M/2
            filter（ji）=0;
        elseif j<2*（i-N/2）+M/2&j<（-2）*（i-N/2）+M/2
            filter（ji）=0;
        end
    end
end
figure（3）
subplot（131）;
imagesc（kffilter）
set（gca‘YDir‘‘normal‘）
xlabel（‘k‘）;
ylabel（‘f（Hz）‘）;
%% f-k滤波  
Y=yy.*filter;  %频域相乘
FY=yy.*（~filter）; 
figure（3）;
subplot（132）;
imagesc（kfabs（Y））;
set（gca‘YDir‘‘normal‘）;
xlabel（‘k‘）;
ylabel（‘f（Hz）‘）;
figure（3）;
subplot（133）;
imagesc（kfabs（FY））;
set（gca‘YDir‘‘normal‘）;
xlabel（‘k‘）;
yla