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