get_p_q.m

我认为很不错的语音处理的ｍａｔｌａｂ源代码

function [Fsn, p, q, errors]=get_p_q(Fs, Fsd, tol)
% % get_p_q: Calculate p and q with p <= 3 for fast resampling
% % 
% % Syntax:
% % 
% % [Fsn, p, q, errors]=get_p_q(Fs, Fsd, tol);
% % 
% % *********************************************************************
% %
% % Description
% % 
% % [Fsn, p, q, errors]=get_p_q(Fs, Fsd, tol);
% % Returns the resampling multiples p and q and the resamping rate Fsn (Hz).  
% % The inputs are Fs (Hz) the sampling rate, Fsd (Hz) the desired 
% % sampling rate, and the sampling rate tolerance tol (Hz).  
% % 
% % Teh purpose of the program is to optimize p and q to reduce memory
% % usage and maximize processing speed.  The resample program is very 
% % useful for resampling, however, optimization is required. 
% % 
% % The resample program upsamples by the multiple p then downsamples by 
% % the mutliple q.  The integers p and q must be optimized.
% % 
% % The idea is to select p and q such p is minimized and q can be any
% % integer.  p is limited to being less than 4 and q can be any 
% % counting number.  A tolerance requirement is imposed and the first
% % selection goes to teh padn q combination with minimum p.  
% % 
% % If the tolerance requirement cannot be satisfied then the combination 
% % of p and q that is the closest to the desired sampling rate is 
% % selected.  The errors variable returns 1 when the tolerance 
% % requirement is not satisfied. 
% % 
% % 
% % *********************************************************************
% % 
% % Input Variables
% % 
% % Fs is the sampling rate (Hz);
% % 
% % Fsd is the desired sampling rate (Hz);
% % 
% % tol is the tolerance for the new sampling rate.  
% %        If tolerance requirement is not achievable then errors=1;
% % 
% % *********************************************************************
% % 
% % Output Variables	
% % 
% % Fsn the resampling rate.  See the formula Fsn=Fs*p/q.
% % 
% % p is the integer multiple for upsampling.
% % 
% % q is the integer multiple for downsampling.
% % 
% % errors indicated whether the range tolerance could be satisfied.
% %           0 if satisfied, 1 if not satisfied.  
% % 
% % *********************************************************************
% 
% Example='1';
% 
% Fs=80000;
% Fsd=50000;
% tol=5000;
% 
% [Fsn, p, q, errors]=get_p_q(Fs, Fsd, tol);
% 
% 
% % *********************************************************************
% %
% % Program Written by Edward L. Zechmann 
% % 
% %     date 10 December    2008    
% % 
% % modified 11 December    2008   Updated Comments
% % 
% % *********************************************************************
% % 
% % Please feel free to modify this code.
% % 
% % See also: rat, rats, resample
% % 


if nargin < 1 || isempty(Fs) || ~isnumeric(Fs)
    Fs=50000;
end

if nargin < 2 || isempty(Fsd) || ~isnumeric(Fsd)
    Fsd=5000;
end

if nargin < 3 || isempty(tol) || ~isnumeric(tol)
    tol=5000;
end



% Get the number of sampling rates.
num_sr=length(Fs);

% Initialize the output arrays.
p=zeros(num_sr, 1);
q=zeros(num_sr, 1);
Fsn=zeros(num_sr, 1);


for e1=1:num_sr;

    if Fs(e1) > Fsd-tol && logical(Fs(e1) < Fsd+tol)

        p1=1;
        q1=1;
        errors=0;

    else

        qa=ones(3, 9);
        Fsa=ones(3, 9);

        for e2=1:3;

            low1=Fs(e1)*e2/(Fsd-tol);
            mid1=Fs(e1)*e2/(Fsd);
            high1=Fs(e1)*e2/(Fsd+tol);

            ratios=[low1 mid1 high1];

            qa(e2, :)=[ceil(ratios) round(ratios) floor(ratios)];
            Fsa(e2, :)=Fs(e1).*e2./qa(e2, :);

        end


        [errora ix]=min(abs(Fsa-Fsd), [], 2);


        [iy]=find(errora <= tol);

        if ~isempty(iy)
            p1=iy(1);
            q1=qa(iy(1), ix(iy(1)));
            errors=0;
        else
            [buf iy]=min(errora);
            p1=iy(1);
            q1=qa(iy(1), ix(iy(1)));
            % If the tolerance requirement cannot be satisfied then the errors
            % variable returns 1.  
            errors=1;
        end

    end

    [p1, q1]=rat(p1/q1);

    % Set the final values
    p(e1)=p1;
    q(e1)=q1;

    Fsn(e1)=Fs(e1)*p1/q1;


end