er_spectrans.m

dsp工具箱

function g = er_spectrans(f, origf, newfs, s, m)
% function g = er_spectrans(f, origf, newfs, s, m)
%   Function to modify digital LPF coefficient vectors given in struct F 
%   using a spectral transformation described by original frequency pts. in
%   ORIGFPT (0.25 for prototype LPF), desired frequency pts. in NEWFPTS, 
%   shift parameter M, and scalar S (S = +/- 1).  Returns transfer function 
%   coefficient vectors in G.
%
%   Author:  Evan Ruzanski, CU-Boulder, ECEN5632 MATLAB assignment, FA2004

% Unpack input struct
num1 = f(1).tf_complete;
den1 = f(2).tf_complete;

% Unnormalize frequency points
origfpt = origf*pi;
newfpts = newfs*pi;

% Find parameter vectors
if (length(origfpt) == 1) & (length(newfpts) == 1) % Lowpass case
    % z^(-2) coefficient is A(2)
    diffa = origfpt - newfpts;
    suma = origfpt + newfpts;
    numa = sin(diffa/2);
    dena = sin(suma/2);
    A(2) = numa/dena;
    A(2) = -A(2);
    A(1) = 1;
elseif (length(origfpt) == 1) & (length(newfpts) >= 2) % Bandpass case
    a = cos((newfpts(2) + newfpts(1))/2) / cos((newfpts(2) - newfpts(1))/2);
    b = cot((newfpts(2) - newfpts(1))/2)*tan(origfpt/2);
    % z^(-2) coefficient is ALPHA(2)
    A(3) = (b - 1)/(b + 1);
    A(2) = -2*a*b/(b + 1);
    A(1) = 1;
else
    error('Check size of frequency vectors!')
    return
end

% Perform spectral transformation
n = length(num1) - 1; % Order of numerator, from given filter, assume order of numerator = order of denominator
for k = 1:n + 1
    A = s*A;
    % Step 1: Handle exponents, first allpass function in equation
    if (k == n + 1)
        t1 = 1;
    elseif (k == n)
        t1 = A;
    else
        t1 = A;
        for q = 1:(n - k)
            tt1 = conv(t1,A);
            t1 = tt1;
        end
    end
    A = A/s;
    % Step 2: Handle exponents, second reversed allpass function in
    % equation
    if (k == 1)
        t2 = 1;
    elseif (k == 2)
        t2 = fliplr(A);
    else
        t2 = fliplr(A);
        for p = 1:(k - 2)
            tt2 = conv(t2,fliplr(A));
            t2 = tt2;
        end
    end
    % Step 3: Combine convolution results
    AAA = conv(t1,t2);
    % Step 4: Shift results and zeropad accordingly
    BBB = [zeros(1,m*((n + 1) - k)) AAA];
    AAA = [zeros(1,m*((n + 1) - k)) AAA];
    % Step 5: Multiply by appropriate coefficient
    F(k,1:length(BBB)) = num1(k)*BBB;
    G(k,1:length(AAA)) = den1(k)*AAA;
end
% Setp 6: Create polynomials
ss = size(G);
for k = 1:ss(2)
    bspec(k) = sum(F(:,k));
    aspec(k) = sum(G(:,k));
end
bspec = bspec/aspec(1);
aspec = aspec/aspec(1);
% Repack output struct
f1(1).tf_complete = bspec;
f1(2).tf_complete = aspec; 
g = f1;