octave_packages/signal-1.1.3/fracshift.m

   1 ## Copyright (C) 2008 Eric Chassande-Mottin, CNRS (France) <ecm@apc.univ-paris7.fr>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File}  {[@var{y} @var{h}]=} fracshift(@var{x},@var{d})
  18 ## @deftypefnx {Function File} {@var{y} =} fracshift(@var{x},@var{d},@var{h})
  19 ## Shift the series @var{x} by a (possibly fractional) number of samples @var{d}.
  20 ## The interpolator @var{h} is either specified or either designed with a Kaiser-windowed sinecard.
  21 ## @end deftypefn
  22 ## @seealso{circshift}
  23
  24 ## Ref [1] A. V. Oppenheim, R. W. Schafer and J. R. Buck,
  25 ## Discrete-time signal processing, Signal processing series,
  26 ## Prentice-Hall, 1999
  27 ##
  28 ## Ref [2] T.I. Laakso, V. Valimaki, M. Karjalainen and U.K. Laine
  29 ## Splitting the unit delay, IEEE Signal Processing Magazine,
  30 ## vol. 13, no. 1, pp 30--59 Jan 1996
  31
  32 function  [y, h] = fracshift( x, d, h )
  33
  34   if nargchk(2,3,nargin)
  35     print_usage;
  36   endif;
  37
  38   ## if the delay is an exact integer, use circshift
  39   if d==fix(d)
  40     y=circshift(x,d);
  41     return
  42   endif;
  43
  44   ## filter design if required
  45
  46   if (nargin < 4)
  47
  48     ## properties of the interpolation filter
  49
  50     log10_rejection = -3.0;
  51     stopband_cutoff_f = 1.0 / 2.0;
  52     roll_off_width = stopband_cutoff_f / 10;
  53
  54     ## determine filter length
  55     ## use empirical formula from [1] Chap 7, Eq. (7.63) p 476
  56
  57     rejection_dB = -20.0*log10_rejection;
  58     L = ceil((rejection_dB-8.0) / (28.714 * roll_off_width));
  59
  60     ## ideal sinc filter
  61
  62     t=(-L:L)';
  63     ideal_filter=2*stopband_cutoff_f*sinc(2*stopband_cutoff_f*(t-(d-fix(d))));
  64
  65     ## determine parameter of Kaiser window
  66     ## use empirical formula from [1] Chap 7, Eq. (7.62) p 474
  67
  68     if ((rejection_dB>=21) && (rejection_dB<=50))
  69       beta = 0.5842 * (rejection_dB-21.0)^0.4 + 0.07886 * (rejection_dB-21.0);
  70     elseif (rejection_dB>50)
  71       beta = 0.1102 * (rejection_dB-8.7);
  72     else
  73       beta = 0.0;
  74     endif
  75
  76     ## apodize ideal (sincard) filter response
  77
  78     m = 2*L;
  79     t = (0 : m)' - (d-fix(d));
  80     t = 2 * beta / m * sqrt (t .* (m - t));
  81     w = besseli (0, t) / besseli (0, beta);
  82     h = w.*ideal_filter;
  83
  84   endif
  85
  86   ## check if input is a row vector
  87   isrowvector=false;
  88   if ((rows(x)==1) && (columns(x)>1))
  89      x=x(:);
  90      isrowvector=true;
  91   endif
  92
  93   ## check if filter is a vector
  94   if ~isvector(h)
  95     error("fracshift.m: the filter h should be a vector");
  96   endif
  97
  98   Lx = length(x);
  99   Lh = length(h);
 100   L = ( Lh - 1 )/2.0;
 101   Ly = Lx;
 102
 103   ## pre and postpad filter response
 104   hpad = prepad(h,Lh);
 105   offset = floor(L);
 106   hpad = postpad(hpad,Ly + offset);
 107
 108   ## filtering
 109   xfilt = upfirdn(x,hpad,1,1);
 110   y = xfilt(offset+1:offset+Ly,:);
 111
 112   y=circshift(y,fix(d));
 113
 114   if isrowvector,
 115      y=y.';
 116   endif
 117
 118 endfunction
 119
 120 %!test
 121 %! N=1024;
 122 %! d=1.5;
 123 %! t=(0:N-1)-N/2;
 124 %! tt=t-d;
 125 %! err=zeros(N/2,1);
 126 %! for n = 0:N/2-1,
 127 %!   phi0=2*pi*rand;
 128 %!   f0=n/N;
 129 %!   sigma=N/4;
 130 %!   x=exp(-t'.^2/(2*sigma)).*sin(2*pi*f0*t' + phi0);
 131 %!   [y,h]=fracshift(x,d);
 132 %!   xx=exp(-tt'.^2/(2*sigma)).*sin(2*pi*f0*tt' + phi0);
 133 %!   err(n+1)=max(abs(y-xx));
 134 %! endfor;
 135 %! rolloff=.1;
 136 %! rejection=10^-3;
 137 %! idx_inband=1:ceil((1-rolloff)*N/2)-1;
 138 %! assert(max(err(idx_inband))<rejection);
 139
 140 %!test
 141 %! N=1024;
 142 %! d=7/6;
 143 %! t=(0:N-1)-N/2;
 144 %! tt=t-d;
 145 %! err=zeros(N/2,1);
 146 %! for n = 0:N/2-1,
 147 %!   phi0=2*pi*rand;
 148 %!   f0=n/N;
 149 %!   sigma=N/4;
 150 %!   x=exp(-t'.^2/(2*sigma)).*sin(2*pi*f0*t' + phi0);
 151 %!   [y,h]=fracshift(x,d);
 152 %!   xx=exp(-tt'.^2/(2*sigma)).*sin(2*pi*f0*tt' + phi0);
 153 %!   err(n+1)=max(abs(y-xx));
 154 %! endfor;
 155 %! rolloff=.1;
 156 %! rejection=10^-3;
 157 %! idx_inband=1:ceil((1-rolloff)*N/2)-1;
 158 %! assert(max(err(idx_inband))<rejection);
 159
 160 %!test
 161 %! N=1024;
 162 %! p=6;
 163 %! q=7;
 164 %! d1=64;
 165 %! d2=d1*p/q;
 166 %! t=128;
 167 %! n=zeros(N,1);
 168 %! n(N/2+(-t:t))=randn(2*t+1,1);
 169 %! [b a]=butter(10,.25);
 170 %! n=filter(b,a,n);
 171 %! n1=fracshift(n,d1);
 172 %! n1=resample(n1,p,q);
 173 %! n2=resample(n,p,q);
 174 %! n2=fracshift(n2,d2);
 175 %! err=abs(n2-n1);
 176 %! rejection=10^-3;
 177 %! assert(max(err)<rejection);