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: 
bleeck@3: function filtered_sig=lowpass(sig,frequency,stopband,ripple,stopbandatten)
bleeck@3: % hack for an phase true lowpassfilter with cutoff at frequency
bleeck@3: % used is a ButterworthFilter
bleeck@3: 
bleeck@3: if nargin < 5
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 < 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: nyquist=getsr(sig)/2;
bleeck@3: % fre_low=2;
bleeck@3: fre_high=frequency;
bleeck@3: 
bleeck@3: % Finde raus, wieviel Punkte der Filter dafür haben muss
bleeck@3: Wpass=fre_high/nyquist;
bleeck@3: Wstop=(fre_high+stopband)/nyquist;
bleeck@3: Wstop=min(Wstop,0.999999);
bleeck@3: % try
bleeck@3:     [n,Wn] = buttord(Wpass,Wstop,ripple,stopbandatten);
bleeck@3:     % Berechne den IIR-Filter
bleeck@3:     [b,a] = butter(n,Wn);
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 to avoid corner effects. this is probably not clever in all
bleeck@3: %     % cases...
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:     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:     filtered_sig=sig;	% a copy of the old one
bleeck@3:     newname=sprintf('Lowpass filterd (%3.2fkHz) Signal: %s',frequency/1000,getname(sig));
bleeck@3:     filtered_sig=setname(filtered_sig,newname);
bleeck@3:     filtered_sig.werte=nvals';
bleeck@3: % 
bleeck@3: % catch
bleeck@3: %     disp('error: cant do the low pass filtering');
bleeck@3: %     filtered_sig=sig;
bleeck@3: % end
bleeck@3: 
bleeck@3: % figure(235423)
bleeck@3: % plot(sig);
bleeck@3: % hold on
bleeck@3: % plot(filtered_sig,'g');
bleeck@3: % s=0;