octave_packages/tsa-4.2.4/invfdemo.m

   1 % invfdemo      demonstrates Inverse Filtering
   2
   3 %       $Id: invfdemo.m 5090 2008-06-05 08:12:04Z schloegl $\r%  Copyright (C) 1997-2002,2008 by Alois Schloegl <a.schloegl@ieee.org>
   4 %
   5 %    This program is free software: you can redistribute it and/or modify
   6 %    it under the terms of the GNU General Public License as published by
   7 %    the Free Software Foundation, either version 3 of the License, or
   8 %    (at your option) any later version.
   9 %
  10 %    This program is distributed in the hope that it will be useful,
  11 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 %    GNU General Public License for more details.
  14 %
  15 %    You should have received a copy of the GNU General Public License
  16 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  17
  18 load eeg8s.mat;     % load signal
  19 ly=length(eeg8s);
  20 Fs=128;
  21 %a=earpyw(eeg8s',11); % Calculates AR(11) parameters with Yule-Walker method
  22 a=lattice(eeg8s',11); % Calculates AR(11) parameters with Yule-Walker method
  23                     % The AR parameters are the weight taps of IIR Filter
  24 isig=filter([1 -a],1,eeg8s);   % Inverse filtering
  25
  26
  27 subplot(221);
  28 plot((1:ly)/Fs,eeg8s);
  29 title('Signal');
  30 xlabel('t [sec]')
  31
  32 subplot(223);
  33 plot((1:ly)/Fs,isig);
  34 xlabel('t [sec]')
  35 title('Inverse filtered process');
  36
  37 subplot(222);
  38 H=abs(fft(eeg8s,128)/ly).^2;
  39 plot(1:Fs,H);
  40 %plot([H mean(H)*ones(Fs,1)]);
  41 ylabel('S(f)')
  42 xlabel('f [Hz]')
  43 title('Spectrum of original signal');
  44
  45 subplot(224);
  46 H=abs(fft(isig,128)/ly).^2;
  47 plot(1:Fs,H);
  48 %plot([H mean(H)*ones(Fs,1)]);
  49 ylabel('S(f)')
  50 xlabel('f [Hz]')
  51 title('Spectrum of inverse filtered signal');