1 ## Copyright (C) 2008 Eric Chassande-Mottin, CNRS (France) <ecm@apc.univ-paris7.fr>
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
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
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/>.
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.
22 ## @seealso{circshift}
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
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
32 function [y, h] = fracshift( x, d, h )
34 if nargchk(2,3,nargin)
38 ## if the delay is an exact integer, use circshift
44 ## filter design if required
48 ## properties of the interpolation filter
50 log10_rejection = -3.0;
51 stopband_cutoff_f = 1.0 / 2.0;
52 roll_off_width = stopband_cutoff_f / 10;
54 ## determine filter length
55 ## use empirical formula from [1] Chap 7, Eq. (7.63) p 476
57 rejection_dB = -20.0*log10_rejection;
58 L = ceil((rejection_dB-8.0) / (28.714 * roll_off_width));
63 ideal_filter=2*stopband_cutoff_f*sinc(2*stopband_cutoff_f*(t-(d-fix(d))));
65 ## determine parameter of Kaiser window
66 ## use empirical formula from [1] Chap 7, Eq. (7.62) p 474
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);
76 ## apodize ideal (sincard) filter response
79 t = (0 : m)' - (d-fix(d));
80 t = 2 * beta / m * sqrt (t .* (m - t));
81 w = besseli (0, t) / besseli (0, beta);
86 ## check if input is a row vector
88 if ((rows(x)==1) && (columns(x)>1))
93 ## check if filter is a vector
95 error("fracshift.m: the filter h should be a vector");
103 ## pre and postpad filter response
106 hpad = postpad(hpad,Ly + offset);
109 xfilt = upfirdn(x,hpad,1,1);
110 y = xfilt(offset+1:offset+Ly,:);
112 y=circshift(y,fix(d));
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));
137 %! idx_inband=1:ceil((1-rolloff)*N/2)-1;
138 %! assert(max(err(idx_inband))<rejection);
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));
157 %! idx_inband=1:ceil((1-rolloff)*N/2)-1;
158 %! assert(max(err(idx_inband))<rejection);
168 %! n(N/2+(-t:t))=randn(2*t+1,1);
169 %! [b a]=butter(10,.25);
171 %! n1=fracshift(n,d1);
172 %! n1=resample(n1,p,q);
173 %! n2=resample(n,p,q);
174 %! n2=fracshift(n2,d2);
177 %! assert(max(err)<rejection);