filter_data_choice

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

  function d = filter_data(d,switchdata,channel,filter,f1,f2,f3,f4,plots,t_start,t_stop, plotOption,verbose)

  Filters the specified data using a Tukey band stop filter or a Low Pass
  filter (only needs f1,f2)

  -- takes fft of data, applies filter in frequency domain and fft's back
  -----          ------ 
      |\        /|      
      | \      / |      
      |  ------  |      
      |  |    |  |      
     f1 f2   f3 f4       
  Notes:
  
  Takes sample length by looking at difference between adjacent timestamps
  N.B. need to check this as I have found slightly different values
  depending on the timestamps selected. Might be better to take length of
  data/no. samples
 
  For power spectrum, it plots y on log scale and x as linear
  h is the handle to the plot

 Inputs   d: complete data set
          switchdata: =1 for switched data, 0 for combined dataT
          channel: channel number
          filter: Pick type of filter ' Tukey' or 'LowPass'
 (Optional)
          f1,f2,f3,f4 : frequency positions of the tukey filter
           plots: 1 if you want to see the plots
          t_start:  start time (in case you only want to plot a subset of
          the data) 
          t_stop  : stop_time
          plotOption:  string of color, etc.
          verbose : switch = 1 if you want to see verbosity

0001 function d = filter_data(d,switchdata,channel,filter,f1,f2,f3,f4,plots,t_start,t_stop, plotOption,verbose)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function d = filter_data(d,switchdata,channel,filter,f1,f2,f3,f4,plots,t_start,t_stop, plotOption,verbose)
0006 %
0007 %  Filters the specified data using a Tukey band stop filter or a Low Pass
0008 %  filter (only needs f1,f2)
0009 %
0010 %  -- takes fft of data, applies filter in frequency domain and fft's back
0011 %  -----          ------
0012 %      |\        /|
0013 %      | \      / |
0014 %      |  ------  |
0015 %      |  |    |  |
0016 %     f1 f2   f3 f4
0017 %  Notes:
0018 %
0019 %  Takes sample length by looking at difference between adjacent timestamps
0020 %  N.B. need to check this as I have found slightly different values
0021 %  depending on the timestamps selected. Might be better to take length of
0022 %  data/no. samples
0023 %
0024 %  For power spectrum, it plots y on log scale and x as linear
0025 %  h is the handle to the plot
0026 %
0027 % Inputs   d: complete data set
0028 %          switchdata: =1 for switched data, 0 for combined dataT
0029 %          channel: channel number
0030 %          filter: Pick type of filter ' Tukey' or 'LowPass'
0031 % (Optional)
0032 %          f1,f2,f3,f4 : frequency positions of the tukey filter
0033 %           plots: 1 if you want to see the plots
0034 %          t_start:  start time (in case you only want to plot a subset of
0035 %          the data)
0036 %          t_stop  : stop_time
0037 %          plotOption:  string of color, etc.
0038 %          verbose : switch = 1 if you want to see verbosity
0039 
0040 % Outputs  filtered data is written out in either d...switchData or
0041 %            d...dataT
0042 %            (Need to change this when we decide how to use it in the pipeline)
0043 %
0044 %
0045 %  act 3/9/2010
0046 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0047 
0048 
0049 
0050 if(nargin<9)
0051   plotOption = 'b';
0052   t_start = 1;
0053   t_stop = length(d.antenna0.receiver.utc);
0054   plots = 1
0055   verbose = 0
0056 end
0057 
0058 if(nargin<10)
0059   plotOption = 'b';
0060   t_start = 1;
0061   t_stop = length(d.antenna0.receiver.utc);
0062   verbose = 0
0063 end
0064 
0065 if(nargin<12)
0066     plotOption = 'b';
0067     verbose = 0
0068 end
0069 
0070 
0071 if (switchdata == 1)
0072     data = d.antenna0.receiver.switchData(:,channel);
0073 else
0074     data = d.antenna0.receiver.dataT(:,channel);
0075 end
0076     
0077 % First plot the power spectrum of the exisiting data
0078 
0079 t_sec = (d.antenna0.receiver.utc(t_start:t_stop)-d.antenna0.receiver.utc(t_start))*24*60*60;
0080 t_sec = t_sec - t_sec(1);
0081 Fs = 1/(t_sec(3)-t_sec(2));
0082 L= length(data(t_start:t_stop));
0083 NFFT = 2^nextpow2(L);
0084 y_power_spec = (fft(data(t_start:t_stop),NFFT));
0085 x_power_spec = Fs/2*linspace(0,1,NFFT/2+1);
0086 
0087 if (plots==1)
0088 eval(sprintf('h = semilogy(x_power_spec, 2*abs(y_power_spec(1:NFFT/2+1)),''%s'');', plotOption));
0089 title('Power spectrum of original data set')
0090 end 
0091 
0092 % Do different things depending on what filter you want.
0093 
0094 if (filter == 'Tukey  ')
0095     
0096     % Now ask what frequencies you want to filter
0097     if(verbose==1)
0098     disp(' ')
0099     disp('filter_data_choice:: A Tukey band stop filter will be applied')
0100     disp('filter_data_choice::  ')
0101     disp('filter_data_choice::  -----          ------ ')
0102     disp('filter_data_choice::      |\        /|      ')
0103     disp('filter_data_choice::      | \      / |      ')
0104     disp('filter_data_choice::      |  ------  |      ')
0105     disp('filter_data_choice::      |  |    |  |      ')
0106     disp('filter_data_choice::      f1 f2   f3 f4      ') 
0107     disp('filter_data_choice::  ')
0108     end
0109     if (nargin==8)
0110         disp('filter_data_choice:: Using the following frequencies: ')
0111         f1
0112         f2
0113         f3
0114         f4
0115     end
0116     if (nargin<8)
0117         f1 = input('filter_data_choice:: Please enter values for f1: ');
0118         f2 = input('filter_data_choice:: Please enter values for f2: ');
0119         f3 = input('filter_data_choice:: Please enter values for f3: ');
0120         f4 = input('filter_data_choice:: Please enter values for f4: ');
0121     end
0122 
0123     % Now call the tukey_bandstop function to create the filter
0124 
0125     filter = Tukey_filter(f1,f2,f3,f4, Fs, NFFT, x_power_spec);
0126 
0127 
0128 else
0129     if (filter =='LowPass')
0130     disp(' ')
0131     disp('filter_data_choice::  A low pass filter will be applied')
0132     disp('filter_data_choice::  ')
0133     disp('filter_data_choice::  -----         ')
0134     disp('filter_data_choice::      |\        ')
0135     disp('filter_data_choice::      | \       ')
0136     disp('filter_data_choice::      |  -----  ')
0137     disp('filter_data_choice::      |  |      ')
0138     disp('filter_data_choice::      f1 f2     ') 
0139     disp('filter_data_choice::  ')
0140 
0141     if (nargin==6)
0142         disp('filter_data_choice:: Using the following frequencies: ')
0143         f1
0144         f2
0145 
0146     end
0147     if (nargin<6)
0148         f1 = input('filter_data_choice:: Please enter values for f1: ');
0149         f2 = input('filter_data_choice:: Please enter values for f2: ');
0150 
0151     end
0152 
0153     % Now call the lowpass filter function to create the filter
0154 
0155     filter = LowPass_filter(f1,f2, Fs, NFFT, x_power_spec);
0156     end
0157 end
0158 
0159     
0160 
0161 
0162 % Now multiply the y_power_spectrum by the filter
0163 
0164 filtered_fft = y_power_spec.* filter';
0165 
0166 
0167 if (plots==1)
0168     figure
0169     eval(sprintf('h = semilogy(x_power_spec, 2*abs(filtered_fft(1:NFFT/2+1)),''%s'');', plotOption));
0170     title('Power spectrum of the filtered data')
0171 end
0172 
0173 % Now go back to the time domain
0174 
0175 dummy = real(ifft(filtered_fft,NFFT));
0176 dummy = dummy(1:L);
0177 
0178 
0179 % Now plot the original and filtered data
0180 
0181 if (plots ==1)
0182 figure
0183 ax(1) = subplot(3,1,1);
0184 plot(data)
0185 title('Original data stream')
0186 ax(2) = subplot(3,1,2);
0187 plot(dummy)
0188 title('Filtered data stream')
0189 ax(3) = subplot(3,1,3);
0190 plot(data - dummy)
0191 title('Residuals')
0192 linkaxes(ax,'x');
0193 figure
0194 %spectrogram(dummy,1024,16,2048,100)
0195 end
0196 
0197 % Finally write back to the data array d
0198 
0199 if (switchdata==1)
0200     disp('filter_data_choice:: writing back to to d....switchData')
0201     d.antenna0.receiver.switchData(:,channel) = dummy;
0202 else
0203      disp('filter_data_choice:: writing back to d....dataT')
0204     d.antenna0.receiver.dataT(:,channel) = dummy;
0205 end
0206 
0207 clear dummy;
0208