time   calls  line
                  1 function d = filter_data(d,switchdata,channel,f1,f2,f3,f4,plots,t_start,t_stop, plotOption,verbose)
                  2 
                  3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                  4 %
                  5 %  function d = filter_data(d,switchdata,channel,f1,f2,f3,f4,plots,t_start,t_stop, plotOption)
                  6 %
                  7 %  Filters the specified data using a Tukey band stop filter
                  8 %
                  9 %  -- takes fft of data, applies filter in frequency domain and fft's back
                 10 %  -----          ------ 
                 11 %      |\        /|      
                 12 %      | \      / |      
                 13 %      |  ------  |      
                 14 %      |  |    |  |      
                 15 %     f1 f2   f3 f4       
                 16 %  Notes:
                 17 %  
                 18 %  Takes sample length by looking at difference between adjacent timestamps
                 19 %  N.B. need to check this as I have found slightly different values
                 20 %  depending on the timestamps selected. Might be better to take length of
                 21 %  data/no. samples
                 22 % 
                 23 %  For power spectrum, it plots y on log scale and x as linear
                 24 %  h is the handle to the plot
                 25 %
                 26 % Inputs   d: complete data set
                 27 %          switchdata: =1 for switched data, 0 for combined dataT
                 28 %          channel: channel number
                 29 % (Optional)
                 30 %          f1,f2,f3,f4 : frequency positions of the tukey filter
                 31 %           plots: 1 if you want to see the plots
                 32 %          t_start:  start time (in case you only want to plot a subset of
                 33 %          the data) 
                 34 %          t_stop  : stop_time
                 35 %          plotOption:  string of color, etc.
                 36 %          verbose: =0 if you don't want to see the ascii art
                 37 % 
                 38 % Outputs  filtered data is written out in either d...switchData or
                 39 %            d...dataT
                 40 %            (Need to change this when we decide how to use it in the pipeline)
                 41 %           
                 42 %
                 43 %  act 3/9/2010
                 44 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 45 
            72   46 if(nargin<11) 
            72   47     plotOption = 'b'; 
            72   48     verbose=0; 
            72   49 end 
                 50 
            72   51 if(nargin<9) 
            72   52   plotOption = 'b'; 
            72   53   t_start = 1; 
            72   54   t_stop = length(d.antenna0.receiver.utc); 
            72   55   verbose =0; 
            72   56 end 
                 57 
            72   58 if(nargin<8) 
                 59   plotOption = 'b';
                 60   t_start = 1;
                 61   t_stop = length(d.antenna0.receiver.utc);
                 62   plots = 1
                 63   verbose = 0;
                 64 end
                 65 
            72   66 if (switchdata == 1) 
  0.33      72   67     data = d.antenna0.receiver.switchData(:,channel); 
                 68 else
                 69     data = d.antenna0.receiver.dataT(:,channel);
                 70 end
                 71     
                 72 % First plot the power spectrum of the exisiting data
                 73 
  1.02      72   74 t_sec = (d.antenna0.receiver.utc(t_start:t_stop)-d.antenna0.receiver.utc(t_start))*24*60*60; 
  0.22      72   75 t_sec = t_sec - t_sec(1); 
            72   76 Fs = 1/(t_sec(3)-t_sec(2)); 
  0.61      72   77 L= length(data(t_start:t_stop)); 
  0.02      72   78 NFFT = 2^nextpow2(L); 
 11.36      72   79 y_power_spec = (fft(data(t_start:t_stop),NFFT)); 
  1.55      72   80 x_power_spec = Fs/2*linspace(0,1,NFFT/2+1); 
                 81 
            72   82 if (plots==1) 
                 83 eval(sprintf('h = semilogy(x_power_spec, 2*abs(y_power_spec(1:NFFT/2+1)),''%s'');', plotOption));
                 84 title('Power spectrum of original data set')
                 85 end 
                 86 
                 87 % Now ask what frequencies you want to filter
            72   88 if(verbose==1) 
                 89 disp(' ')
                 90 disp(' A Tukey band stop filter will be applied')
                 91 disp(' ')
                 92 disp(' -----          ------ ')
                 93 disp('     |\        /|      ')
                 94 disp('     | \      / |      ')
                 95 disp('     |  ------  |      ')
                 96 disp('     |  |    |  |      ')
                 97 disp('     f1 f2   f3 f4      ') 
                 98 disp(' ')
                 99 end
            72  100 if (nargin==6) 
                101     disp('Using the following frequencies: ')
                102     f1
                103     f2
                104     f3
                105     f4
                106 end
  0.01      72  107 if (nargin<6) 
                108     f1 = input('Please enter values for f1: ');
                109     f2 = input('Please enter values for f2: ');
                110     f3 = input('Please enter values for f3: ');
                111     f4 = input('Please enter values for f4: ');
                112 end
                113 
                114 % Now call the tukey_bandstop function to create the filter
                115 
  5.55      72  116 filter = Tukey_filter(f1,f2,f3,f4, Fs, NFFT, x_power_spec); 
                117 
                118 % Now multiply the y_power_spectrum by the filter
                119 
  1.30      72  120 filtered_fft = y_power_spec.* filter'; 
                121 
                122 
            72  123 if (plots==1) 
                124     figure
                125     eval(sprintf('h = semilogy(x_power_spec, 2*abs(filtered_fft(1:NFFT/2+1)),''%s'');', plotOption));
                126     title('Power spectrum of the filtered data')
                127 end
                128 
                129 % Now go back to the time domain 
                130 
 22.11      72  131 dummy = real(ifft(filtered_fft,NFFT)); 
  0.66      72  132 dummy = dummy(1:L); 
                133 
                134 
                135 % Now plot the original and filtered data 
                136 
            72  137 if (plots ==1) 
                138 figure
                139 ax(1) = subplot(3,1,1);
                140 plot(data)
                141 title('Original data stream')
                142 ax(2) = subplot(3,1,2);
                143 plot(dummy)
                144 title('Filtered data stream')
                145 ax(3) = subplot(3,1,3);
                146 plot(data - dummy)
                147 title('Residuals')
                148 
                149 linkaxes(ax,'x');
                150 end
                151 
                152 % Finally write back to the data array d
                153 
            72  154 if (switchdata==1) 
                155 %    disp('writing back to to d....switchData')
  0.33      72  156     d.antenna0.receiver.switchData(:,channel) = dummy; 
                157 else
                158 %     disp('writing back to d....dataT')
                159     d.antenna0.receiver.dataT(:,channel) = dummy;
                160 end
                161 
  0.03      72  162 clear dummy; 

Other subfunctions in this file are not included in this listing.