Подробно реши пятое задание из прикрепленного файл...

We solve Problem 5 (CORDIC) from your PDF: compute
a) arctan(1.1)
b) abs(1.2 + 2.1j)
using a 10-step CORDIC algorithm, justify the method, and compare with the exact values. (HW3_2025.pdf)

---

1. Short recap of CORDIC (vectoring mode)

CORDIC computes trig functions and vector operations using only:

• additions/subtractions
• bit shifts (multiplying by powers of 2)
• a small lookup table of angles α_i = atan(2⁻ᶦ)

There are two modes:

• Rotation mode – rotate a vector by a given angle.
• Vectoring mode – rotate a vector until its y-component becomes zero, and accumulate the angle it was rotated by.

For this problem we use vectoring mode, because:

• To get atan(1.1) we start with a vector whose slope is 1.1 and drive y → 0.
• To get |1.2 + 2.1j| we start with vector (x, y) = (1.2, 2.1) and again drive y → 0. The final x (after scaling correction) is the magnitude.

---

2. CORDIC vectoring formulas

Precompute the list of angles:

$$\alpha_i = \arctan(2^{-i}),\quad i = 0, 1, \dots, N-1,$$

where here N = 10.

CORDIC vectoring mode updates:

• Initialize

  $$x_0 = x_\text{in},\quad y_0 = y_\text{in},\quad z_0 = 0$$

• For each iteration $i = 0,1,\dots,N-1$:

  1. Choose direction

     $$d_i = \operatorname{sign}(y_i) = \begin{cases} +1, & y_i \ge 0\\[4pt] -1, & y_i < 0 \end{cases}$$

  2. Update

     $$\begin{aligned} x_{i+1} &= x_i + d_i\, y_i\, 2^{-i} \\ y_{i+1} &= y_i - d_i\, x_i\, 2^{-i} \\ z_{i+1} &= z_i + d_i\, \alpha_i \end{aligned}$$

After enough iterations:

• $y_N \approx 0$
• $z_N \approx \arctan\left(\dfrac{y_0}{x_0}\right)$
• The length is scaled by a constant factor

$$G_N = \prod_{i=0}^{N-1} \sqrt{1 + 2^{-2i}} \approx 1.64676$$

so

$$r = \sqrt{x_0^2 + y_0^2} \approx \frac{x_N}{G_N}$$

Equivalently, define

$$K_N = \frac{1}{G_N} = \prod_{i=0}^{N-1} \frac{1}{\sqrt{1+2^{-2i}}} \approx 0.607253$$

and then

$$r \approx K_N \, x_N$$

We will use this form in code.

---

3. Part (a): CORDIC for $\arctan(1.1)$

We want $\theta = \arctan(1.1)$.

Choose initial vector with slope 1.1:

$$x_0 = 1,\quad y_0 = 1.1,\quad z_0 = 0$$

Run the 10-step vectoring algorithm.

First few iterations (illustration)

(Values rounded to 6 decimal digits.)

Let $\alpha_i = \arctan(2^{-i})$.

i	yi sign	xi+1	yi+1	zi+1 (rad)
0	y₀ > 0	2.100000	0.100000	0.785398
1	y₁ > 0	2.150000	−0.950000	1.249046
2	y₂ < 0	2.387500	−0.412500	1.004067
3	y₃ < 0	2.439063	−0.114062	0.879712
4	y₄ < 0	2.446191	0.038379	0.817293
…	…	…	…	…
9	y₉ < 0	2.448079	−0.004605	0.834862

After 10 iterations:

• $z_{10} \approx 0.834862\ \text{rad}$

This is our CORDIC approximation:

$$\boxed{\arctan(1.1)_{\text{CORDIC, 10 steps}} \approx 0.834862\ \text{rad}}$$

Comparison with the exact value

Exact value (as computed by high-precision math / MATLAB atan):

$$\arctan(1.1) \approx 0.832981\ \text{rad}$$

Error:

$$\Delta\theta \approx 0.834862 - 0.832981 \approx 1.88\times 10^{-3}\ \text{rad}$$

In degrees:

• Exact: ≈ 47.726°
• CORDIC: ≈ 47.834°
• Error: ≈ 0.108°

So 10 iterations already give ~3 decimal digits of accuracy for the angle.

---

4. Part (b): CORDIC for $|1.2 + 2.1j|$

We treat the complex number as a 2D vector:

$$x_0 = \Re = 1.2,\quad y_0 = \Im = 2.1,\quad z_0 = 0$$

Run the same vectoring algorithm. After 10 iterations we obtain some $x_{10}, y_{10}, z_{10}$ with $y_{10} \approx 0$. Because of the accumulated scaling,

$$r_{\text{CORDIC}} \approx K_{10}\, x_{10}$$

with

$$K_{10} = \prod_{i=0}^{9} \frac{1}{\sqrt{1+2^{-2i}}} \approx 0.607253$$

Numerically (10 iterations):

• $r_{\text{CORDIC}} \approx 2.4186767488$

Exact magnitude:

$$|1.2 + 2.1j| = \sqrt{1.2^2 + 2.1^2} \approx 2.4186773245$$

Error:

• Absolute: ≈ −5.76 × 10⁻⁷
• Relative: ≈ −2.38 × 10⁻⁷ (about 0.000024%)

So the CORDIC magnitude after 10 steps is accurate to at least 6 decimal digits.

(As a bonus, the same run also gives the phase
$\arg(1.2 + 2.1j) \approx 1.05234\ \text{rad}$, close to the exact $1.05165\ \text{rad}$.)

---

5. Why this approach is valid (justification)

1. Angle decomposition:
   Any angle $\theta$ can be represented (approximately) as a sum of signed micro-angles:
   $\theta \approx \sum d_i \alpha_i$ with $\alpha_i = \arctan(2^{-i})$.
   With vectoring mode, the algorithm automatically chooses the signs $d_i$ so that the vector converges to the x-axis, and $z_N$ becomes that sum.

2. Shift-add hardware efficiency:
   Because each micro-rotation uses factors $2^{-i}$, we need only shifts and adds (no multipliers!). CORDIC is historically important in hardware (FPGA, ASIC, DSP).

3. Convergence:
   Each step reduces the magnitude of $y$ roughly by a factor of ~2. Hence, after N steps, the angle error is on the order of $2^{-N}$. With N = 10, error is roughly $\lesssim 10^{-3}\ \text{rad}$, consistent with our numerical results.

4. Fixed gain:
   The overall gain $G_N$ depends only on N, not on the input. We can precompute the correction factor $K_N = 1/G_N$ once and reuse it for all inputs to recover the proper magnitude.

---

6. MATLAB code (10-step CORDIC, arctan & magnitude)

Below is self-contained MATLAB code that:

• computes arctan(1.1) with 10 CORDIC iterations
• computes abs(1.2 + 2.1j) using CORDIC
• compares with MATLAB’s built-in atan and abs

matlab
%% HW3 – Problem 5: CORDIC implementation (10 steps)
% Computes:
%  (a) arctan(1.1)
%  (b) abs(1.2 + 2.1j)
% using 10-step CORDIC in vectoring mode.

clear; clc;

N = 10;                 % number of CORDIC iterations

%% (a) CORDIC approximation of arctan(1.1)

y_over_x = 1.1;
[theta_cordic, xN_atan, yN_atan] = cordic_arctan_ratio(y_over_x, N);

theta_true = atan(y_over_x);

fprintf('--- Part (a): arctan(1.1) ---\n');
fprintf('CORDIC (10 steps): theta = %.12f rad\n', theta_cordic);
fprintf('True value       : theta = %.12f rad\n', theta_true);
fprintf('Absolute error   : %.12e rad\n', theta_cordic - theta_true);
fprintf('Degrees (CORDIC) : %.6f deg\n', theta_cordic * 180/pi);
fprintf('Degrees (true)   : %.6f deg\n\n', theta_true * 180/pi);

%% (b) CORDIC approximation of abs(1.2 + 2.1j)

x0 = 1.2;              % real part
y0 = 2.1;              % imaginary part

[r_cordic, theta_vec, xN_vec, yN_vec] = cordic_polar(x0, y0, N);

r_true     = abs(1.2 + 2.1j);
theta_true_vec = atan2(y0, x0);

fprintf('--- Part (b): abs(1.2 + 2.1j) ---\n');
fprintf('CORDIC (10 steps): r     = %.12f\n', r_cordic);
fprintf('True value       : r     = %.12f\n', r_true);
fprintf('Absolute error   : %.12e\n', r_cordic - r_true);
fprintf('Relative error   : %.12e\n\n', (r_cordic - r_true)/r_true);

fprintf('CORDIC angle     : theta = %.12f rad (%.6f deg)\n', ...
    theta_vec, theta_vec * 180/pi);
fprintf('True angle       : theta = %.12f rad (%.6f deg)\n', ...
    theta_true_vec, theta_true_vec * 180/pi);

%% --- Local functions ---

function [theta, xN, yN] = cordic_arctan_ratio(y_over_x, N)
% CORDIC in vectoring mode to compute atan(y_over_x)
% We choose x0 = 1, y0 = y_over_x, so theta = atan(y0/x0) = atan(y_over_x).

    % Precompute elementary angles alpha_i = atan(2^-i)
    alpha = atan(2.^-(0:N-1));

    x = 1;              % x0
    y = y_over_x;       % y0
    z = 0;              % z0

    for i = 0:N-1
        if y >= 0
            d = 1;
        else
            d = -1;
        end

        % Pseudo-rotation in vectoring mode
        x_new = x + d * y * 2^(-i);
        y_new = y - d * x * 2^(-i);
        z_new = z + d * alpha(i+1);

        x = x_new;
        y = y_new;
        z = z_new;
    end

    theta = z;  % approximate atan(y_over_x)
    xN = x;
    yN = y;
end

function [r, theta, xN, yN] = cordic_polar(x0, y0, N)
% CORDIC in vectoring mode to convert (x0, y0) to (r, theta)
% Output:
%   r     ~ sqrt(x0^2 + y0^2)
%   theta ~ atan2(y0, x0)

    alpha = atan(2.^-(0:N-1));

    x = x0;
    y = y0;
    z = 0;

    for i = 0:N-1
        if y >= 0
            d = 1;
        else
            d = -1;
        end

        x_new = x + d * y * 2^(-i);
        y_new = y - d * x * 2^(-i);
        z_new = z + d * alpha(i+1);

        x = x_new;
        y = y_new;
        z = z_new;
    end

    % CORDIC gain correction: K_N = prod(1/sqrt(1 + 2^(-2i)))
    K = prod(1 ./ sqrt(1 + 2.^(-2*(0:N-1))));

    r     = K * x;    % correct the scale to recover magnitude
    theta = z;        % accumulated angle
    xN    = x;
    yN    = y;
end

You can paste this into a file like hw3_problem5_cordic.m, run it in MATLAB, and it will print the CORDIC results and their errors compared to MATLAB’s exact functions.

We’re solving Problem 3 from your PDF: (HW3_2025.pdf)

A filter has the transfer function
$H(z) = 1 + 2z^{-1} + 4z^{-2} + 8z^{-3}$.
Determine the impulse response of the filter with the modified frequency response
$F(\omega) = H(\omega - 3\pi/8)$.

---

1. Original filter in time and frequency domains

The transfer function of an LTI discrete-time system is

$$H(z) = \sum_{n=-\infty}^{\infty} h[n]\,z^{-n}.$$

Given

$$H(z) = 1 + 2z^{-1} + 4z^{-2} + 8z^{-3},$$

we can directly read off the impulse response $h[n]$:

$$\begin{aligned} h[0] &= 1, \\ h[1] &= 2, \\ h[2] &= 4, \\ h[3] &= 8, \\ h[n] &= 0 \quad \text{for all other } n. \end{aligned}$$

The frequency response of this system is the DTFT of $h[n]$:

$$\begin{aligned} H(\omega) &= \sum_{n=0}^{3} h[n] e^{-j\omega n} \\ &= 1 + 2e^{-j\omega} + 4e^{-j2\omega} + 8e^{-j3\omega}. \end{aligned}$$

---

2. Modified frequency response $F(\omega) = H(\omega - 3\pi/8)$

We’re told the new system has frequency response

$$F(\omega) = H\big(\omega - 3\pi/8\big).$$

So the new system is just the old system, but with its frequency response shifted in frequency by $\omega_0 = 3\pi/8$.

To find the impulse response $f[n]$ of this new filter, we use a standard DTFT property.

---

3. Frequency-shift property of the DTFT

Recall the DTFT pair:

$$x[n] \xleftrightarrow{\text{DTFT}} X(\omega).$$

The frequency-shift (modulation) property says:

$$e^{j\omega_0 n} x[n] \xleftrightarrow{\text{DTFT}} X(\omega - \omega_0).$$

In words:

Multiplying a time-domain signal by $e^{j\omega_0 n}$ shifts its spectrum to the right by $\omega_0$.

Apply this to our case:

• Original: $h[n] \xleftrightarrow{\text{DTFT}} H(\omega)$.
• We want a new sequence $f[n]$ whose DTFT is $F(\omega) = H(\omega - 3\pi/8)$.

By the property above, we simply take

$$f[n] = e^{j(3\pi/8)n} h[n].$$

Then

$$\text{DTFT}\{f[n]\} = H(\omega - 3\pi/8) = F(\omega),$$

exactly as required.

---

4. Explicit impulse response $f[n]$

We already know

$$h[0]=1,\; h[1]=2,\; h[2]=4,\; h[3]=8,\; h[n]=0 \text{ otherwise}.$$

So

$$f[n] = e^{j(3\pi/8)n} h[n].$$

Compute term by term:

For $n = 0$:

$$f[0] = e^{j(3\pi/8)\cdot 0} \cdot h[0] = 1 \cdot 1 = 1.$$

For $n = 1$:

$$\begin{aligned} f[1] &= e^{j(3\pi/8)\cdot 1} \cdot 2 \\ &= 2 e^{j3\pi/8} \\ &= 2\big(\cos(3\pi/8) + j\sin(3\pi/8)\big). \end{aligned}$$

Numerically,

$$\cos(3\pi/8) \approx 0.382683,\quad \sin(3\pi/8) \approx 0.923880,$$

so

$$f[1] \approx 0.76537 + j1.84776.$$

For $n = 2$:

$$\begin{aligned} f[2] &= e^{j(3\pi/8)\cdot 2} \cdot 4 \\ &= 4 e^{j3\pi/4} \\ &= 4\big(\cos(3\pi/4) + j\sin(3\pi/4)\big). \end{aligned}$$

We know $\cos(3\pi/4) = -\sqrt{2}/2$, $\sin(3\pi/4) = \sqrt{2}/2$, so

$$f[2] = 4\left(-\frac{\sqrt{2}}{2} + j\frac{\sqrt{2}}{2}\right) = -2\sqrt{2} + j2\sqrt{2} \approx -2.82843 + j2.82843.$$

For $n = 3$:

$$\begin{aligned} f[3] &= e^{j(3\pi/8)\cdot 3} \cdot 8 \\ &= 8 e^{j9\pi/8} \\ &= 8\big(\cos(9\pi/8) + j\sin(9\pi/8)\big). \end{aligned}$$

Note that $9\pi/8 = \pi + \pi/8$, so

$$\cos(9\pi/8) = -\cos(\pi/8), \quad \sin(9\pi/8) = -\sin(\pi/8).$$

Hence

$$f[3] = -8\cos(\pi/8) - j8\sin(\pi/8).$$

Numerically,

$$\cos(\pi/8) \approx 0.980785,\quad \sin(\pi/8) \approx 0.195090,$$

so

$$f[3] \approx -7.39104 - j1.56072$$

(if you compute with higher precision you might get slightly different last digits; another equivalent approximate value is $-7.39104 - j3.06147$ depending on rounding).

(You don’t really need the decimal values in an exam — the exponential form is usually enough.)

---

5. Final answer

The impulse response of the filter with modified frequency response
$F(\omega) = H(\omega - 3\pi/8)$ is

$$\boxed{ f[n] = e^{j(3\pi/8)n} h[n] }$$

where $h[n]$ is the original impulse response:

$$h[0]=1,\; h[1]=2,\; h[2]=4,\; h[3]=8,\; h[n]=0\ \text{otherwise}.$$

Explicitly,

$$f[n] = \begin{cases} 1, & n=0,\\[4pt] 2e^{j3\pi/8}, & n=1,\\[4pt] 4e^{j3\pi/4}, & n=2,\\[4pt] 8e^{j9\pi/8}, & n=3,\\[4pt] 0, & \text{otherwise}. \end{cases}$$

This sequence is the impulse response of the new, frequency-shifted filter.

Here’s MATLAB code for Problem 2 and Problem 4 as two separate pieces, each building only the plots that the assignment asks for (impulse and frequency responses). (HW3_2025.pdf)

---

✅ Problem 2 — Analog vs Digital (bilinear transform, Fs = 20 Hz)

Assuming the analog transfer function is

$$H(s)=\frac{s+2}{s^2 + s + 1}$$

and the digital version is obtained by the bilinear transform with sampling frequency Fs = 20 Hz.

matlab
%% Problem 2 – Impulse and Frequency Responses (Analog vs Digital)

clear; clc; close all;

%% Analog filter H(s) = (s + 2) / (s^2 + s + 1)
numA = [1 2];      % numerator in s
denA = [1 1 1];    % denominator in s

Fs  = 20;          % sampling frequency [Hz]
T   = 1/Fs;

%% (a) Impulse response (analog vs digital)

% --- Analog impulse response (continuous time) ---
t_end = 3;                          % seconds (choose as you like)
t     = linspace(0, t_end, 2000);   % time vector for analog response
[hA, tA] = impulse(numA, denA, t);  % analog impulse response

% --- Digital filter via bilinear transform ---
% Bilinear transform: s = 2*Fs * (1 - z^-1) / (1 + z^-1)
[bz, az] = bilinear(numA, denA, Fs);   % digital filter coefficients

N  = 60;                         % number of digital samples to show
[hd, n] = impz(bz, az, N);       % digital impulse response

% Plot impulse responses
figure;
subplot(2,1,1);
plot(tA, hA, 'LineWidth', 1.5);
grid on;
xlabel('Time t [s]');
ylabel('h_A(t)');
title('Analog Impulse Response');

subplot(2,1,2);
stem(n*T, hd, 'filled');         % plot vs time in seconds for comparison
grid on;
xlabel('Time t = nT [s]');
ylabel('h_D[n]');
title('Digital Impulse Response (Bilinear, Fs = 20 Hz)');

%% (b) Frequency response (analog vs digital)

% --- Analog frequency response using freqs (continuous-time) ---
% Angular frequency vector (rad/s)
wA = linspace(0, 2*pi*Fs, 2048);       % up to some convenient max freq
HA = freqs(numA, denA, wA);            % analog frequency response
fA = wA / (2*pi);                      % convert rad/s to Hz

% --- Digital frequency response using freqz ---
% H_D(e^{j\omega}) over [0, Fs/2] Hz
[HD, fD] = freqz(bz, az, 2048, Fs);    % fD is in Hz

% Plot magnitude responses (in dB)
figure;
subplot(2,1,1);
plot(fA, 20*log10(abs(HA) + eps), 'LineWidth', 1.5);
grid on;
xlabel('Frequency [Hz]');
ylabel('|H_A(j\omega)| [dB]');
title('Analog Magnitude Response');

subplot(2,1,2);
plot(fD, 20*log10(abs(HD) + eps), 'LineWidth', 1.5);
grid on;
xlabel('Frequency [Hz]');
ylabel('|H_D(e^{j\omega})| [dB]');
title('Digital Magnitude Response (Bilinear, Fs = 20 Hz)');

You can put this as a single script (e.g., problem2_plots.m). It:

• Builds the analog filter from numA, denA
• Applies bilinear to get the digital filter
• Plots analog & digital impulse responses
• Plots analog & digital magnitude responses (in dB)

---

✅ Problem 4 — Discrete-time System $H(z) = \dfrac{z + 1}{z^3 + 2z^2 + 3z + 4}$

We rewrite $H(z)$ in terms of powers of $z^{-1}$:

$$H(z) = \frac{z + 1}{z^3 + 2z^2 + 3z + 4} = \frac{z^{-2} + z^{-3}}{1 + 2z^{-1} + 3z^{-2} + 4z^{-3}}.$$

So the discrete-time filter coefficients (for use with impz, freqz, filter) are:

• Numerator: $b[0..3] = [0,\ 0,\ 1,\ 1]$
• Denominator: $a[0..3] = [1,\ 2,\ 3,\ 4]$

Here is the MATLAB code to plot the impulse and frequency responses.

matlab
%% Problem 4 – Impulse and Frequency Responses of H(z) = (z + 1)/(z^3 + 2z^2 + 3z + 4)

clear; clc; close all;

% Given transfer function:
% H(z) = (z + 1) / (z^3 + 2 z^2 + 3 z + 4)
% Rewrite in powers of z^-1:
% H(z) = (z^-2 + z^-3) / (1 + 2 z^-1 + 3 z^-2 + 4 z^-3)

b = [0 0 1 1];   % numerator coefficients in z^-1
a = [1 2 3 4];   % denominator coefficients in z^-1

%% (a) Impulse response

N = 50;                      % number of samples of impulse response
[h, n] = impz(b, a, N);      % discrete-time impulse response

figure;
stem(n, h, 'filled');
grid on;
xlabel('n');
ylabel('h[n]');
title('Impulse Response of H(z)');

%% (b) Frequency response

Nfft = 1024;                 % FFT points for freqz
[H, w] = freqz(b, a, Nfft);  % w in rad/sample

figure;
subplot(2,1,1);
plot(w/pi, abs(H), 'LineWidth', 1.5);
grid on;
xlabel('Normalized Frequency (\omega/\pi)');
ylabel('|H(e^{j\omega})|');
title('Magnitude Response of H(z)');

subplot(2,1,2);
plot(w/pi, angle(H), 'LineWidth', 1.5);
grid on;
xlabel('Normalized Frequency (\omega/\pi)');
ylabel('Phase [rad]');
title('Phase Response of H(z)');

You can save this as something like problem4_plots.m.

If you also need the difference equation or block diagrams (Direct-Form I / II) later, I can write those out too, but this code already produces all the plots the assignment asks for in Problem 4.

We need to design a Chebyshev lowpass digital filter using the impulse invariance method, with specs: (HW3_2025.pdf)

• Passband edge: $f_p = 18\,\text{MHz}$
• Stopband edge: $f_s = 20\,\text{MHz}$
• Passband ripple: $R_p = 0.3\ \text{dB}$
• Stopband attenuation: $R_s = 60\ \text{dB}$
• Sampling frequency: $F_s = 50\,\text{MHz}$

Then:

1. Plot impulse responses of analog and digital filters.
2. Plot magnitude responses of analog and digital filters in the frequency domain.
3. Provide MATLAB code.

---

1. Step 1 – Design the analog Chebyshev lowpass filter

We first design a continuous-time (analog) Chebyshev Type I lowpass filter that meets the given specs.

1.1. Convert frequency specs to rad/s

$$\begin{aligned} F_s &= 50\cdot 10^6\ \text{Hz}, \quad T = \frac{1}{F_s} \\ f_p &= 18\cdot 10^6\ \text{Hz}, \quad f_s = 20\cdot 10^6\ \text{Hz} \\ \Omega_p &= 2\pi f_p \\ \Omega_s &= 2\pi f_s \end{aligned}$$

These $\Omega_p, \Omega_s$ are used in the analog filter design.

1.2. Chebyshev Type I parameters

Passband ripple $R_p$ and stopband attenuation $R_s$:

$$\begin{aligned} \epsilon &= \sqrt{10^{R_p/10} - 1} \\ D_1 &= 10^{R_p/10} - 1 = \epsilon^2 \\ D_2 &= 10^{R_s/10} - 1 \end{aligned}$$

The minimum filter order $N$ for a Chebyshev Type I lowpass is

$$N \ge \frac{\operatorname{acosh}\left( \sqrt{\dfrac{D_2}{D_1}} \right)} {\operatorname{acosh}\left( \dfrac{\Omega_s}{\Omega_p} \right)}$$

For the given numbers, this gives approximately

$$N \approx 19.1 \Rightarrow N = 20$$

(in practice, MATLAB’s cheb1ord will compute the exact minimum integer order.)

In MATLAB, we don’t have to compute these by hand; we can use:

matlab
[Na, Wp_analog] = cheb1ord(Wp, Ws, Rp, Rs, 's');

This returns the required analog filter order Na (which should be about 20) and an adjusted analog passband edge Wp_analog.

1.3. Analog transfer function $H_a(s)$

We then design the analog Chebyshev lowpass transfer function:

$$H_a(s) = \frac{B_a(s)}{A_a(s)}$$

with MATLAB:

matlab
[ba, aa] = cheby1(Na, Rp, Wp_analog, 'low', 's');

• ba and aa are vectors of numerator and denominator coefficients of $H_a(s)$.
• This filter meets the analog specs at 18–20 MHz.

---

2. Step 2 – Impulse invariance: analog → digital

The impulse invariance method defines the digital impulse response as:

$$h_d[n] = T\, h_a(nT), \quad n = 0, 1, 2, \dots$$

where:

• $h_a(t)$ is the analog impulse response of $H_a(s)$,
• $T = 1/F_s$ is the sampling period.

If the analog transfer function has partial fraction expansion

$$H_a(s) = \sum_{k} \frac{R_k}{s - p_k},$$

then:

• Analog impulse response: $$h_a(t) = \sum_{k} R_k e^{p_k t} u(t)$$
• Digital impulse response: $$h_d[n] = T \sum_{k} R_k e^{p_k nT} u[n]$$

Taking the Z-transform:

$$H_d(z) = \sum_{n=0}^{\infty} h_d[n] z^{-n} = \sum_k \frac{T R_k}{1 - e^{p_k T} z^{-1}}$$

So every analog pole $p_k$ maps to a digital pole at $z = e^{p_k T}$, preserving the shape of the impulse response between 0 and Nyquist (but causing aliasing of higher-frequency components).

In MATLAB, this mapping is implemented by:

matlab
[bd, ad] = impinvar(ba, aa, Fs);

This directly produces the digital filter:

$$H_d(z) = \frac{B_d(z)}{A_d(z)} = \frac{\sum_k b_{d,k} z^{-k}}{\sum_k a_{d,k} z^{-k}}.$$

---

3. Step 3 – Impulse responses (analog and digital)

• Analog impulse response $h_a(t)$ is obtained by:

  matlab
  [ha, t] = impulse(ba, aa, t_vec);

• Digital impulse response $h_d[n]$ is obtained by:

  matlab
  [hd, n] = impz(bd, ad, N_imp);

If designed correctly, you should see that

$$h_d[n] \approx T\, h_a(nT)$$

for the plotted time window (this is the core idea of impulse invariance).

---

4. Step 4 – Magnitude responses (analog and digital)

4.1. Analog frequency response

For the analog filter, the frequency response is

$$H_a(j\Omega) = H_a(s)\big|_{s = j\Omega}$$

We compute it over a set of angular frequencies $\Omega$ (rad/s):

matlab
wA = linspace(0, 2*pi*25e6, 2000);   % up to Nyquist (25 MHz)
HA = freqs(ba, aa, wA);              % analog frequency response
fA = wA / (2*pi);                    % convert to Hz

We can then plot $|H_a(j\Omega)|$ in linear scale or dB.

4.2. Digital frequency response

The digital filter’s DTFT is $H_d(e^{j\omega})$. Using freqz with sampling frequency Fs:

matlab
[HD, fD] = freqz(bd, ad, 2048, Fs);  % fD in Hz, range 0..Fs/2

This computes the digital magnitude response on $[0, F_s/2]$, directly comparable with the analog response.

---

5. Complete MATLAB code for Problem 1

This script:

• Designs the analog Chebyshev lowpass.
• Converts it to digital using impulse invariance.
• Plots impulse responses for analog and digital.
• Plots magnitude responses for analog and digital on a frequency (Hz/MHz) axis.

matlab
%% Problem 1 – Chebyshev LPF via Impulse Invariance
% Specs:
%   Passband edge: fp = 18 MHz
%   Stopband edge: fs = 20 MHz
%   Passband ripple: Rp = 0.3 dB
%   Stopband attenuation: Rs = 60 dB
%   Sampling freq: Fs = 50 MHz
%
% Tasks:
%   (a) Plot analog and digital impulse responses
%   (b) Plot analog and digital magnitude responses

clear; clc; close all;

%% Given specifications
Fp = 50e6;          % sampling frequency [Hz]
T  = 1/Fp;          % sampling period [s]

fp = 18e6;          % analog passband edge [Hz]
fs = 20e6;          % analog stopband edge [Hz]

Rp = 0.3;           % passband ripple [dB]
Rs = 60;            % stopband attenuation [dB]

% Convert to rad/s for analog design
Wp = 2*pi*fp;       % passband edge [rad/s]
Ws = 2*pi*fs;       % stopband edge [rad/s]

%% Step 1: Analog Chebyshev Type I lowpass design

% Find minimum order and passband edge for analog Chebyshev filter
[Na, Wp_analog] = cheb1ord(Wp, Ws, Rp, Rs, 's');

fprintf('Analog Chebyshev Type I order Na = %d\n', Na);

% Design analog lowpass filter H_a(s)
[ba, aa] = cheby1(Na, Rp, Wp_analog, 'low', 's');

%% Step 2: Impulse invariance to obtain digital filter

% Digital filter via impulse invariance
[bd, ad] = impinvar(ba, aa, Fp);

%% (a) Impulse responses: analog vs digital

% Choose a time window for impulse response comparison
N_imp = 200;                       % number of samples (discrete-time)
t_imp = (0:N_imp-1)*T;            % corresponding continuous times

% Analog impulse response h_a(t)
[ha, t_cont] = impulse(ba, aa, t_imp);  % use same time grid as digital

% Digital impulse response h_d[n]
[hd, n] = impz(bd, ad, N_imp);    % n = 0,1,...,N_imp-1

% Plot impulse responses
figure;
subplot(2,1,1);
plot(t_cont*1e6, ha, 'LineWidth', 1.5);
grid on;
xlabel('Time t [\mus]');
ylabel('h_a(t)');
title('Analog Impulse Response (Chebyshev LPF)');

subplot(2,1,2);
stem(n*T*1e6, hd, 'filled');
grid on;
xlabel('Time t = nT [\mus]');
ylabel('h_d[n]');
title('Digital Impulse Response (Impulse Invariance)');

%% (b) Magnitude responses: analog vs digital

% --- Analog magnitude response ---
fA_max = 25e6;                             % plot up to Nyquist of digital system
wA = linspace(0, 2*pi*fA_max, 2048);       % angular frequency [rad/s]
HA = freqs(ba, aa, wA);
fA = wA/(2*pi);                            % Hz

% --- Digital magnitude response ---
[HD, fD] = freqz(bd, ad, 2048, Fp);        % fD in Hz, 0..Fs/2

% Convert to MHz for nicer axis labels
fA_MHz = fA / 1e6;
fD_MHz = fD / 1e6;

figure;

subplot(2,1,1);
plot(fA_MHz, 20*log10(abs(HA) + eps), 'LineWidth', 1.5);
grid on;
xlabel('Frequency [MHz]');
ylabel('|H_a(j\Omega)| [dB]');
title('Analog Chebyshev LPF Magnitude Response');
xline(fp/1e6, '--', 'f_p', 'LabelVerticalAlignment','bottom');
xline(fs/1e6, '--', 'f_s', 'LabelVerticalAlignment','bottom');

subplot(2,1,2);
plot(fD_MHz, 20*log10(abs(HD) + eps), 'LineWidth', 1.5);
grid on;
xlabel('Frequency [MHz]');
ylabel('|H_d(e^{j\omega})| [dB]');
title('Digital LPF Magnitude Response (Impulse Invariance)');
xline(fp/1e6, '--', 'f_p', 'LabelVerticalAlignment','bottom');
xline(fs/1e6, '--', 'f_s', 'LabelVerticalAlignment','bottom');

%% Optional: check the impulse-invariance relation h_d[n] ≈ T * h_a(nT)
% err = hd - T*ha(:);
% fprintf('Max abs error between hd and T*ha(nT): %e\n', max(abs(err)));

What this code does

• Designs the analog Chebyshev LPF using cheb1ord + cheby1 with the given passband/stopband specs.
• Converts that analog filter to a digital IIR filter using impinvar (impulse invariance).
• Plots:
  • Analog impulse response $h_a(t)$ and digital impulse response $h_d[n]$.
  • Analog and digital magnitude responses on the same frequency scale (in MHz), so you can visually verify that in the band up to Nyquist (0–25 MHz), the digital magnitude closely follows the analog one.

You can put this into problem1_impulse_invariance.m or similar and run it in MATLAB.