comparison utilities/UTIL_Butterworth.m @ 38:c2204b18f4a2 tip

End nov big change

+function x=UTIL_Butterworth (x, dt, fl, fu, order)
+% UTIL_Butterworth (x, dt, fu, fl, order)
+% Taken from Yuel and beauchamp page 261
+% NB error in their table for K (see their text)
+% x is original signal
+% x is the filtered output (approx 3 dB down at cutoff for first order filter)
+% fu, fl upper and lower cutoff
+% order is the number of times the filter is applied
+% (approx 6 dB attenuation is corrected)
+
+sampleRate=1/dt;
+if 4*fu>sampleRate  
+    error(['UTIL_Butterworth: sample rate ' num2str(sampleRate) ' is too low.  Sampling rate should be >' num2str(4*fu)])
+end
+
+
+q=(pi*dt*(fu-fl));
+J=1/(1+ cot(q));
+K= (2*cos(pi*dt*(fu+fl)))/((1+tan(q))*cos(q));
+L= (tan(q)-1)/(tan(q)+1);
+b=[J 0 -J];
+a=[1 -K  -L];
+
+[numChannels nDataPoints]=size(x);
+for channelNumber=1:numChannels
+	for j=1:order
+		% x2 compensate for 6 dB attenuation
+		x(channelNumber,:)=2*filter(b, a, x(channelNumber,:));
+	end
+end
+