% Filter Tool (L Braile, October 2001; revised January 2002)
%
% Filters (LP, HP, BP or Notch) using a Butterworth filter.  Filtering is performed
% using the filtfilt routine that results in a zero phase filter.  The order
% of the resulting filter is 6 ("6-pole" filter; notch filter is of order 16).  
% Requires Matlab "Signal 
% Processing Toolbox".
%
% define signals each 201 points long, dt=0.01
dt=0.01;
t=0:dt:2;
%
% b1 short boxcar
b1=[zeros(1,75) ones(1,51) zeros(1,75)];
% b2 long boxcar
b2=[zeros(1,50) ones(1,101) zeros(1,50)];
% sq square wave
sq=[zeros(1,75) ones(1,25) -1*ones(1,25) zeros(1,76)];
% sp spike
sp=[zeros(1,100) 1 zeros(1,100)];
% st step function
st=[zeros(1,50) ones(1,151)];
% ra random number
randn('state',1);
ra=randn(1,201);
ramax=max(abs(ra));
ra=ra/ramax;
% un uniformly distributed random number
un=rand(1,201);
un=2*un - 1;
% c1 low frequency cosine
%c1=cos(2*pi*1*t);
% Sum of 3, 8, 25 Hz sins and noise
sn=sin(2*pi*3*t)+sin(2*pi*8*t)+sin(2*pi*20*t)+ra;
smax=max(abs(sn));
sn=sn/smax;
% c2 high frequency cosine
%c2=cos(2*pi*5*t);
% bn boxcar + random noise (ra)
bn=b2+0.5*ra;
bmax=max(abs(bn));
bn=bn/bmax;
% si sin x/x or sinc function
%sit=sin(2*pi*4*t(2:101))./(2*pi*4*t(2:101));
%si=[fliplr(sit) 1 sit];
% bs boxcar + 20 Hx sin
bs=b2+0.2*sin(2*pi*20*t);
bmax=max(abs(bs));
bs=bs/bmax;
% df differentiation operator
%df=[zeros(1,99) 1 -1 zeros(1,100)];
% ts two sin (1 + 20 Hz)
ts=sin(2*pi*t) + sin(2*pi*20*t);
tmax=max(abs(ts));
ts=ts/tmax;
% t1 short triangle
tr1=[1:1:26]/26;
tr2=[25:-1:1]/26;
t1=[zeros(1,75) tr1 tr2 zeros(1,75)];
% t2 long triangle
%tr1=[1:1:51]/51;
%tr2=[50:-1:1]/51;
%t2=[zeros(1,50) tr1 tr2 zeros(1,50)];
% ga gaussian
ga1=exp(-(8*t(2:101).*t(2:101)));
ga=[fliplr(ga1) 1 ga1];
% gn gaussian + noise
gn=ga+0.5*ra;
gmax=max(abs(gn));
gn=gn/gmax;
% cb cosine bell
nsw=100;
   for i = 1:nsw
      j = nsw-i;
      fa1(i) = 0.5*(cos(pi*j/nsw) + 1);
   end
cb=[fa1 1 fliplr(fa1)];
% ex exponential
thalf=0.25;
x1=t(2:151)*0.696/thalf;
ex1=exp(-x1);
ex=[zeros(1,50) 1 ex1];
% w1 Ricker wavelet
[rk]=ricker(201,1,dt,0,25);
rmax=max(abs(rk));
w1=rk/rmax;
% w2 Ricker wavelet
[rk]=ricker(201,5,dt,0,25);
rmax=max(abs(rk));
w2=rk/rmax;
% s1 spike series 1
s1=[zeros(1,10) 1 zeros(1,80) -0.7 zeros(1,100) 0.4 zeros(1,8)];
% s2 spike series 2
s2=[zeros(1,10) 1 zeros(1,45) .4 zeros(1,34) -0.7 zeros(1,55) ...
      -.3 zeros(1,44) 0.4 zeros(1,8)];
% s3 spike series 3
s3=[zeros(1,50) 1 zeros(1,8) -0.2 zeros(1,9) 0.7 zeros(1,8) -0.8 zeros(1,5) 0.3...
      zeros(1,12) -0.4 zeros(1,8) 0.25 zeros(1,5) 0.1 zeros(1,9) -0.1 zeros(1,6)...
      0.1 zeros(1,71)];
% rp ramp
rp1=[1:1:101]/101;
rp=[zeros(1,50) rp1 zeros(1,50)];
% sw sweep 1-5 Hz, 1000 points, dt=0.01
for i=1:183
   ti=(i-1)*dt;
   f=2 + 1*(ti);
   sw1(i)=sin(2*pi*f*ti);
end
sw=[zeros(1,9) sw1 zeros(1,9)];
% di Dirac comb
%di=[zeros(1,10) 1 zeros(1,19) 1 zeros(1,19) 1 zeros(1,19) 1 zeros(1,19) 1 ...
%  zeros(1,19) 1 zeros(1,19) 1 zeros(1,19) 1 zeros(1,19) 1 zeros(1,19) 1 ...
%  zeros(1,10)];
% gs gaussian + sin
gs=ga+0.5*sin(2*pi*20*t);
gmax=max(abs(gs));
gs=gs/gmax;
% w3 wavelet - signal n=2
tm=0.21; n=2; npts=201; na=25; thalf=0;
[w3]=signal(tm,dt,n,0,npts,na,thalf);
% w4 wavelet - signal n=4
%tm=0.41; n=4; npts=201; na=25; thalf=0;
%[w4]=signal(tm,dt,n,0,npts,na,thalf);
%wn wavelet + noise
wn=w2+0.5*ra;
wmax=max(abs(wn));
wn=wn/wmax;
% w5 wavelet - signal n=5
tm=0.51; n=5; npts=201; na=25; thalf=0;
[w5]=signal(tm,dt,n,0,npts,na,thalf);
% w6 wavelet - signal n=7
tm=0.71; n=7; npts=201; na=25; thalf=0;
[w6]=signal(tm,dt,n,0,npts,na,thalf);
% w7 wavelet - signal n=7 thalf=0.15
%tm=0.71; n=7; npts=201; na=25; thalf=0.15;
%[w7]=signal(tm,dt,n,0,npts,na,thalf);
% hu humps + noise
x1=-0.4:0.01:1.6;
hu=0.3*ra+humps(x1)/100;
hmax=max(abs(hu));
hu=hu/hmax;
% w8 wavelet - signal n=9 thalf=0.2
tm=0.91; n=9; npts=201; na=25; thalf=0.2;
[w8]=signal(tm,dt,n,0,npts,na,thalf);
%
%
% plot signals
figure
subplot(5,6,1), plot(t,b1,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'b1, boxcar','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
%
subplot(5,6,2), plot(t,b2,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'b2, boxcar','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,3), plot(t,sq,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
title('                                     Functions and Signals for Filtering Tool',...
   'fontweight','bold','fontsize',12)
text(0.1,-0.9,'sq, sq. wave','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,4), plot(t,sp,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'sp, spike','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,5), plot(t,bs,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'bs, b2 + 20 Hz sin','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,6), plot(t,st,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'st, step func.','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,7), plot(t,ra,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'ra, random nor.','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
%
subplot(5,6,8), plot(t,un,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'un, random unif.','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,9), plot(t,sn,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'sn, 3,8,20 Hz + ra','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,10), plot(t,bn,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'bn, b2 + noise','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,11), plot(t,ts,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'ts, 1&20 Hz sin','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,12), plot(t,t1,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'t1, triangle','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,13), plot(t,gn,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'gn, ga + noise','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
ylabel('Relative Amplitude','FontWeight','bold')
%
subplot(5,6,14), plot(t,ga,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'ga, gaussian','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,15), plot(t,cb,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'cb, cosine bell','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,16), plot(t,ex,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'ex, exponential','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,17), plot(t,rp,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'rp, ramp','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,18), plot(t,gs,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'gs, ga + 20Hz sin','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,19), plot(t,s1,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'s1, spike series','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
%
subplot(5,6,20), plot(t,s2,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'s2, spike series','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,21), plot(t,s3,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'s3, spike series','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,22), plot(t,sw,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'sw, sweep','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,23), plot(t,w1,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'w1, ricker wavelet','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,24), plot(t,w2,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'w2, ricker wavelet','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,25), plot(t,w3,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'w3, wavelet','FontWeight','bold')
xlabel('Time (s)','FontWeight','bold')
%
subplot(5,6,26), plot(t,wn,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'wn, w2 + noise','FontWeight','bold')
set(gca,'YTickLabelMode','manual','YTickLabel',[])
xlabel('Time (s)','FontWeight','bold')
%
subplot(5,6,27), plot(t,w5,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'w5, wavelet','FontWeight','bold')
set(gca,'YTickLabelMode','manual','YTickLabel',[])
xlabel('Time (s)','FontWeight','bold')
%
subplot(5,6,28), plot(t,w6,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'w6, wavelet','FontWeight','bold')
set(gca,'YTickLabelMode','manual','YTickLabel',[])
xlabel('Time (s)','FontWeight','bold')
%
subplot(5,6,29), plot(t,hu,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'hu, humps+noise','FontWeight','bold')
set(gca,'YTickLabelMode','manual','YTickLabel',[])
xlabel('Time (s)','FontWeight','bold')
%
subplot(5,6,30), plot(t,w8,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'w8, wavelet','FontWeight','bold')
set(gca,'YTickLabelMode','manual','YTickLabel',[])
xlabel('Time (s)','FontWeight','bold')
%
%
pause
%

fnq=50;
%
%
while 1
%
% Select signal to be filtered
sig=input('Type 2 character name for signal to be filtered or qq to quit ');
%
%
  if strcmp(sig,'qq') | strcmp(sig,'QQ'), break, end
%
% Select filter type
n1=input('Type 1 for LP, 2 for HP, 3 for BP, 4 for Notch filtering ','s');
nc=str2num(n1);
%
if nc==1
% Low Pass Filtering
   fcut=input('Low Pass Filter; Type filter cutoff (0<fc<50 Hz) ','s');
   fc=str2num(fcut);
   [b,a]=butter(3,fc/fnq,'low');
   y=filtfilt(b,a,sig);
end
%
if nc==2
% High Pass Filtering
   fcut=input('High Pass Filter; Type filter cutoff (0<fc<50 Hz) ','s');
   fc=str2num(fcut);
   [b,a]=butter(3,fc/fnq,'high');
   y=filtfilt(b,a,sig);
end
%
if nc==3
% Band Pass Filtering
   fcut1=input('Band Pass Filter; Type first filter cutoff, f1 (0<f1<f2<50 Hz) ','s');
   fcut2=input('Band Pass Filter; Type second filter cutoff, f2 (0<f1<f2<50 Hz) ','s');
   f1=str2num(fcut1);
   f2=str2num(fcut2);
   [b,a]=butter(3,f1/fnq,'high');
   x=filtfilt(b,a,sig);
   [b,a]=butter(3,f2/fnq,'low');
   y=filtfilt(b,a,x);
end
%
if nc==4
% Notch Filtering
   fcut1=input('Notch Filter; Type first filter cutoff, f1 (0<f1<f2<50 Hz) ','s');
   fcut2=input('Notch Filter; Type second filter cutoff, f2 (0<f1<f2<50 Hz) ','s');
   f1=str2num(fcut1);
   f2=str2num(fcut2);
   [b,a]=butter(8,f1/fnq,'low');
   x1=filtfilt(b,a,sig);
   [b,a]=butter(8,f2/fnq,'high');
   x2=filtfilt(b,a,sig);
   y=x1+x2;
end
%
% Plot input signal and filtered output
figure
subplot(2,2,1), plot(t,sig,'r-','linewidth', 2), axis([0 2 -1.1 1.1])
title('Filtering (LP, HP, BP or Notch) Tool -- Butterworth Filter',...
   'FontWeight','bold','fontsize',12)
text(0.1,-0.9,'Input Signal','FontWeight','bold')
xlabel('Time (s)','FontWeight','bold')
ylabel('Relative Amplitude','FontWeight','bold')
subplot(2,2,3), plot(t,y,'r-','linewidth', 2), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'Filtered Signal','FontWeight','bold')
xlabel('Time (s)','FontWeight','bold')
ylabel('Relative Amplitude','FontWeight','bold')
%
end