mmfilter.m

基于matlab的探地雷达信号处理开源程序

function  dmf = mmfilter( d, fmode )
%
%   MFILT : Applies a 2-D spatial smoothing filter to the 2-D data in d, depending on 'FMODE'
%
%   Usage : dmf  = mfilt( d, fmode )
%
%   fmode : MEDIFILT : Applies median filter by sliding a user defined 2-D window over a 2-D 
%                      data wall. The data corresponding to the central element of the window 
%                      data is substituted for the median of the data in the window. 
%           MEANFILT : Applies mean filter by sliding a user defined 2-D window over a 2-D 
%                      data wall. The data corresponding to the central element of the window 
%                      data is substituted for the mean of the data in the window. 
%
%      ==>  Zero padding equal to 1/2 the window length is applied, therefore the edges of the 
%           output data are expected to be distorted. 
%      ==>  An external Fortran 90 program does the numerical work for the sake of speed but
%           the (much slower) native MATLAB code is also included (commented out)
%
%  Author    : Andreas Tzanis,
%              Dept. of Geophysics,
%              University of Athens
%              [e-mail]
%  Created   : September 2003
%

%%%%%   Get data size
[ns,ntr]=size(d);

%%% Test filter mode
ftest = strcmp(fmode,'meanfilt') || strcmp(fmode,'MEANFILT') || ...
    strcmp(fmode,'medifilt') || strcmp(fmode,'MEDIFILT');
if ~ftest,
    erh = errordlg('Wrong filter mode - Aborting!','MFILTER : ERROR');
    uiwait(erh);
    dmf = [];
    return
end

%%%%%  Get size of filter window 
clb    = cell(2,1);                                        
clb(1) = cellstr(['Dimension along the VERTICAL   AXIS in Number of samples, >=1']);
clb(2) = cellstr(['Dimension along the HORIZONTAL AXIS in Number of traces,  >=1']);
answer = inputdlg(clb,'Give Filter Dimensions',1 );          
if isempty(answer),                                        % operation was canceled
    dmf = [];
    return
end
fdim = str2num(char(answer));                               
%%%%%   Check filter parameters
if fdim(1) < 1 | fdim(2) < 1,
    erh = errordlg(['Window size cannot be less than 1 element long! ' ...
            'Please try again!'],'MFILTER : ERROR');
    uiwait(erh);
    dmf = [];
    return
end
if fdim(1) > ns/2 | fdim(2) > ntr/2,
    erh = errordlg(['Window size cannot be more than half the data size!' ...
            ' Please try again!'],'MFILTER : ERROR');
    uiwait(erh);
    dmf = [];
    return
end
%%%%%    Define window size around the center
if rem(fdim(1),2) ~=0,
    nslow  = ceil(fdim(1)/2);
    nshigh = floor(fdim(1)/2);
else
    nslow  = fdim(1)/2;
    nshigh = fdim(1)/2;
end
if rem(fdim(2),2) ~=0,
    ntrlow  = ceil(fdim(2)/2);
    ntrhigh = floor(fdim(2)/2);
else
    ntrlow  = fdim(2)/2;
    ntrhigh = fdim(2)/2;
end
clear fdim;
%%%%%   Pad data with zeros
dpad = [ zeros(ns,ntrlow)    d    zeros(ns,ntrhigh)];
dpad = [ zeros(nslow,ntrlow+ntr+ntrhigh); dpad ;  zeros(nshigh,ntrlow+ntr+ntrhigh)];
%%% Compute image using the external Fortran 90 program "mfilt.exe" 
%%% through the driver subfunction "external_f90"
dmf = external_f90(dpad,ns,nslow,nshigh,ntr,ntrlow,ntrhigh,fmode);   %%
%%%%%   Restore original data  size and wrap up
dmf = dmf(nslow+1:ns+nslow,ntrlow+1:ntr+ntrlow);            
return
%
function dmf = external_f90(dpad, ns, nslow, nshigh, ntr, ntrlow, ntrhigh, fmode)
%
% Driver for the external Fortran 90 program MFILT.EXE, which will  
% do the number crunching 
%
%  Author     : Andreas Tzanis,
%               Dept. of Geophysics,   University of Athens
%

%%%%  Export the padded data 
OutDir = tempdir;
fid    = fopen([OutDir 'formfilt.dat'],'w');
fwrite(fid,fmode,'char');
fwrite(fid,ns,'int');      
fwrite(fid,nslow,'int');
fwrite(fid,nshigh,'int');
fwrite(fid,ntr,'int');
fwrite(fid,ntrlow,'int');
fwrite(fid,ntrhigh,'int');
for i=1:nslow+ns+nshigh; 
    fwrite(fid,dpad(i,:),'float'); 
end
fclose(fid);
%%%%%%   Run external program to apply filter 
h = helpdlg('Please wait - External MFILT running','InfoBox');
mfiltpath = which('mfilt.exe');
dos([mfiltpath ' ' OutDir 'formfilt.dat  ' OutDir 'frommfilt.dat'],'-echo')
%%%%  Import filtered data 
fid=fopen([OutDir 'frommfilt.dat'],'r');
dmf=fread(fid,[nslow+ns+nshigh,ntrlow+ntr+ntrhigh],'float');
fclose(fid);
delete([OutDir 'frommfilt.dat']);
close(h);
return