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}]=} resample(@var{x},@var{p},@var{q})
18 ## @deftypefnx {Function File} {@var{y} =} resample(@var{x},@var{p},@var{q},@var{h})
19 ## Change the sample rate of @var{x} by a factor of @var{p}/@var{q}. This is
20 ## performed using a polyphase algorithm. The impulse response @var{h} of the antialiasing
21 ## filter is either specified or either designed with a Kaiser-windowed sinecard.
23 ## Ref [1] J. G. Proakis and D. G. Manolakis,
24 ## Digital Signal Processing: Principles, Algorithms, and Applications,
25 ## 4th ed., Prentice Hall, 2007. Chap. 6
27 ## Ref [2] A. V. Oppenheim, R. W. Schafer and J. R. Buck,
28 ## Discrete-time signal processing, Signal processing series,
29 ## Prentice-Hall, 1999
32 function [y, h] = resample( x, p, q, h )
34 if nargchk(3,4,nargin)
36 elseif any([p q]<=0) || any([p q]~=floor([p q])),
37 error("resample.m: p and q must be positive integers");
40 ## simplify decimation and interpolation factors
42 great_common_divisor=gcd(p,q);
43 if (great_common_divisor>1)
44 p=p/great_common_divisor;
45 q=q/great_common_divisor;
48 ## filter design if required
52 ## properties of the antialiasing filter
54 log10_rejection = -3.0;
55 stopband_cutoff_f = 1.0/(2.0 * max(p,q));
56 roll_off_width = stopband_cutoff_f / 10.0;
58 ## determine filter length
59 ## use empirical formula from [2] Chap 7, Eq. (7.63) p 476
61 rejection_dB = -20.0*log10_rejection;
62 L = ceil((rejection_dB-8.0) / (28.714 * roll_off_width));
67 ideal_filter=2*p*stopband_cutoff_f*sinc(2*stopband_cutoff_f*t);
69 ## determine parameter of Kaiser window
70 ## use empirical formula from [2] Chap 7, Eq. (7.62) p 474
72 if ((rejection_dB>=21) && (rejection_dB<=50))
73 beta = 0.5842 * (rejection_dB-21.0)^0.4 + 0.07886 * (rejection_dB-21.0);
74 elseif (rejection_dB>50)
75 beta = 0.1102 * (rejection_dB-8.7);
80 ## apodize ideal filter response
82 h=kaiser(2*L+1,beta).*ideal_filter;
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("resample.m: the filter h should be a vector");
103 ## pre and postpad filter response
105 nz_pre = floor(q-mod(L,q));
106 hpad = prepad(h,Lh+nz_pre);
108 offset = floor((L+nz_pre)/q);
110 while ceil( ( (Lx-1)*p + nz_pre + Lh + nz_post )/q ) - offset < Ly
113 hpad = postpad(hpad,Lh + nz_pre + nz_post);
116 xfilt = upfirdn(x,hpad,p,q);
117 y = xfilt(offset+1:offset+Ly,:);
136 %! x=sin(2*pi*f0*t' + phi0);
137 %! [y,h]=resample(x,p,q);
138 %! xx=sin(2*pi*f0/r*tt' + phi0);
139 %! t0=ceil((length(h)-1)/2/q);
141 %! err(n+1)=max(abs(y(idx)-xx(idx)));
145 %! idx_inband=1:ceil((1-rolloff/2)*r*N/2)-1;
146 %! assert(max(err(idx_inband))<rejection);
155 %! reject=zeros(N/2,1);
159 %! x=sin(2*pi*f0*t' + phi0);
160 %! [y,h]=resample(x,p,q);
161 %! xx=sin(2*pi*f0/r*tt' + phi0);
162 %! t0=ceil((length(h)-1)/2/q);
164 %! reject(n+1)=max(abs(y(idx)));
168 %! idx_stopband=ceil((1+rolloff/2)*r*N/2)+1:N/2;
169 %! assert(max(reject(idx_stopband))<=rejection);
182 %! x=sin(2*pi*f0*t' + phi0);
183 %! [y,h]=resample(x,p,q);
184 %! xx=sin(2*pi*f0/r*tt' + phi0);
185 %! t0=ceil((length(h)-1)/2/q);
187 %! err(n+1)=max(abs(y(idx)-xx(idx)));
191 %! idx_inband=1:ceil((1-rolloff/2)*r*N/2)-1;
192 %! assert(max(err(idx_inband))<rejection);