annotate aim-mat/tools/@signal/lowpass.m @ 4:537f939baef0 tip

various bug fixes and changed copyright message
author Stefan Bleeck <bleeck@gmail.com>
date Tue, 16 Aug 2011 14:37:17 +0100
parents 20ada0af3d7d
children
rev   line source
bleeck@3 1 % method of class @signal
bleeck@3 2 %
bleeck@3 3 % INPUT VALUES:
bleeck@3 4 %
bleeck@3 5 % RETURN VALUE:
bleeck@3 6 %
bleeck@3 7 %
bleeck@3 8 % This external file is included as part of the 'aim-mat' distribution package
bleeck@3 9 % (c) 2011, University of Southampton
bleeck@3 10 % Maintained by Stefan Bleeck (bleeck@gmail.com)
bleeck@3 11 % download of current version is on the soundsoftware site:
bleeck@3 12 % http://code.soundsoftware.ac.uk/projects/aimmat
bleeck@3 13 % documentation and everything is on http://www.acousticscale.org
bleeck@3 14
bleeck@3 15
bleeck@3 16 function filtered_sig=lowpass(sig,frequency,stopband,ripple,stopbandatten)
bleeck@3 17 % hack for an phase true lowpassfilter with cutoff at frequency
bleeck@3 18 % used is a ButterworthFilter
bleeck@3 19
bleeck@3 20 if nargin < 5
bleeck@3 21 stopbandatten=60; % in dB - how many dB the signal is reduced in the stopband at least
bleeck@3 22 end
bleeck@3 23 if nargin < 4
bleeck@3 24 ripple=1; % in dB = ripple in the passband
bleeck@3 25 end
bleeck@3 26 if nargin <3
bleeck@3 27 stopband=frequency*2; % eine Oktave drüber
bleeck@3 28 end
bleeck@3 29
bleeck@3 30 nyquist=getsr(sig)/2;
bleeck@3 31 % fre_low=2;
bleeck@3 32 fre_high=frequency;
bleeck@3 33
bleeck@3 34 % Finde raus, wieviel Punkte der Filter dafür haben muss
bleeck@3 35 Wpass=fre_high/nyquist;
bleeck@3 36 Wstop=(fre_high+stopband)/nyquist;
bleeck@3 37 Wstop=min(Wstop,0.999999);
bleeck@3 38 % try
bleeck@3 39 [n,Wn] = buttord(Wpass,Wstop,ripple,stopbandatten);
bleeck@3 40 % Berechne den IIR-Filter
bleeck@3 41 [b,a] = butter(n,Wn);
bleeck@3 42
bleeck@3 43 vals=sig.werte';
bleeck@3 44
bleeck@3 45 % % fill the part behind the signal and in front of the signal with
bleeck@3 46 % % values to avoid corner effects. this is probably not clever in all
bleeck@3 47 % % cases...
bleeck@3 48 % firstval=vals(1);
bleeck@3 49 % lastval=vals(end);
bleeck@3 50 % nr_vals=length(vals);
bleeck@3 51 % vals=[ones(1,nr_vals)*firstval vals ones(1,nr_vals)*lastval];
bleeck@3 52
bleeck@3 53 nvals = filtfilt(b,a,vals);
bleeck@3 54 % extract the values back
bleeck@3 55 % nvals=nvals(nr_vals+1:2*nr_vals);
bleeck@3 56
bleeck@3 57 filtered_sig=sig; % a copy of the old one
bleeck@3 58 newname=sprintf('Lowpass filterd (%3.2fkHz) Signal: %s',frequency/1000,getname(sig));
bleeck@3 59 filtered_sig=setname(filtered_sig,newname);
bleeck@3 60 filtered_sig.werte=nvals';
bleeck@3 61 %
bleeck@3 62 % catch
bleeck@3 63 % disp('error: cant do the low pass filtering');
bleeck@3 64 % filtered_sig=sig;
bleeck@3 65 % end
bleeck@3 66
bleeck@3 67 % figure(235423)
bleeck@3 68 % plot(sig);
bleeck@3 69 % hold on
bleeck@3 70 % plot(filtered_sig,'g');
bleeck@3 71 % s=0;