一、目的
掌握信号处理的基本思想，理解采样信号的频谱特性,加强信号采样与重建的有关基本概念的理解，深入理解线性时不变系统输出与输入的关系，了解数字信号采样率转换前后信号频谱的特征。
二、内容与设计思想
1、给定序列，绘制其图像并分析其混叠。
设计思路：
首先定义一个由两个不同频率正弦波组成的序列，然后绘制这个序列的图像。接着，将这个序列扩展为两个不同周期的版本，分别是16周期和10周期，并通过绘图来观察这些扩展序列。最后，通过分析这些序列的频谱来探讨混叠现象。
代码分析：
结构包括初始化序列、绘制原始序列、实现周期性扩展、绘制扩展序列。在初始化阶段，定义了时间序列和原始序列的值。周期性扩展是通过循环和mod函数实现的，它将原始序列的值重复以形成新的周期性序列。
混叠分析：
式子①由两个正弦波组成，一个频率为f1 = 0.125Hz，另一个频率为f2 = 0.25 Hz。根据奈奎斯特定理，采样频率应至少为最高频率的两倍，即2*0.25 = 0.5Hz。由于采样频率为1 Hz，高于0.5 Hz，因此原始序列不会发生混叠。
式子②被扩展为16周期，因此可以看成以1/16Hz的频率采样原始序列。由于原始序列的最高频率为0.25 Hz，而周期性扩展的等效采样频率为1/16Hz，远低于0.5 Hz，这可能导致混叠。
式子③被扩展为10周期，可以视为以1/10Hz 的频率采样原始序列，与16周期扩展类似，10周期扩展的等效采样频率为1/10Hz，也远低于0.5 Hz，这同样可能导致混叠，且更有可能与低频分量重叠。
输出结果截图：
2、设计序列的移位和累加函数
设计思路：
根据给定的原始序列x(n)，计算并绘制两个新的序列 x1(n)和 x2(n)。x1(n)是通过原始序列 x(n)的移位、累加和缩放操作得到的，而x2(n)则是基于x(n)的线性组合与累加求和。由于n仅四位，因此要限制n的范围为0到3。
主要代码分析：
n1 = 0:3 //设定 x1(n)的索引范围
x1(i) = 2x(mod(n1(i)+2, 4)+1) - x(mod(n1(i)-1, 4)+1) - 2x(n1(i)+1);
//根据题目的公式计算 x1(n)x1(n) 的每个元素，并使用 mod 函数处理可能的越界索引
x2(i) = x2(i) + n2(i) * x(mod(n2(i)-k, 4)+1);
//循环累加，计算 x2(n)的每个元素，并同样使用 mod 函数处理越界索引
输出结果截图：
由图可知，
序列 x1(n)：
x1(0)=-1；x1(1)=11；x1(2)=-3；x1(3)=-15
序列 x2(n)：
x2(0)=0；x2(1)=9；x2(2)=14；x2(3)=33
3、绘制因果系统的幅频响应和相频响应
设计思路：
首先定义给定因果系统的分子和分母系数分别为b和a，这些系数将用于描述系统的传递函数，再使用freqz函数来计算系统的频率响应，将计算得到的频率响应绘制成幅频响应和相频响应图。通过观察绘制的图形，分析系统的频率特性。
代码分析：
b = [1, sqrt(2), 1];
a = [1, -0.67, 0.9];
//定义了系统的分子和分母系数，这些系数直接从系统的传递函数中提取。
[H, w] = freqz(b, a, ‘half’, 1000);
//使用freqz函数计算频率响应。half参数指定只计算从0到π的频率响应，1000指定计算1000个频率点。
plot(w/pi, 20log10(abs(H)));
//绘制幅频响应，20log10(abs(H))将幅值转换为分贝单位，以便于观察系统的增益特性。
plot(w/pi, unwrap(angle(H)));
//绘制相频响应，unwrap(angle(H))用于处理相位计算中的不连续性，确保相位变化是平滑的。
输出结果截图：
幅频响应分析
在低频段，幅频响应接近0 dB，表明系统在低频段对信号的增益接近1，即信号几乎不受影响。在大约0.4到0.5的归一化频率处，幅频响应开始显著下降，这可能是系统的截止频率区域，表明系统开始衰减高频信号。在高频段（接近0.8），幅频响应迅速下降到-300 dB以下，这表明系统在高频段对信号有很强的衰减作用。在0.7到0.8之间有一个突然的下降，这可能是由于系统的极点导致的不稳定行为，或者是数值计算的不稳定性。
相频响应分析
在低频段，相位接近0弧度，随着频率的增加，相位逐渐下降，表明信号在通过系统时会有一定的相位延迟。在大约0.7处，相位有一个突然的跳变，这可能是系统在该频率附近有一个极点或零点，导致相位特性发生突变。相位响应在-3到0弧度之间变化，这表明系统对信号的相位影响主要表现为相位延迟。
4、
设计思路：
首先创建一个复数序列，实部和虚部在指定区间内均匀分布，再使用FFT算法计算原始序列的DTFT，通过翻转序列来模拟时间域的反转。再次使用FFT算法计算折叠后序列的DTFT。
最后绘制并比较原始序列和折叠后序列的DTFT幅度响应。
代码分析：
x = rand(1, length(n)) + 1i * rand(1, length(n));
//生成了一个复随机序列x，其实部和虚部在 [0,1] 区间内均匀分布。
y = fft(x);
//使用fft函数计算序列x的DTFT。
y_flipped = fft(fliplr(x));
//使用fliplr函数将序列x左右翻转，然后计算折叠后序列的DTFT。
plot(abs(y));
plot(abs(y_flipped));
//绘制图形
输出结果分析：
两个结果图展示了原始序列和折叠后序列的DTFT幅度响应。可以看到，两个图的幅度响应在低频区域有较高的幅度，随着频率的增加，幅度逐渐减小。这是典型的随机序列的频谱特性。尽管两图的幅度响应在视觉上看起来相似，但折叠特性应该表现为频率域的对称性。在本实验中，由于使用了随机序列，因此这种对称性可能不太明显。
5、绘制给定系统的稳态响应
设计思路：
本实验模拟并分析一个线性时不变系统的稳态响应。系统由一阶差分方程 y(n)=0.8y(n-1)+x(n)描述，其中x(n)是输入信号，y(n)是输出信号。输入信号为余弦波 x(n)=cos(0.05πn)u(n)，其中 u(n)是单位阶跃函数。实验目的是计算系统的输出并绘制输入和输出信号的图形。
代码分析：
x = cos(0.05 * pi * n) .* (n >= 0);
//使用余弦函数生成周期性信号，并通过 (n >= 0) 确保信号从 n=0开始。
y(i) = a * y(i-1) + x(i);
//通过迭代差分方程计算输出信号。从 n=2开始，因为y(1) 需要 y(0)的值。
输出结果截图：
上图显示了周期性余弦波输入信号x(n)，其幅度在[−1,1]之间，周期性地变化。下图显示了系统的输出信号y(n)。由于系统的反馈系数a=0.8小于1，系统是稳定的，输出信号在输入信号的基础上逐渐达到稳态，显示出衰减的振荡特性。输出信号的幅度随着时间的推移而减小，这是由于反馈系数小于1导致的。在足够长的时间后，输出信号将趋于一个稳定的周期性波形，其幅度和相位将受到系统特性和输入信号的影响。
三、环境
Matlab
四、小结
1、实验中遇到的问题及解决过程
(1)最初可能不理解为什么在特定的采样率下会发生混叠。
解决方法：具体计算采样频率和信号最高频率的比值，确认是否满足奈奎斯特准则，再绘制不同采样率下的频谱图，直观比较混叠现象。
(2)在实现信号的周期性扩展时，难以正确地重复信号值。
解决方法：使用repmat函数来重复信号值，创建周期性扩展的序列，确保在重复信号时，正确地处理了信号的起始和结束索引，避免不必要的重复或遗漏。
2、实验中产生的错误及原因分析
(1)在实现信号移位和累加函数时，发现数组越界。
原因分析：没有正确处理数组的边界条件，导致越界错误。通过添加边界检查和使用 mod 函数可以解决这个问题。
3、实验体会和收获
本次实验中，我通过动手实践，加深了对数字信号处理核心概念的理解，特别是在信号采样、混叠现象、系统频率响应分析以及线性时不变系统稳定性等方面。通过编写和调试MATLAB代码，我不仅提升了编程技能，还学会了如何运用数学工具来分析和解决实际问题。实验过程中，我通过观察不同采样率下的信号混叠现象，直观地理解了奈奎斯特采样定理的重要性。同时，通过分析系统对不同输入信号的响应，我更加熟悉了如何使用MATLAB进行信号处理和系统分析。这些经验对于我未来在信号处理领域的学习和研究工作具有重要的价值。
源代码附件
1、
n = 0:100;
x = 3cos(0.125pin + 0.2pi) + 2sin(0.25pin + 0.1pi);
figure;
plot(n, x);
title(‘x(n)’);
xlabel(‘n’);
ylabel(‘x(n)’);
n1 = 16;
x1 = repmat(x(1:n1), 1, 4);
figure;
plot(0:4n1-1, x1);
title(‘周期序列 x_{16}(n) 周期为16’);
xlabel(‘n’);
ylabel(‘x_{16}(n)’);
n2 = 10;
x2 = repmat(x(1:n2), 1, 4);
figure;
plot(0:4n2-1, x2);
title(‘周期序列 x_{10}(n) 周期为10’);
xlabel(‘n’);
ylabel(‘x_{10}(n)’);
2、
n = 0:3;
x = [1, -1, 3, 5];
n1 = 0:3;
x1 = zeros(size(n1));
for i = 1:4
x1(i) = 2x(mod(n1(i)+2, 4)+1) - x(mod(n1(i)-1, 4)+1) - 2x(n1(i)+1);
end
n2 = 0:3;
x2 = zeros(size(n2));
for i = 1:4
for k = 1:5
x2(i) = x2(i) + n2(i) * x(mod(n2(i)-k, 4)+1);
end
end
disp(['x1(n) = ', num2str(x1)]);
disp(['x2(n) = ', num2str(x2)]);
figure;
subplot(2, 1, 1);
stem(n1, x1);
title(‘x1(n)’);
xlabel(‘n’);
ylabel(‘x1(n)’);
subplot(2, 1, 2);
stem(n2, x2);
title(‘x2(n)’);
xlabel(‘n’);
ylabel(‘x2(n)’);
3、
b = [1, sqrt(2), 1];
a = [1, -0.67, 0.9];
[H, w] = freqz(b, a, ‘half’, 1000);
figure;
subplot(2,1,1);
plot(w/pi, 20*log10(abs(H)));
title(‘幅频响应’);
xlabel(‘归一化频率’);
ylabel(‘幅度’);
grid on;
subplot(2,1,2);
plot(w/pi, unwrap(angle(H)));
title(‘相频响应’);
xlabel(‘归一化频率’);
ylabel(‘相位’);
grid on;
4、
n = -10:20;
x = rand(1, length(n)) + 1i * rand(1, length(n));
y = fft(x);
y_flipped = fft(fliplr(x));
figure;
subplot(2,1,1);
plot(abs(y));
title(‘原始序列的DTFT幅度响应’);
xlabel(‘频率点’);
ylabel(‘幅度’);
subplot(2,1,2);
plot(abs(y_flipped));
title(‘折叠后序列的DTFT幅度响应’);
xlabel(‘频率点’);
ylabel(‘幅度’);
5、
a = 0.8;
n = 0:99;
x = cos(0.05 * pi * n) .* (n >= 0);
y = zeros(1, 100);
for i = 2:100
y(i) = a * y(i-1) + x(i);
end
figure;
subplot(2,1,1);
plot(n, x);
title(‘x(n)’);
xlabel(‘n’);
ylabel(‘x(n)’);
subplot(2,1,2);
plot(n, y);
title(‘y(n)’);
xlabel(‘n’);
ylabel(‘y(n)’);