Home > reduc > filter_data.m

filter_data

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

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

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

  Filters the specified data using a Tukey band stop filter

  -- 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 data
          channel: channel number
 (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: =0 if you don't want to see the ascii art
 
 Outputs  filtered data is written out in either d...switchData or
            d...data
            (Need to change this when we decide how to use it in the pipeline)
           

  act 3/9/2010
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Sun 14-Jun-2015 17:12:45 by m2html © 2005