f_rateconv.m

digital signal processing常用工具箱

function [y,b] = f_rateconv (x,fs,L,M,m,f_type,alpha)

%F_RATECONV: Convert sampling rate by rational factor L/M.
%
% Usage: [y,b] = f_rateconv (x,fs,L,M,m,f_type,alpha)
%
% Inputs: 
%         x      = a vector of length P containing the 
%                  input samples
%         fs     = sampling frequency of x
%         L      = an integer specifying the numerator 
%                  of the conversion factor (L >= 1).
%         M      = an integer specifying the denominator
%                  of the conversion factor (M >= 1)
%         m      = the order of the lowpass FIR anti-
%                  aliasing anti-imaging filter. 
%         f_type = the FIR filer type to be used:
%
%                  0 = windowed (rectangular)
%                  1 = windowed (Hanning)
%                  2 = windowed (Hamming)
%                  3 = windowed (Blackman)
%                  4 = frequency-sampled
%                  5 = least-squares
%                  6 = equiripple
%
%          alpha = an optional scaling factor for the cutoff
%                  frequency of the FIR filter.  Default:
%                  alpha = 1.  If present, the cutoff 
%                  frequency used for the anti-aliasing 
%                  anti-imaging filter H_0(z) is
% 
%                  F_c = alpha*min(fs/(2L),fs/(2M))
% Outputs: 
%          y   = a 1 by N vector containing the output 
%                samples. Here N = floor(LP/M).
%          b   = a 1 by (m+1) vector containing FIR 
%                filter coefficients 
%
% Notes: If L or M are relatively large (e.g. 10 or 
%        higher), then it is the responsibilty of the user
%        to perform the rate conversion in stages using 
%        multiple calls to f_rateconv.  Otherwise, the 
%        required value for m can be very large.
%
% See also: F_DECIMATE, F_INTERPOL

% Initialize

L = f_clip (L,1,L);
M = f_clip (M,1,M);
m = f_clip (m,2,m);
m = 2*floor(m/2);
f_type = f_clip (f_type,0,6);
P = length (x);
T = 1/fs;
if M == L
    L = 1;
    M = 1;
end

% Design lowpass FIR filter

if nargin < 7
    alpha = 1;
end
F_0 = alpha*min(fs/(2*L),fs/(2*M));
p = [0,F_0,F_0,0];
sym = 0;
switch f_type
    
    case {0,1,2,3},
        sym = 0;
        win = f_type;
        b = f_firideal (0,F_0,m,fs,win);        
    case 4,
        q = floor(m/2)+1;
        F = linspace (0,fs/2,q);
        A = f_firamp (F,fs,p);
        b = f_firsamp (A,m,fs,sym);
    case 5,
        q = 2*m;
        F = linspace (0,fs/2,q);
        A = f_firamp (F,fs,p);
        b = f_firls (F,A,m,fs);
    case 6,
        F_p = F_0;
        F_s = F_0+fs/m;
        delta = 0.05;
        [b,m0] = f_firparks (m,F_p,F_s,delta/2,delta,0,fs);
end
b = L*b;        

% Perform sampling rate conversion

N = floor(L*P/M);
y = zeros(N,1);
hstr = sprintf ('Performing rate conversion by %.4f',L/M);
h = waitbar (0,hstr);
for k = 0 : N-1
   waitbar (k/(N-1),h)
   for i = 0 : min(m,M*k)
      k0 = M*k-i;
      if mod(k0,L) == 0
         i0 = floor(k0/L);
         if i0 < P;
            y(k+1) = y(k+1) + b(i+1)*x(i0+1);
         end
      end
   end
end
close (h);