2-20 线性因果系统用下面的差分方程描述: y(n)-2ry(n-1)cosθ+r2y(n-2)=x(n) 式中,x(n)=anu(n),0<a<1;0<r<1,0=常数。试求系统的输出 y(n)用Matlab代码实现

该线性因果系统是一个二阶差分方程，可以表示为状态空间形式或者通过z变换求解。给定输入 \( x(n) = a^n u(n) \)，其中 \( u(n) \) 是单位样点序列，\( 0 < a < 1 \) 表示衰减系数，\( 0 < r < 1 \) 是系统极点所在的位置，\( θ \) 是极点相对于实轴的角度。

对于这个系统的输出 \( y(n) \)，我们可以先找到其Z变换Y(z)。由于 \( x(n) \) 的Z变换是 \( X(z) = \frac{a}{1 - az} \)，代入差分方程：

\[ Y(z) - 2r z^{-1} \cdot Y(z) \cdot e^{j\theta} + r^2 z^{-2} \cdot Y(z) = X(z) \]

然后我们可以通过求解这个齐次线性微分方程并结合初始条件来得到Y(z)。不过，为了直接给出Matlab代码实现，我们需要先计算出系统函数H(z)（即系统的零极点图），然后利用`filter`函数。

系统函数 H(z) 可以这样表示：
\[ H(z) = \frac{Y(z)}{X(z)} = \frac{1}{1 - 2rz\cos(\theta) + rz^2} \]

在Matlab中，假设已知\( a \), \( r \), 和 \( \theta \)，可以这样编写代码：

% 定义系统参数
a = 0.7; % 衰减系数
r = 0.4; % 极点位置
theta = pi / 4; % 极点角度 (弧度)

% 计算系统函数 H(z)
num = [1, -2*r*cos(theta), r^2];
den = [1, -2*r, r^2];

% 初始化输出和滤波器
y = zeros(1, N); % 输出序列，N是你想要计算的样本数
sys = tf(num, den); % 创建传递函数模型

% 使用filter函数得到输出
for n = 1:N
    if n == 1
        % 初始条件 y(0) = 0
        y(n) = filter(sys, 1, zeros(1, 1));
    else
        y(n) = filter(sys, 1, y(n-1)); % 前一时刻的输出作为当前时刻的输入
    end
end

注意：上述代码中，你需要替换 `N` 成实际需要的采样点数目。另外，如果初始条件不是零，记得相应修改 `y(n=1)` 的值。