processradar.m

％ BackgroundRemoval=[true],false ％ Gain=[tsquare],linear ％ BandPass=[paul],fircls ％ CenterFrequen

function A=processradar(twt,A,varargin)
%
% A=processradar(twt,A [,parameters])
%
%Parameters:
%    BackgroundRemoval=[true],false
%    Gain=[tsquare],linear
%    BandPass=[paul],fircls
%    CenterFrequency, auto   (determined using pburg)
%    BandWidth=auto          (a fraction of the CenterFrequency default=0.25)
%    ContrastStretch=[true],false
%    HilbertAmplitude=[true],false
%    HorizontalStacking=1    (a number of traces)
%
% 
% Aslak Grinsted 2006


Args=struct('BackgroundRemoval',1, ...
    'Gain','tsquare', ...
    'BandPass','paul', ...
    'CenterFrequency', 'auto', ...
    'BandWidth', 'auto', ...
    'ContrastStretch',1, ...
    'HilbertAmplitude',1, ...
    'HorizontalStacking',1);


Args=parseArgs(varargin,Args,{'BackgroundRemoval','Amplitude'});

scans=size(A,2);
N=size(A,1);
dt=mean(diff(twt));

switch lower(Args.Gain)
    case 'linear'
        gain=twt/max(twt);
    case 'tsquare'
        gain=(twt/max(twt)).^2;
    otherwise
        error('Unknown gain method.')
end

bg=0;
if Args.BackgroundRemoval
    bg=mean(A,2);
end

for ii=1:scans
    A(:,ii)=(A(:,ii)-bg).*gain;
end

if strcmpi(Args.CenterFrequency,'auto')
    [P,F]=pburg(detrend(bg.*gain),2,512,1/dt);
    [mx,mxi]=max(P);
    Args.CenterFrequency=F(mxi);
end

if strcmpi(Args.BandWidth,'auto')
    Args.BandWidth=0.25;
end
Args.BandWidth=1-Args.BandWidth;


Npow2=2.^nextpow2(N*1.1); %ensure some zero padding!

bandpass=[];
k = [1:fix(Npow2/2)];
k = k.*((2.*pi)/(Npow2*dt));
k = [0., k, -k(fix((Npow2-1)/2):-1:1)];

complexfilter=0;
switch lower(Args.BandPass)
    case 'paul'
        m=4;
        fourier_factor = 4*pi/(2*m+1);
        scale=1./(Args.CenterFrequency*fourier_factor);
        expnt = -(scale.*k).*(k > 0.);
        norm = sqrt(scale*k(2))*(2^m/sqrt(m*prod(2:(2*m-1))))*sqrt(Npow2);
        bandpass = norm*((scale.*k).^m).*exp(expnt);
        bandpass = bandpass.*(k > 0.);     % Heaviside step function
        %see wave_bases by torrence and compo
        complexfilter=1;
    case 'fircls'
        nfilt=floor(2.5/(Args.CenterFrequency*dt));
        fq=Args.CenterFrequency*dt;
        f=[0 fq*Args.BandWidth fq/Args.BandWidth 1];
        a=[0 1 0];
        up=[0.01 1.02 0.01];
        lo=[-0.01 0.98 -0.01];
        bandpass=fircls(nfilt,f,a,up,lo);
        bandpass(Npow2)=0;
        bandpass=fft(circshift(bandpass(:),-floor(nfilt*.5)));
    case 'none'
    otherwise
        error('Unknown BandPass method.')
end

if ~isempty(bandpass)
    bandpass=bandpass(:);
    for ii=1:scans
        F=bandpass.*fft(A(:,ii),Npow2);
        if Args.HilbertAmplitude&(~complexfilter)
            F(2:Npow2/2)=F(2:Npow2/2)*2; %manual hilbert transform
            F(Npow2/2+2:end)=0;
        end
        F=ifft(F);
        F(N+1:end)=[];
        if Args.HilbertAmplitude
            A(:,ii)=abs(F);
        else
            if ~isreal(F)
                F=real(F);
            end
            A(:,ii)=F;
        end
    end
else
    if Args.HilbertAmplitude
        error('You Should BandPass prior to hilbert amplitude')
    end
end

if Args.HorizontalStacking>1
    scans=floor(scans/Args.HorizontalStacking)*Args.HorizontalStacking;
    A(:,scans+1:end)=[];
    A=permute(mean(reshape(A,size(A,1),Args.HorizontalStacking,[]),2),[1 3 2]);
end

if Args.ContrastStretch
    Q=fix(N/10);ix=[1 (fix(Q/2):Q:N) N];
    mA=mean(abs(A),2);
    f=fit(twt(ix),mA(ix),'smoothingspline');
    regain=1./feval(f,twt);
    for ii=1:size(A,2)
        A(:,ii)=A(:,ii).*regain;
    end
end