Butterworth Filters

This set of functions is simply four built-in Matlab functions, repackaged for ease of use (Signal Processing Toolbox is required). If you don't want to go through the rigmarole of designing and implementing a filter with normalized frequencies and so forth every time you filter a signal, this package may be for you. If you are a Matlab pro and an expert in digital signal processing, you will probably not be impressed.

Each function takes the form [filtered_signal,filtb,filta] = bandstop_butterworth(inputsignal,cutoff_freqs,Fs,order)

INPUTS:
inputsignal = input time series
cutoff_freqs = filter corner frequencies in the form [f1 f2]
Fs = data sampling frequency
order = order of Butterworth filter

OUTPUTS:
filtered_signal = the filtered time series
filtb, filta = filter numerator and denominator (optional)

EXAMPLE 1:
load train
t = (1:length(y))/Fs;
y_filt = bandstop_butterworth(y,[800 1000],Fs,4); % cut off between 800 Hz and 1000 Hz

figure
plot(t,y,'b',t,y_filt,'r')
xlabel('time in seconds')
box off
legend('unfiltered','filtered')
sound(y,Fs) % play original time series
pause(2) % pause two seconds
sound(y_filt,Fs) % play filtered time series

EXAMPLE 2:
load train
t = (1:length(y))/Fs;
[y_filt,filtb,filta] =bandstop_butterworth(y,[800 1000],Fs,4); % cut off between 800 Hz and 1000 Hz
[h1,f1] = freqz(filtb,filta,256,Fs);

figure
subplot(3,1,1)
plot(t,y,'b',t,y_filt,'r')
xlabel('time in seconds')
box off
text(0,.1,' time series','units','normalized')

subplot(3,1,2)
AX = plotyy(f1,10*log10(abs(h1)),f1,angle(h1),'semilogx');
set(get(AX(1),'ylabel'),'string','gain (dB)')
set(get(AX(2),'ylabel'),'string','phase (rad)')
xlim(AX(1),[min(f1) max(f1)])
xlim(AX(2),[min(f1) max(f1)])
text(0,.1,' filter response','units','normalized')
box off

[Pxx,f] = pwelch(y,512,256,[],Fs,'onesided');
[Pxxf,f_f]= pwelch(y_filt,512,256,[],Fs,'onesided');
subplot(3,1,3)
semilogx(f,10*log10(Pxx))
hold on
semilogx(f_f,10*log10(Pxxf),'r')
xlabel('frequency (Hz)')
ylabel('PSD (dB)')
xlim([min(f1) max(f1)])
box off
legend('unfiltered','filtered','location','northwest')
legend boxoff