bleeck@3: % method of class @signal
bleeck@3: % 
bleeck@3: %   INPUT VALUES:
bleeck@3: %  
bleeck@3: %   RETURN VALUE:
bleeck@3: %
bleeck@3: % 
bleeck@3: % This external file is included as part of the 'aim-mat' distribution package
bleeck@3: % (c) 2011, University of Southampton
bleeck@3: % Maintained by Stefan Bleeck (bleeck@gmail.com)
bleeck@3: % download of current version is on the soundsoftware site: 
bleeck@3: % http://code.soundsoftware.ac.uk/projects/aimmat
bleeck@3: % documentation and everything is on http://www.acousticscale.org
bleeck@3: 
bleeck@3: function filtered_sig=bandpass(sig,lowfrequency,highfrequency,stopbandwidth_low,stopbandwidth_high,ripple,stopbandatten)
bleeck@3: % hack for an phase true bandpassfilter with cutoff at frequency
bleeck@3: % used is a ButterworthFilter
bleeck@3: % this is all crap, but it does not work properly otherwise. 
bleeck@3: 
bleeck@3: grafix=0;
bleeck@3: 
bleeck@3: if nargin < 7
bleeck@3:     stopbandatten=60; % in dB - how many dB the signal is reduced in the stopband at least
bleeck@3: end
bleeck@3: if nargin < 6
bleeck@3:     ripple=3; % in dB = ripple in the passband
bleeck@3: end
bleeck@3: if nargin < 5
bleeck@3:     stopbandwidth_high=highfrequency*2; % one octave above
bleeck@3: end
bleeck@3: if nargin < 4
bleeck@3:     stopbandwidth_low=lowfrequency/2; % one octave above
bleeck@3: end
bleeck@3: 
bleeck@3: % I dont understand enough to make it work with a butterwort filer.
bleeck@3: % easy solution: first perform a lowpass, then a highpass filter...
bleeck@3: % sorry, Thanks for any help!
bleeck@3: if grafix
bleeck@3: 	figure(235340);
bleeck@3:     clf 
bleeck@3: 	plot(powerspectrum(sig),[100,getsr(sig)/2],'b')
bleeck@3: 	hold on
bleeck@3: end
bleeck@3: 
bleeck@3: if highfrequency>0 && highfrequency<getsr(sig)/2
bleeck@3:     filtered_sig_low=lowpass(sig,highfrequency,highfrequency+stopbandwidth_high,ripple,stopbandatten);
bleeck@3: else
bleeck@3:     filtered_sig_low=sig;
bleeck@3: end
bleeck@3: 
bleeck@3: if grafix
bleeck@3: 	figure(235340)
bleeck@3: 	plot(powerspectrum(filtered_sig_low),[100,getsr(filtered_sig_low)/2],'r')
bleeck@3: end
bleeck@3: 
bleeck@3: 
bleeck@3: if lowfrequency>0 && lowfrequency<getsr(sig)/2
bleeck@3:     filtered_sig=highpass(filtered_sig_low,lowfrequency,lowfrequency-stopbandwidth_low,ripple,stopbandatten);
bleeck@3: else
bleeck@3:     filtered_sig=filtered_sig_low;
bleeck@3: end
bleeck@3: if grafix
bleeck@3: 	figure(235340)
bleeck@3: 	plot(powerspectrum(filtered_sig),[100,getsr(filtered_sig_low)/2],'g')
bleeck@3: end
bleeck@3: 
bleeck@3: 
bleeck@3: return
bleeck@3: 
bleeck@3: 
bleeck@3: % if nargin < 4
bleeck@3: %     ripple=1; % in dB = ripple in the passband
bleeck@3: % end
bleeck@3: % if nargin <3
bleeck@3: %     stopband=frequency*2; % eine Oktave drüber
bleeck@3: % end
bleeck@3: 
bleeck@3: ripple=3;
bleeck@3: 
bleeck@3: nyquist=getsr(sig)/2;
bleeck@3: % fre_low=2;
bleeck@3: % fre_high=frequency;
bleeck@3: 
bleeck@3: Wpass=[lowfrequency highfrequency]/nyquist;
bleeck@3: Wstop=[lowfrequency-stopbandwidth highfrequency+stopbandwidth]/nyquist;
bleeck@3: [n,Wn] = buttord(Wpass,Wstop,ripple,stopbandatten);
bleeck@3: % [n,Wn] = buttord(fre_high/nyquist,(fre_high+stopband)/nyquist,ripple,stopbandatten);
bleeck@3: % Berechne den IIR-Filter
bleeck@3: [b,a] = butter(n,Wn);
bleeck@3: 
bleeck@3: % zum testen:
bleeck@3: freqz(b,a,128,getsr(sig))
bleeck@3: 
bleeck@3: vals=sig.werte';
bleeck@3: 
bleeck@3: % fill the part behind the signal and in front of the signal with
bleeck@3: % values
bleeck@3: % firstval=vals(1);
bleeck@3: % lastval=vals(end);
bleeck@3: % nr_vals=length(vals);
bleeck@3: % vals=[ones(1,nr_vals)*firstval vals ones(1,nr_vals)*lastval];
bleeck@3: 
bleeck@3: plot(vals)
bleeck@3: hold on
bleeck@3: 
bleeck@3: nvals = filter(b,a,vals);
bleeck@3: % nvals = filtfilt(b,a,vals);
bleeck@3: % extract the values back 
bleeck@3: % nvals=nvals(nr_vals+1:2*nr_vals);
bleeck@3: 
bleeck@3: plot(nvals,'r')
bleeck@3: 
bleeck@3: % filtered_sig=0;
bleeck@3: % return
bleeck@3: 
bleeck@3: filtered_sig=sig;	% a copy of the old one
bleeck@3: newname=sprintf('Bandpass filterd (%3.2fkHz - %3.2fkHz) Signal: %s',lowfrequency/1000,highfrequency/1000,getname(sig));
bleeck@3: filtered_sig=setname(filtered_sig,newname);
bleeck@3: filtered_sig.werte=nvals';
bleeck@3: 
bleeck@3: % figure(235423)
bleeck@3: % plot(sig);
bleeck@3: % hold on
bleeck@3: % plot(filtered_sig,'g');
bleeck@3: % s=0;