% EE3221 W7D3 1/29/2021 Dr. Durant

% z_freq_scaling

% There was a question today about whether z-domain transfer functions can be
% scaled in frequency. They can, with the same caveats about complex signals
% that we had in EE3032 in continuous time.

% From Table 7-6.5, the relevant property is
% a^n x[n] u[n] <--> X(z/a)

% That is, if you multiply x[n] by a geometric series, you are also affecting
% its z-transform by replacing all the z's with z/a.

% Let a = e^(jw0), where w0 is the shift frequency in radians per sample.

% This accomplishes a shift of the frequency response by w0, and rotates all
% the poles and zeros by w0. But, it also takes a real impulse response and
% makes it complex and, correspondingly, breaks the conjugate symmetry of the
% transfer function, which only exists for real impulse responses:
% H(e^(jw) =? H*(e^-(jw)).

% The same real-to-complex issues occur in continuous time. See Table 5-7, Property 5 (Frequency shift).

% There are many approaches to filter design and transformation. To shift a
% frequency response (roughly) up by w0 radians/sample, you could rotate all the
% poles at non-negative frequencies (0<=w<=pi) by +w0 and all the conjugate poles
% at non-positive frequencies by -w0. This has the complication that you won't get
% the intended effect if you rotate anything past w=0 or w=pi=Nyquist. I write
% "roughly" since the part of the response due to the distance to the conjugate poles
% and zeros changes in a non-trivial way since H(z) is being evaluated on a circle.

% For example, we looked at the following (stable) filter today:

clearvars
close all

a = [0.9157    0.7922    0.9595    0.6557]; % a0..a3
b = [0.9649    0.1576    0.9706    0.9572]; % b0..b3

po = roots(a)
ze = roots(b)

% Work in polar form...
po_mag = abs(po) % observe they're < 1, so stable H(z), h(n)
po_angle = angle(po)/pi % fraction between 0 and 1 (Nyquist) - corresponds to frequency in rad/sample

ze_mag = abs(ze)
ze_angle = angle(ze)/pi

% Shift frequency up, maintaining conjugate symmetry
w_shift = 0.2;
po_angle = [po_angle(1)+w_shift; po_angle(2)-w_shift; po_angle(3)]
ze_angle = [ze_angle(1)+w_shift; ze_angle(2)-w_shift; ze_angle(3)]

% Put mag and angle together into a complex number
po_new = po_mag .* exp(1j*pi*po_angle)
ze_new = ze_mag .* exp(1j*pi*ze_angle)

% poly takes roots and gives a monic polynomial; the leading coefficient's
% scaling was lost when we called roots, so we restore it here...
a_new = poly(po_new) * a(1);
b_new = poly(ze_new) * b(1);

% Floating point roundoff on the rotated poles and zeros will cause tiny
% imaginary portions in new a and b coefficients. Confirm they're
% negligibly small and then discard them.
assert(norm(imag(a_new))/norm(real(a_new)) < 10*eps, 'Unexpected imaginary component in new a')
assert(norm(imag(b_new))/norm(real(b_new)) < 10*eps, 'Unexpected imaginary component in new b')
a_new = real(a_new) % Keep just the real part
b_new = real(b_new)

if false
    % Simple approach with separate figures since freqz doesn't work with subplot.
    figure, zplane(b,a), title('Original')
    figure, zplane(b_new,a_new), title('Shifted')
    figure, freqz(b,a), title('Original')
    figure, freqz(b_new,a_new), title('Shifted')
else
    % More complicated approach that uses subplots. freqz is used just to return
    % values (not plot) in this case and then user code handles the plotting in subplots.
    [H,w] = freqz(b,a);
    [H_new,w_new] = freqz(b_new,a_new);
    figure
    subplot(221),zplane(b,a),title('Original H(z) Roots')
    subplot(222),plot(w/pi,mag2db(abs(H)),w_new/pi,mag2db(abs(H_new))),grid
    legend('Original', 'Shifted')
    ylabel('Magnitude (dB)')
    subplot(223),zplane(b_new,a_new),title('Shifted H(z) Roots')
    subplot(224),plot(w/pi,rad2deg(angle(H)),w_new/pi,rad2deg(angle(H_new))),grid
    xlabel('Normalized Frequency (\times\pi rad/sample)')
    ylabel('Phase (degrees')
end