Description of freq

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 %   e.g. fresp = @(f) (1+(f/10).^6).^(-0.5);  % Butter 2*3rd order
0010 % FR  = frame rate (or use FR parameter on data structures)
0011 %    or FR is a structure with fields
0012 %     .FR = Frame Rate
0013 %     .padding = length of padding (in seconds)
0014 % dim = dimension along which to filter (default is 2)
0015 %     .padding_type = zero|extend|linear
0016 %       pad with zeros, extend the last values (default), or linearly fit between start/end
0017 
0018 % (C) Andy Adler 2019. License: GPL v2 or v3.
0019 % $Id: freq_filt.m 6977 2024-11-03 15:23:53Z aadler $
0020 
0021 % if FR not given, see if it's a parameter
0022 
0023 if ischar(vin) && strcmp(vin,'UNIT_TEST'); do_unit_test; return; end
0024 
0025 p.padding = 30; % seconds default;
0026 p.padding_type = 'extend';
0027 if nargin<3; 
0028    p.FR = vin.FR;
0029 else
0030    if isnumeric(FR)
0031      p.FR = FR;
0032    elseif isstruct(FR);
0033      for ff= fieldnames(FR)'; % wish matlab could do this easily
0034        p.(ff{1}) = FR.(ff{1});
0035      end
0036    else
0037      error('Don''t understand parameter FR');
0038    end
0039      
0040       
0041 end
0042 if nargin>=4
0043    p.dim = dim;
0044 elseif ~isfield(p,'dim');
0045    p.dim = 2;
0046 end
0047 
0048 if isstruct(vin) && isfield(vin,'type');
0049    switch vin.type
0050      case 'data'
0051         vin.meas = do_freq_filt(vin.meas, fresp, p);
0052         vfilt = vin;
0053      case 'image'
0054         vin.elem_data = do_freq_filt(vin.elem_data, fresp, p);
0055         vfilt = vin;
0056      otherwise
0057         error('Can''t process object of type %',vin.type);
0058    end
0059 elseif isnumeric(vin)
0060    vfilt = do_freq_filt(vin, fresp, p);
0061 else 
0062    error('can''t process object');
0063 end
0064       
0065 
0066 
0067 % Filter in the frequency direction
0068 function s = do_freq_filt(s,fresp, p)
0069   rs = isreal(s);
0070   s = padsignal(s,p);
0071   f = fft(s,[],p.dim);
0072   fax = freq_axis_filt(p.FR,size(s,p.dim),fresp);
0073   f = f .* fshape(p,fax);
0074   s= ifft(f,[],p.dim);
0075   if rs % is signal is real
0076      if norm(imag(s(:))) / norm(real(s(:))) > 1e-12
0077         error('FFT filter has imag output');
0078      end
0079      s = real(s);
0080   end
0081   s = unpadsignal(s,p);
0082 
0083 function s = fshape(p,s);
0084   fsh = ones(1,length(size(s)));
0085   fsh(p.dim) = prod(size(s));
0086   s= reshape(s,fsh);
0087 
0088 function plen = pad_len(p);
0089   plen = round(p.FR * p.padding);
0090 
0091 % Add padding connecting the last to the first sample
0092 function s = padsignal(s,p);
0093   lsup = linspace(0,1,pad_len(p)); lsup=fshape(p,lsup);
0094   lsdn = linspace(1,0,pad_len(p)); lsdn=fshape(p,lsdn);
0095   switch p.padding_type
0096     case 'zero';
0097       pad= calc_padding(s,0*lsdn,0*lsup,p);
0098     case 'extend';
0099       pad= calc_padding(s,lsdn,lsup,p);
0100     case 'linear';
0101       mn  = mean(s,p.dim); ss = s-mn;
0102       xx  = linspace(-1,1,size(ss,p.dim));
0103       perm = ones(1,5); perm(p.dim) = length(xx); % do we ever want d>5??
0104       xx  = reshape(xx,perm);
0105       a   = mean(ss.*xx,p.dim) ./mean(xx.^2,p.dim);
0106       pad = (mn-a).*lsdn + (mn+a).*lsup;
0107     otherwise; 
0108       error('padding type not understood');
0109   end
0110   
0111   s = cat(p.dim,s,pad);
0112 
0113 function pad= calc_padding(s,lsdn,lsup,p)
0114   switch p.dim % wish I knew how to generalize in m*lab
0115     case 1; pad = s(1,:,:).*lsdn + s(end,:,:).*lsup;
0116     case 2; pad = s(:,1,:).*lsdn + s(:,end,:).*lsup;
0117     case 3; pad = s(:,:,1).*lsdn + s(:,:,end).*lsup;
0118     otherwise; error('Problem with dim above 3');
0119   end
0120 
0121 % Add padding connecting the last to the first sample
0122 function s = unpadsignal(s,p);
0123   plen = round(p.FR * p.padding);
0124   switch p.dim % wish I knew how to generalize in m*lab
0125     case 1; s(end+1-(1:plen),:,:) = [];
0126     case 2; s(:,end+1-(1:plen),:) = [];
0127     case 3; s(:,:,end+1-(1:plen)) = [];
0128     otherwise; error('Problem with dim above 3');
0129   end
0130   
0131 function fax = freq_axis_filt( FR, lD, fresp);
0132   fax = linspace(0,FR,lD+1);
0133   fax(end)=[];
0134   fax(fax>FR/2) = fax(fax>FR/2) - FR;
0135   fax = feval(fresp, abs(fax));
0136 
0137 function do_unit_test
0138   clf; subplot(3,1,1);
0139   FR = 100;
0140   t = (0:1.1e3)/FR;
0141   s = sin(2*pi*5*t) + 2*cos(2*pi*15*t) +  3*sin(2*pi*0.1*t) + 1;
0142   subplot(211); plot(t,s); hold on;
0143   subplot(212); plot(t,s); hold on; xlim(max(t)-[0.5,0]);
0144   pp.FR = FR;
0145   pp.padding_type = 'extend';
0146   pp.padding = 1;
0147 
0148   fnum=0; while true; fnum=fnum+1; switch fnum
0149      case 1; fresp = @(f) f<10;
0150              sf= freq_filt(s,fresp, pp);
0151              sf12=[1.879618111164496, 1.275026363303550];
0152 
0153      case 2; fresp = @(f) (f<1) + (f>=1)./(f+eps);
0154              p.padding_type = 'extend';
0155              p.FR = FR; p.padding = 1;
0156              sf= freq_filt(s,fresp, p);
0157              sf12=[2.033878107482741, 1.790784258437183];
0158      case 3; fresp = @(f) f<1;
0159              p.padding_type = 'extend';
0160              p.FR = FR; p.padding = 0;
0161              sf= freq_filt(s,fresp, p);
0162              sf12=[0.930430975008685, 0.910716555382052];
0163      case 4; fresp = @(f) f<1;
0164              p.padding_type = 'extend';
0165              p.FR = FR; p.padding = 2;
0166              sf= freq_filt(s,fresp, p);
0167              sf12=[1.977191625314283, 1.912923819876781];
0168      case 5; fresp = @(f) f<1;
0169              p.padding_type = 'linear';
0170 % TODO- debug for padding size differences
0171              p.FR = FR; p.padding = 3;
0172              sf= freq_filt(s,fresp, p);
0173              sf12=[1.231446973826249, 1.277407267967623]-2;
0174      case 6; fresp = @(f) f<1;
0175              p.padding_type = 'zero';
0176 % TODO- debug for padding size differences
0177              p.FR = FR; p.padding = 1;
0178              sf= freq_filt(s,fresp, p);
0179              sf12=[1.818358872667770, 1.846165834456450]-2;
0180      case 7; fresp = @(f) f<10;
0181              sf= freq_filt(s,fresp, FR);
0182              sf12=[1.884878325803839, 1.271527154193287];
0183      otherwise; break
0184      end
0185      unit_test_cmp(sprintf('ff%02d',fnum),sf(1:2)-1, sf12,1e-13)
0186      subplot(211); plot(t,sf);
0187      subplot(212); plot(t,sf);
0188   end
0189   legend('0','1','2','3','4','5','6');
0190   hold off;
0191 
0192   sv = s.*[2;1;3;4];
0193   fresp = @(f) f<10;
0194   sf= freq_filt(sv,fresp, pp);
0195   sf12=[1.879618111164496, 1.275026363303550]+1;
0196   unit_test_cmp('ff10',sf(2,1:2), sf12,1e-13)
0197 
0198   sf= freq_filt(sv,fresp, pp, 2);
0199   unit_test_cmp('ff11',sf(2,1:2), sf12,1e-13)
0200 
0201   sf= freq_filt(sv',fresp, pp, 1);
0202   unit_test_cmp('ff12',sf(1:2,2), sf12',1e-13)
0203 
0204   so = struct('type','data','meas',sv);
0205   sf= freq_filt(so,fresp, pp, 2);
0206   unit_test_cmp('ff21',sf.meas(2,1:2), sf12,1e-13)
0207   
0208   so = struct('type','image','elem_data',sv);
0209   sf= freq_filt(so,fresp, pp, 2);
0210   unit_test_cmp('ff21',sf.elem_data(2,1:2), sf12,1e-13)
0211 
0212   % TODO add test for complex input
freq_filt

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
