annotate utilities/UTIL_Butterworth.m @ 38:c2204b18f4a2 tip

End nov big change
author Ray Meddis <rmeddis@essex.ac.uk>
date Mon, 28 Nov 2011 13:34:28 +0000
parents
children
rev   line source
rmeddis@38 1 function x=UTIL_Butterworth (x, dt, fl, fu, order)
rmeddis@38 2 % UTIL_Butterworth (x, dt, fu, fl, order)
rmeddis@38 3 % Taken from Yuel and beauchamp page 261
rmeddis@38 4 % NB error in their table for K (see their text)
rmeddis@38 5 % x is original signal
rmeddis@38 6 % x is the filtered output (approx 3 dB down at cutoff for first order filter)
rmeddis@38 7 % fu, fl upper and lower cutoff
rmeddis@38 8 % order is the number of times the filter is applied
rmeddis@38 9 % (approx 6 dB attenuation is corrected)
rmeddis@38 10
rmeddis@38 11 sampleRate=1/dt;
rmeddis@38 12 if 4*fu>sampleRate
rmeddis@38 13 error(['UTIL_Butterworth: sample rate ' num2str(sampleRate) ' is too low. Sampling rate should be >' num2str(4*fu)])
rmeddis@38 14 end
rmeddis@38 15
rmeddis@38 16
rmeddis@38 17 q=(pi*dt*(fu-fl));
rmeddis@38 18 J=1/(1+ cot(q));
rmeddis@38 19 K= (2*cos(pi*dt*(fu+fl)))/((1+tan(q))*cos(q));
rmeddis@38 20 L= (tan(q)-1)/(tan(q)+1);
rmeddis@38 21 b=[J 0 -J];
rmeddis@38 22 a=[1 -K -L];
rmeddis@38 23
rmeddis@38 24 [numChannels nDataPoints]=size(x);
rmeddis@38 25 for channelNumber=1:numChannels
rmeddis@38 26 for j=1:order
rmeddis@38 27 % x2 compensate for 6 dB attenuation
rmeddis@38 28 x(channelNumber,:)=2*filter(b, a, x(channelNumber,:));
rmeddis@38 29 end
rmeddis@38 30 end
rmeddis@38 31