hilo_filter

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

  function d = hilo_filter(d,f1, f2, type, taper, plotflag)

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