freq_filt

FREQ_FILT: frequency filter data
 vfilt = freq_filt(vin, fresp, FR, dim)
 vin = matrix of data values
     = data structure (with field vin.meas)
     = image structure (with field vin.elem_data)
 fresp = function of freq filter
      e.g. fresp = @(f) f<10;  % values in Hz
 FR  = frame rate (or use FR parameter on data structures)
    or FR is a structure with fields
     .FR = Frame Rate
     .padding = length of padding (in seconds)
 dim = dimension along which to filter (default is 2)

0001 function vfilt = freq_filt(vin, fresp, FR, dim)
0002 %FREQ_FILT: frequency filter data
0003 % vfilt = freq_filt(vin, fresp, FR, dim)
0004 % vin = matrix of data values
0005 %     = data structure (with field vin.meas)
0006 %     = image structure (with field vin.elem_data)
0007 % fresp = function of freq filter
0008 %      e.g. fresp = @(f) f<10;  % values in Hz
0009 % FR  = frame rate (or use FR parameter on data structures)
0010 %    or FR is a structure with fields
0011 %     .FR = Frame Rate
0012 %     .padding = length of padding (in seconds)
0013 % dim = dimension along which to filter (default is 2)
0014 
0015 % (C) Andy Adler 2019. License: GPL v2 or v3.
0016 % $Id: freq_filt.m 6066 2020-03-12 23:52:06Z aadler $
0017 
0018 % if FR not given, see if it's a parameter
0019 
0020 if ischar(vin) && strcmp(vin,'UNIT_TEST'); do_unit_test; return; end
0021 
0022 p.padding = 1; % seconds;
0023 if nargin<3; 
0024    p.FR = vin.FR;
0025 else
0026    if isnumeric(FR)
0027      p.FR = FR;
0028    elseif isstruct(FR);
0029      for ff= fieldnames(FR)'; % wish matlab could do this easily
0030        p.(ff{1}) = FR.(ff{1});
0031      end
0032    else
0033      error('Don''t understand parameter FR');
0034    end
0035      
0036       
0037 end
0038 if nargin>=4
0039    p.dim = dim;
0040 elseif ~isfield(p,'dim');
0041    p.dim = 2;
0042 end
0043 
0044 if isstruct(vin) && isfield(vin,'type');
0045    switch vin.type
0046      case 'data'
0047         vin.meas = do_freq_filt(vin.meas, fresp, p);
0048         vfilt = vin;
0049      case 'image'
0050         vin.elem_data = do_freq_filt(vin.elem_data, fresp, p);
0051         vfilt = vin;
0052      otherwise
0053         error('Can''t process object of type %',vin.type);
0054    end
0055 elseif isnumeric(vin)
0056    vfilt = do_freq_filt(vin, fresp, p);
0057 else 
0058    error('can''t process object');
0059 end
0060       
0061 
0062 
0063 % Filter in the frequency direction
0064 function s = do_freq_filt(s,fresp, p)
0065   rs = isreal(s);
0066   s = padsignal(s,p);
0067   f = fft(s,[],p.dim);
0068   fax = freq_axis_filt(p.FR,size(s,p.dim),fresp);
0069   f = f .* fshape(p,fax);
0070   s= ifft(f,[],p.dim);
0071   if rs % is signal is real
0072      if norm(imag(s(:))) > 1e-11
0073         error('FFT filter has imag output');
0074      end
0075      s = real(s);
0076   end
0077   s = unpadsignal(s,p);
0078 
0079 function s = fshape(p,s);
0080   fsh = ones(1,length(size(s)));
0081   fsh(p.dim) = prod(size(s));
0082   s= reshape(s,fsh);
0083 
0084 function plen = pad_len(p);
0085   plen = round(p.FR * p.padding);
0086 
0087 % Add padding connecting the last to the first sample
0088 function s = padsignal(s,p);
0089   lsup = linspace(0,1,pad_len(p)); lsup=fshape(p,lsup);
0090   lsdn = linspace(1,0,pad_len(p)); lsdn=fshape(p,lsdn);
0091   switch p.dim % wish I knew how to generalize in m*lab
0092     case 1; pad = s(1,:,:).*lsdn + s(end,:,:).*lsup;
0093     case 2; pad = s(:,1,:).*lsdn + s(:,end,:).*lsup;
0094     case 3; pad = s(:,:,1).*lsdn + s(:,:,end).*lsup;
0095     otherwise; error('Problem with dim above 3');
0096   end
0097   
0098   s = cat(p.dim,s,pad);
0099 
0100 % Add padding connecting the last to the first sample
0101 function s = unpadsignal(s,p);
0102   plen = round(p.FR * p.padding);
0103   switch p.dim % wish I knew how to generalize in m*lab
0104     case 1; s(end+1-(1:plen),:,:) = [];
0105     case 2; s(:,end+1-(1:plen),:) = [];
0106     case 3; s(:,:,end+1-(1:plen)) = [];
0107     otherwise; error('Problem with dim above 3');
0108   end
0109   
0110 function fax = freq_axis_filt( FR, lD, fresp);
0111   fax = linspace(0,FR,lD+1);
0112   fax(end)=[];
0113   fax(fax>FR/2) = fax(fax>FR/2) - FR;
0114   fax = feval(fresp, abs(fax));
0115 
0116 function do_unit_test
0117   clf; subplot(3,1,1);
0118   FR = 100; t = (0:1.1e3)/FR;
0119   s = sin(2*pi*5*t) + 2*cos(2*pi*15*t) +  3*sin(2*pi*0.1*t);
0120   plot(t,s); hold on;
0121 
0122   for fnum = 1:4; switch fnum
0123      case 1; fresp = @(f) f<10;
0124              sf= freq_filt(s,fresp, FR);
0125              sf12=[1.879618111164496, 1.275026363303550];
0126 
0127      case 2; fresp = @(f) (f<1) + (f>=1)./(f+eps);
0128              p.FR = FR; p.padding = 1;
0129              sf= freq_filt(s,fresp, p);
0130              sf12=[2.033878107482741, 1.790784258437183];
0131      case 3; fresp = @(f) f<1;
0132              p.FR = FR; p.padding = 0;
0133              sf= freq_filt(s,fresp, p);
0134              sf12=[0.930430975008685, 0.910716555382052];
0135      case 4; fresp = @(f) f<1;
0136              p.FR = FR; p.padding = 2;
0137              sf= freq_filt(s,fresp, p);
0138              sf12=[1.977191625314283, 1.912923819876781];
0139      end
0140      unit_test_cmp(sprintf('ff%02d',fnum),sf(1:2), sf12,1e-13)
0141      plot(t,sf,'LineWidth',2);
0142   end
0143   hold off;
0144 
0145   sv = s.*[2;1;3;4];
0146   fresp = @(f) f<10;
0147   sf= freq_filt(sv,fresp, FR);
0148   sf12=[1.879618111164496, 1.275026363303550];
0149   unit_test_cmp('ff10',sf(2,1:2), sf12,1e-13)
0150 
0151   sf= freq_filt(sv,fresp, FR, 2);
0152   unit_test_cmp('ff11',sf(2,1:2), sf12,1e-13)
0153 
0154   sf= freq_filt(sv',fresp, FR, 1);
0155   unit_test_cmp('ff12',sf(1:2,2), sf12',1e-13)
0156 
0157   so = struct('type','data','meas',sv);
0158   sf= freq_filt(so,fresp, FR, 2);
0159   unit_test_cmp('ff21',sf.meas(2,1:2), sf12,1e-13)
0160   
0161   so = struct('type','image','elem_data',sv);
0162   sf= freq_filt(so,fresp, FR, 2);
0163   unit_test_cmp('ff21',sf.elem_data(2,1:2), sf12,1e-13)
0164 
0165   % TODO add test for complex input