both_filter

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function d = both_filter(d,f1, f2, taper, plotflag)

  Function which filters the data wiht both a high and low-pass filter
  It just takes the FFT of the data, applies the filter, and FFTS back.
   
  Conforms with the convention in Angela's tukey filter.

 Inputs   d: complete data set
          f1: freq1 (if low-pass, where it starts to taper, if hi, where it
          starts to rise)
          f2: freq2 (if low, where it is 0, if hi, where it's 1);
          taper: 'lin', 'cos', 'cos3' (cos^3 is bad)
          plotflag:  whether to show plots or not

 Outputs  d with the "data" field modified

 Defaults:  hi-pass, f1=4.5Hz, f2=4.6Hz, taper is gaussian, will plot

  sjcm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function d = both_filter(d,f1, f2, taper, plotflag)
0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0003 %
0004 %  function d = both_filter(d,f1, f2, taper, plotflag)
0005 %
0006 %  Function which filters the data wiht both a high and low-pass filter
0007 %  It just takes the FFT of the data, applies the filter, and FFTS back.
0008 %
0009 %  Conforms with the convention in Angela's tukey filter.
0010 %
0011 % Inputs   d: complete data set
0012 %          f1: freq1 (if low-pass, where it starts to taper, if hi, where it
0013 %          starts to rise)
0014 %          f2: freq2 (if low, where it is 0, if hi, where it's 1);
0015 %          taper: 'lin', 'cos', 'cos3' (cos^3 is bad)
0016 %          plotflag:  whether to show plots or not
0017 %
0018 % Outputs  d with the "data" field modified
0019 %
0020 % Defaults:  hi-pass, f1=4.5Hz, f2=4.6Hz, taper is gaussian, will plot
0021 %
0022 %  sjcm
0023 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0024 
0025 if(nargin<2)
0026   f1 = 4.5;
0027   f2 = 4.6;
0028   taper = 'gauss';
0029   plotflag = 1;
0030 end
0031 
0032 
0033 % get time information from data itself.
0034 t  = (d.antenna0.receiver.utc-d.antenna0.receiver.utc(1))*24*60*60;
0035 data = d.antenna0.receiver.data;
0036 numChans = size(data,2);
0037 dt = median(diff(t));
0038 Fs = 1/dt;
0039 L  = size(data,1);
0040 NFFT = 2^nextpow2(L);
0041 
0042 % next we start to loop over the channels
0043 for m=1:numChans
0044   y_power_spec = (fft(data(:,m),NFFT));
0045   x_power_spec = Fs/2*linspace(0,1,NFFT/2+1);
0046 
0047   
0048   if(plotflag)
0049     figure(1)
0050     setwinsize(gcf, 1600, 1200);
0051     subplot(2, ceil(numChans/2), m)
0052     h = semilogy(x_power_spec, 2*abs(y_power_spec(1:NFFT/2+1)), 'b');
0053     hold on
0054   end
0055   
0056   % set up our filter.
0057   % find freq boundaries in index
0058   half_filter = ones(1,NFFT/2);
0059   ff1 = find(abs(x_power_spec - f1) == min(abs(x_power_spec - f1)));
0060   ff2 = find(abs(x_power_spec - f2) == min(abs(x_power_spec - f2)));
0061   transL = ff2 - ff1;
0062   width  = x_power_spec(ff2) - x_power_spec(ff1);
0063   
0064   % defining the hi-pass ones.
0065   switch taper
0066     case 'lin'
0067       filtline = (1)/(x_power_spec(ff2)-x_power_spec(ff1))*(x_power_spec - x_power_spec(ff1)) + 0;
0068     
0069     case 'cos'
0070       filtline = 0.5*(1+cos(pi*x_power_spec/width));
0071       
0072     case 'cos3'
0073       % really a cos^3
0074       filtline =  0.5*(1+(cos(pi*x_power_spec/width)).^3);
0075   end  
0076   half_filter = filtline;
0077   half_filter(1:ff1) = 0;
0078   half_filter(ff2:length(filtline)) = 1;
0079   half_filter(length(filtline)) = [];
0080   
0081   sec_half = fliplr(half_filter);
0082   full_filter = [half_filter sec_half];
0083   
0084   hi_filter = full_filter;
0085   lo_filter = 1-full_filter;
0086   
0087   % now we apply the filter
0088   hi_filtered_fft = y_power_spec.* hi_filter';
0089   lo_filtered_fft = y_power_spec.* lo_filter';
0090   if(plotflag)
0091     subplot(2, ceil(numChans/2), m)
0092     h = semilogy(x_power_spec, 2*abs(hi_filtered_fft(1:NFFT/2+1)), 'r');
0093     title('Power spectrum')
0094     legend('Unfiltered', 'Filtered')
0095   end
0096 
0097   % return to time domain
0098   hiTimeDom = real(ifft(hi_filtered_fft,NFFT));
0099   hiTimeDom = hiTimeDom(1:L);
0100   loTimeDom = real(ifft(lo_filtered_fft,NFFT));
0101   loTimeDom = loTimeDom(1:L);  
0102   
0103 
0104   % plot the two and the residuals
0105   if(plotflag)
0106     figure(2)
0107     setwinsize(gcf, 1600, 1200);
0108     subplot(2, ceil(numChans/2), m)
0109     plot(t, data(:,m));
0110     hold on
0111     plot(t, hiTimeDom, 'r');
0112     txt = sprintf('Channel %d', m);
0113     title(txt);
0114     legend('Original', 'Filtered');
0115 
0116     figure(3)
0117     setwinsize(gcf, 1600, 1200);
0118     subplot(2, ceil(numChans/2), m)
0119     plot(t, data(:,m) - hiTimeDom);
0120     title('Residual');
0121   end
0122 
0123   d.antenna0.receiver.lo(:,m) = loTimeDom;
0124   d.antenna0.receiver.hi(:,m) = hiTimeDom;
0125 end
0126 
0127 return;
0128