octave_packages/signal-1.1.3/resample.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}]=} 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.
  22 ##
  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
  26 ##
  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
  30 ## @end deftypefn
  31
  32 function  [y, h] = resample( x, p, q, h )
  33
  34   if nargchk(3,4,nargin)
  35     print_usage;
  36   elseif any([p q]<=0) || any([p q]~=floor([p q])),
  37     error("resample.m: p and q must be positive integers");
  38   endif
  39
  40   ## simplify decimation and interpolation factors
  41
  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;
  46   endif
  47
  48   ## filter design if required
  49
  50   if (nargin < 4)
  51
  52     ## properties of the antialiasing filter
  53
  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;
  57
  58     ## determine filter length
  59     ## use empirical formula from [2] Chap 7, Eq. (7.63) p 476
  60
  61     rejection_dB = -20.0*log10_rejection;
  62     L = ceil((rejection_dB-8.0) / (28.714 * roll_off_width));
  63
  64     ## ideal sinc filter
  65
  66     t=(-L:L)';
  67     ideal_filter=2*p*stopband_cutoff_f*sinc(2*stopband_cutoff_f*t);
  68
  69     ## determine parameter of Kaiser window
  70     ## use empirical formula from [2] Chap 7, Eq. (7.62) p 474
  71
  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);
  76     else
  77       beta = 0.0;
  78     endif
  79
  80     ## apodize ideal filter response
  81
  82     h=kaiser(2*L+1,beta).*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("resample.m: the filter h should be a vector");
  96   endif
  97
  98   Lx = rows(x);
  99   Lh = length(h);
 100   L = ( Lh - 1 )/2.0;
 101   Ly = ceil(Lx*p/q);
 102
 103   ## pre and postpad filter response
 104
 105   nz_pre = floor(q-mod(L,q));
 106   hpad = prepad(h,Lh+nz_pre);
 107
 108   offset = floor((L+nz_pre)/q);
 109   nz_post = 0;
 110   while ceil( ( (Lx-1)*p + nz_pre + Lh + nz_post )/q ) - offset < Ly
 111       nz_post++;
 112   endwhile
 113   hpad = postpad(hpad,Lh + nz_pre + nz_post);
 114
 115   ## filtering
 116   xfilt = upfirdn(x,hpad,p,q);
 117   y = xfilt(offset+1:offset+Ly,:);
 118
 119   if isrowvector,
 120      y=y.';
 121   endif
 122
 123 endfunction
 124
 125 %!test
 126 %! N=512;
 127 %! p=3; q=5;
 128 %! r=p/q;
 129 %! NN=ceil(r*N);
 130 %! t=0:N-1;
 131 %! tt=0:NN-1;
 132 %! err=zeros(N/2,1);
 133 %! for n = 0:N/2-1,
 134 %!   phi0=2*pi*rand;
 135 %!   f0=n/N;
 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);
 140 %!   idx=t0+1:NN-t0;
 141 %!   err(n+1)=max(abs(y(idx)-xx(idx)));
 142 %! endfor;
 143 %! rolloff=.1;
 144 %! rejection=10^-3;
 145 %! idx_inband=1:ceil((1-rolloff/2)*r*N/2)-1;
 146 %! assert(max(err(idx_inband))<rejection);
 147
 148 %!test
 149 %! N=512;
 150 %! p=3; q=5;
 151 %! r=p/q;
 152 %! NN=ceil(r*N);
 153 %! t=0:N-1;
 154 %! tt=0:NN-1;
 155 %! reject=zeros(N/2,1);
 156 %! for n = 0:N/2-1,
 157 %!   phi0=2*pi*rand;
 158 %!   f0=n/N;
 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);
 163 %!   idx=t0+1:NN-t0;
 164 %!   reject(n+1)=max(abs(y(idx)));
 165 %! endfor;
 166 %! rolloff=.1;
 167 %! rejection=10^-3;
 168 %! idx_stopband=ceil((1+rolloff/2)*r*N/2)+1:N/2;
 169 %! assert(max(reject(idx_stopband))<=rejection);
 170
 171 %!test
 172 %! N=1024;
 173 %! p=2; q=7;
 174 %! r=p/q;
 175 %! NN=ceil(r*N);
 176 %! t=0:N-1;
 177 %! tt=0:NN-1;
 178 %! err=zeros(N/2,1);
 179 %! for n = 0:N/2-1,
 180 %!   phi0=2*pi*rand;
 181 %!   f0=n/N;
 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);
 186 %!   idx=t0+1:NN-t0;
 187 %!   err(n+1)=max(abs(y(idx)-xx(idx)));
 188 %! endfor;
 189 %! rolloff=.1;
 190 %! rejection=10^-3;
 191 %! idx_inband=1:ceil((1-rolloff/2)*r*N/2)-1;
 192 %! assert(max(err(idx_inband))<rejection);