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