0001 function d = filter_data(d,switchdata,channel,f1,f2,f3,f4,plots,t_start,t_stop, plotOption,verbose)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
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
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
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
0115
0116 filter = Tukey_filter(f1,f2,f3,f4, Fs, NFFT, x_power_spec);
0117
0118
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
0130
0131 dummy = real(ifft(filtered_fft,NFFT));
0132 dummy = dummy(1:L);
0133
0134
0135
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
0153
0154 if (switchdata==1)
0155
0156 d.antenna0.receiver.switchData(:,channel) = dummy;
0157 else
0158
0159 d.antenna0.receiver.data(:,channel) = dummy;
0160 end
0161
0162 clear dummy;
0163