%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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