fkdesign.m

这是matlab在地球物理数据处理方面的源码

function [newtraces,mask] = fkdesign(traces,t,offset)

% function [newtraces,mask] = fkdesign(traces,t,offset)
% Program to plot a set of traces in the f-k domain then design a filter
% which is manually picked using the mouse.
% traces are the data in column oriented format
% t is the timebase in seconds, this is a vector of the same length as each trace.
% offset is the offset in metres and is a vector equal to the number of traces
% 
% You will probably need to run this a number of times iteratively updating
% what you have surgically removed.
%
% The f-k image shown may first be zoomed if required, then picking begins
% when any key is hit.
%
% User may begin selection of area immediately then finish polgonal region
% picking by a shift-click,  a right-click, or a double-click adds
% final vertex to the selection and then starts making the mask which
% is then applied to the data.  The program then outputs the filtered
% data plus the mask used.
%
% A hamming taper is used to mask the data.  The data are further padded with 
% zeros to the next highest power of 2 - this implementation runs more 
% quickly than earlier ones.
%
% Studenst should feel free to modify the program to their own desires.
% One improvement centres on the tapering of the filter mask.

% Written by D. Schmitt, last modified January 27, 2000.

[m,n] = size(traces); 
close all
fclose('all');

[x,y] = meshgrid(hamming(n),hamming(m));  tracemask = x.*y; clear x; clear y;
traces = tracemask.*traces;  % smooth out some rough edges for ffting later, applied twice.

M = pow2(nextpow2(m)); N = pow2(nextpow2(n));

delt = t(2)-t(1); fnyq = 1/(2*delt); delf = 2*fnyq/M;
freqs = -fnyq:delf:fnyq - delf;
delx = abs(offset(2)-offset(1)); knyq = 2*pi/(2*delx); delk = knyq/N;
ks = -knyq:delk:knyq-delk;

newtraces = zeros(M,N);  

% moveout = 1000000000 %1850;  % Moveout velocity in m/s if want to flatten
% intercept =  0 % 0.008 % Time intercept for linear moveout% ind = round((intercept + offset/moveout)/delt)+1;
% for i = 1:n
% newtraces(m/2+(1:m-ind(i)+1),n/2+i) = traces(ind(i):m,i);
% end
% newtraces(m/2+1:3*m/2,n/2+1:3*n/2) = tracemask.*newtraces(m/2+1:3*m/2,n/2+1:3*n/2);
% imagesc(newtraces); This section may be reincluded if one wants to preferentially 
% flatten on a given events as may be done in the vsp data.
% clear traces ind
indxy = [(M-m)/2+1  m+(M-m)/2  (N-n)/2+1  n+(N-n)/2 ]; % Index corners of data within newtraces
%newtraces((M-m)/2+1:m+(M-m)/2,(N-n)/2+1:n+(N-n)/2) = traces;
newtraces(indxy(1):indxy(2),indxy(3):indxy(4)) = traces;
clear traces


z = fft2(newtraces); z = fftshift(z); 
clear newtraces

imagesc(ks,freqs(1:M/2+1),abs(z(1:M/2+1,:)));

ylabel('frequency (Hz)'); xlabel('Spatial Wavenumber (radians/m)')
title('Zoom if necessary then hit any key to begin picking, end picking by double click')
pause

[mask,xi,yi] = roipoly;
mask = not(mask);  % change area with 1's to 0's and vice-versa
%mask(1:m+1,n+1:2*n) = 0*mask(1:m+1,n+1:2*n);  % This includes zero frequency
%mask(m+2:2*m,1:n+1) = 0*mask(m+2:2*m,1:n+1); % Also includes zero here.
mask(M/2+1:M,1:N) = fliplr(flipud(mask(1:M/2,1:N)));

kernel = fspecial('average',[30 7]);  %Smooth the edges of the filter
mask = conv2(mask,kernel,'same');

z = mask.*z;  
z = fftshift(z); newtraces = real(ifft2(z)); 
clear z
newtraces = newtraces(indxy(1):indxy(2),indxy(3):indxy(4))./tracemask;