tomwalters@0: % method of class @signal
tomwalters@0: % 
tomwalters@0: %   INPUT VALUES:
tomwalters@0: %  
tomwalters@0: %   RETURN VALUE:
tomwalters@0: %
tomwalters@0: % 
tomwalters@0: % (c) 2003, University of Cambridge, Medical Research Council 
tomwalters@0: % Stefan Bleeck (stefan@bleeck.de)
tomwalters@0: % http://www.mrc-cbu.cam.ac.uk/cnbh/aimmanual
tomwalters@0: % $Date: 2003/07/17 10:56:16 $
tomwalters@0: % $Revision: 1.2 $
tomwalters@0: 
tomwalters@0: function filtered_sig=bandpass(sig,lowfrequency,highfrequency,stopbandwidth_low,stopbandwidth_high,ripple,stopbandatten)
tomwalters@0: % hack for an phase true bandpassfilter with cutoff at frequency
tomwalters@0: % used is a ButterworthFilter
tomwalters@0: % this is all crap, but it does not work properly otherwise. 
tomwalters@0: 
tomwalters@0: grafix=0;
tomwalters@0: 
tomwalters@0: if nargin < 7
tomwalters@0:     stopbandatten=60; % in dB - how many dB the signal is reduced in the stopband at least
tomwalters@0: end
tomwalters@0: if nargin < 6
tomwalters@0:     ripple=3; % in dB = ripple in the passband
tomwalters@0: end
tomwalters@0: if nargin < 5
tomwalters@0:     stopbandwidth_high=highfrequency*2; % one octave above
tomwalters@0: end
tomwalters@0: if nargin < 4
tomwalters@0:     stopbandwidth_low=lowfrequency/2; % one octave above
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % I dont understand enough to make it work with a butterwort filer.
tomwalters@0: % easy solution: first perform a lowpass, then a highpass filter...
tomwalters@0: % sorry, Thanks for any help!
tomwalters@0: if grafix
tomwalters@0: 	figure(23534)
tomwalters@0: 	plot(powerspectrum(sig),[100,getsr(sig)/2],'b')
tomwalters@0: 	hold on
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if highfrequency>0
tomwalters@0:     filtered_sig_low=lowpass(sig,highfrequency,highfrequency+stopbandwidth_high,ripple,stopbandatten);
tomwalters@0: else
tomwalters@0:     filtered_sig_low=sig;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if grafix
tomwalters@0: 	figure(23534)
tomwalters@0: 	plot(powerspectrum(filtered_sig_low),[100,getsr(filtered_sig_low)/2],'r')
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: if lowfrequency>0
tomwalters@0:     filtered_sig=highpass(filtered_sig_low,lowfrequency,lowfrequency-stopbandwidth_low,ripple,stopbandatten);
tomwalters@0: if grafix
tomwalters@0: 	figure(23534)
tomwalters@0: 	plot(powerspectrum(filtered_sig),[100,getsr(filtered_sig_low)/2],'g')
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: return
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % if nargin < 4
tomwalters@0: %     ripple=1; % in dB = ripple in the passband
tomwalters@0: % end
tomwalters@0: % if nargin <3
tomwalters@0: %     stopband=frequency*2; % eine Oktave drüber
tomwalters@0: % end
tomwalters@0: 
tomwalters@0: ripple=3;
tomwalters@0: 
tomwalters@0: nyquist=getsr(sig)/2;
tomwalters@0: % fre_low=2;
tomwalters@0: % fre_high=frequency;
tomwalters@0: 
tomwalters@0: Wpass=[lowfrequency highfrequency]/nyquist;
tomwalters@0: Wstop=[lowfrequency-stopbandwidth highfrequency+stopbandwidth]/nyquist;
tomwalters@0: [n,Wn] = buttord(Wpass,Wstop,ripple,stopbandatten);
tomwalters@0: % [n,Wn] = buttord(fre_high/nyquist,(fre_high+stopband)/nyquist,ripple,stopbandatten);
tomwalters@0: % Berechne den IIR-Filter
tomwalters@0: [b,a] = butter(n,Wn);
tomwalters@0: 
tomwalters@0: % zum testen:
tomwalters@0: freqz(b,a,128,getsr(sig))
tomwalters@0: 
tomwalters@0: vals=sig.werte';
tomwalters@0: 
tomwalters@0: % fill the part behind the signal and in front of the signal with
tomwalters@0: % values
tomwalters@0: % firstval=vals(1);
tomwalters@0: % lastval=vals(end);
tomwalters@0: % nr_vals=length(vals);
tomwalters@0: % vals=[ones(1,nr_vals)*firstval vals ones(1,nr_vals)*lastval];
tomwalters@0: 
tomwalters@0: plot(vals)
tomwalters@0: hold on
tomwalters@0: 
tomwalters@0: nvals = filter(b,a,vals);
tomwalters@0: % nvals = filtfilt(b,a,vals);
tomwalters@0: % extract the values back 
tomwalters@0: % nvals=nvals(nr_vals+1:2*nr_vals);
tomwalters@0: 
tomwalters@0: plot(nvals,'r')
tomwalters@0: 
tomwalters@0: % filtered_sig=0;
tomwalters@0: % return
tomwalters@0: 
tomwalters@0: filtered_sig=sig;	% a copy of the old one
tomwalters@0: newname=sprintf('Bandpass filterd (%3.2fkHz - %3.2fkHz) Signal: %s',lowfrequency/1000,highfrequency/1000,getname(sig));
tomwalters@0: filtered_sig=setname(filtered_sig,newname);
tomwalters@0: filtered_sig.werte=nvals';
tomwalters@0: 
tomwalters@0: % figure(235423)
tomwalters@0: % plot(sig);
tomwalters@0: % hold on
tomwalters@0: % plot(filtered_sig,'g');
tomwalters@0: % s=0;