超详细指南lab1a：MATLAB 信号采样、频谱分析、非线性削波及混叠现象完全解读

Ultra-Detailed Guide: Signal Sampling, Spectrum Analysis, Nonlinear Clipping, and Aliasing in MATLAB

📌 Article Abstract (Directly Searchable)

This ultra-detailed guide provides a “textbook-level” solution to a classic set of digital signal processing exercises, covering:

1. Ideal Sampled Signal Generation: 1 kHz sine wave, 1s duration, 5 kHz sampling, with rigorous time-base verification.
2. Spectrum Analysis & Resolution: FFT computation, 0Hz-centered bilateral spectrum, and derivation of frequency resolution (( \Delta f = 1\text{Hz} )).
3. Nonlinear Distortion (Hard Clipping): 2.25 kHz sine (amp 1.0) sampled at 5 kHz, hard-clipped to ±0.3, with time-domain waveform inspection.
4. Quantitative Aliasing Analysis: Exact calculation of folded positions for 3rd, 5th, 7th, and 9th harmonics (e.g., 1.75 kHz, 1.25 kHz, etc.) and explanation of spurious peaks in the FFT plot.
5. Oversampling (100 kHz) Demonstration: How high sampling rates prevent aliasing.
6. Auxiliary Experiments: Square wave spectrum and different clipping thresholds.

🛠️ Zero: Mathematical Preliminaries

• Nyquist Theorem: ( f_s > 2 f_{max} ), otherwise aliasing.
• Fourier Series of Hard-Clipped Sine: Produces odd harmonics (3f, 5f, 7f…) with decreasing amplitudes.
• Aliasing Formula: ( f_{alias} = | f_{harm} - k \cdot f_s | ), where ( k ) is an integer chosen to bring the frequency into ( [0, f_s/2] ).

📝 Part 1: Sampled Sine Generation (Code Q1)

Task: 1 kHz sine, 1 s, 5 kHz sampling.

Code Analysis:

• Sampling interval ( T_s = 0.2 ) ms.
• Samples per cycle = 5, satisfying ( f_s = 5000 > 2000 ).
• Time vector0:1/fs:duration-1/fsensures exactly 5000 points without boundary overflow.

Result: Plot shows 1000 complete cycles, each with 5 points.

📊 Part 2: FFT Spectrum & Frequency Resolution (Code Q2)

Task: Compute FFT, center at 0 Hz, calculate resolution.

Code:

fft_signal=fft(signal);f_axis=linspace(-fs/2,fs/2,N);plot(f_axis,fftshift(abs(fft_signal))/N);

Resolution Derivation:
[
\Delta f = \frac{f_s}{N} = \frac{5000}{5000} = 1 \text{ Hz}
]
Interpretation: Adjacent frequency bins are 1 Hz apart.
Spectrum: Two peaks at ±1 kHz with magnitude 0.5 (normalized bilateral spectrum).

✂️ Part 3: Hard Clipping (Code Q3)

Task: 2.25 kHz sine (amp 1.0, fs=5 kHz) clipped to ±0.3.

Code:

signal(signal>0.3)=0.3;signal(signal<-0.3)=-0.3;

Time-domain effect: Peaks are flattened, creating sharp corners. The fundamental (2.25 kHz) does NOT alias by itself (( 2 \times 2.25 = 4.5 < 5 )), but harmonics will.

🎯 Part 4: Clipped Spectrum & Quantitative Aliasing (Code Q4) —Core of this Guide

Continuous-time theory: Hard clipping produces odd harmonics (3f, 5f, 7f…).

Aliasing Calculation Table(with ( f_s = 5 ) kHz, Nyquist = 2.5 kHz):

FFT Code:

X=abs(fft(clipped_signal));f_axis=linspace(-fs/2,fs/2,N);plot(f_axis,fftshift(X)/N);xlim([-3000,3000]);

Actual Result: Peaks at 2.25k, 1.75k, 1.25k, 0.75k, 0.25k Hz — exactly matching the table. This proves thatinsufficient sampling rate causes out-of-band harmonics to fold back, corrupting the in-band spectrum.

🚀 Part 5: Oversampling to the Rescue (Code Q6)

Task: Raise sampling rate to 100 kHz, repeat clipping.

Why it works: Nyquist becomes 50 kHz. Even the 21st harmonic (47.25 kHz) remains below Nyquist. No folding occurs. The spectrum consists of clean, distinct harmonic lines. Oversampling provides guard bands for subsequent digital filtering.

🔲 Part 6: Square Wave Verification (Code Q7)

A 1.1 kHz square wave (5 kHz sampling) also shows aliased harmonics (e.g., 5.5 kHz folds to 0.5 kHz). This confirms aliasing is a general phenomenon, not limited to clipping.

🔬 Part 7: Clipping Threshold Effect (Code Q8)

Lowering the threshold to ±0.1 (sampling at 10 kHz) increases harmonic amplitudes. Although fs is higher, the 3rd harmonic (6.75 kHz) still aliases to 1.75 kHz (since 6.75 > 5 kHz Nyquist). This highlights thatboth clipping severity and sampling rate determine final spectral purity.

💎 Final Summary

📌 文章简介（可直接检索到题目）

本文针对数字信号处理中的一组经典实验任务，提供“教科书级别”的详细解答。内容严格按照以下步骤展开：

1. 生成理想采样信号：创建频率为1 kHz、时长为1 秒、采样率为5 kHz的正弦波，并严格验证时间轴生成的正确性。
2. 频谱分析与分辨率计算：对信号进行FFT，绘制以0 Hz 为中心的双边频谱，并推导频率分辨率（( \Delta f = 1\text{Hz} )）的物理含义。
3. 非线性失真（硬削波）：生成频率为2.25 kHz、幅度为1.0的正弦波（仍用 5 kHz 采样），执行硬削波（限幅至±0.3），绘制时域畸变波形。
4. 削波信号的频谱与混叠定量分析：对削波信号做 FFT。通过数学计算，逐一找出3次、5次、7次、9次谐波因采样率不足而“折叠”到低频区的具体位置（如 1.75 kHz、1.25 kHz 等），并解释为什么实际频谱中会出现这些“虚假”峰值。
5. 过采样对比实验：将采样率提升至100 kHz，展示过采样如何为高次谐波提供足够的频率空间，从而避免混叠。
6. 辅助实验：方波频谱分析（验证奇次谐波普遍性）及不同削波阈值（±0.1）对混叠程度的影响。

本文包含完整的 MATLAB 代码、傅里叶级数理论推导以及工程启示，适合作为数字信号处理课程的深度辅导材料。

🔬 零、预备知识与数学工具

在开始编码前，必须明确两个核心数学工具：

• 奈奎斯特采样定理：若信号最高频率为 ( f_{max} )，则采样率 ( f_s ) 必须满足 ( f_s > 2 f_{max} )，否则发生混叠。
• 傅里叶级数（硬削波）：一个幅度为 ( A )、频率为 ( f ) 的正弦波，经过对称硬削波（限幅至 ( \pm V_{clip} )）后，其频谱包含：
  • 基波（频率 ( f )，幅度降低）；
  • 奇次谐波（( 3f, 5f, 7f, \dots )），幅度随阶数增加而递减。
• 离散频谱混叠公式：若某频率成分 ( f_{harm} ) 高于奈奎斯特频率 ( f_s/2 )，则其在 ( [0, f_s/2] ) 内的混叠位置为：
  [
  f_{alias} = | f_{harm} - k \cdot f_s |
  ]
  其中 ( k ) 是使 ( f_{alias} \leq f_s/2 ) 的整数。

📝 第一部分：生成采样正弦波（对应代码 Q1）

任务描述

创建一个长度为1 秒、频率为1 kHz的正弦信号，采样率为5 kHz。绘制该信号，确保时间轴完全正确，使波形在视觉上清晰无误。

代码实现（逐行解析）

clc;clear all;% 清空工作区与命令窗口frequency=1000;% 信号频率 1000 Hzsampling_rate=5000;% 采样率 5000 Hz（即每秒 5000 个样本）duration=1;% 信号时长 1 秒% 计算总采样点数：5000 * 1 = 5000 个点num_samples=duration*sampling_rate;% 生成时间向量（极其重要！）% 从 0 开始，步长为 1/sampling_rate = 0.2ms，到 duration - 1/sampling_rate 结束% 这样保证了向量长度正好为 5000，且最后一个点不会超出 1 秒边界time=0:1/sampling_rate:duration-1/sampling_rate;% 生成正弦波信号：sin(2*pi*f*t)，幅度默认为 1signal=sin(2*pi*frequency*time);% 绘图figure(1);plot(time,signal);xlabel('Time (s)');ylabel('Amplitude');title('1kHz 正弦信号，采样率 5kHz (Sinusoidal Signal at 1kHz, fs=5kHz)');grid on;% 显示网格，便于观察周期

理论验证与图解

• 采样间隔( T_s = 1/5000 = 0.0002 ) 秒（0.2 ms）。
• 信号周期( T = 1/1000 = 0.001 ) 秒（1 ms）。
• 每个周期的采样点数= ( T / T_s = 0.001 / 0.0002 = 5 ) 个点。
• 因为 ( f_s = 5000 > 2 \times 1000 = 2000 )，严格满足奈奎斯特条件，信号可以被完美重建。
• 绘图结果：在 1 秒内恰好出现 1000 个完整的正弦波周期，每个周期由 5 个采样点勾勒出来。

📊 第二部分：FFT 频谱分析与频率分辨率（对应代码 Q2）

任务描述

绘制该采样信号的理论预期频谱（纸上），然后使用 MATLAB 的fft函数计算频谱，要求以0 Hz 为中心显示，并计算该频谱图的频率分辨率。

代码实现（逐行解析）

figure(2);% 计算 FFT（快速傅里叶变换）fft_signal=fft(signal);% 构建频率轴：从 -fs/2 到 fs/2，共 num_samples 个点frequency_axis=linspace(-sampling_rate/2,sampling_rate/2,num_samples);% 绘制频谱% fftshift 将零频分量移到数组中心% 除以 num_samples 是为了归一化，使得谱线幅度等于信号的实际幅度（双边谱）plot(frequency_axis,fftshift(abs(fft_signal))/num_samples);xlabel('Frequency (Hz)');ylabel('Magnitude');title('采样信号的归一化双边频谱 (Normalized Bilateral Spectrum)');grid on;

频率分辨率推导（极其重要）

频率分辨率( \Delta f ) 定义为离散频谱中相邻两个频率点之间的间隔。其计算公式为：
[
\Delta f = \frac{f_s}{N}
]
其中 ( N ) 是 FFT 的点数（这里等于num_samples）。
代入数值：
[
\Delta f = \frac{5000 \text{ Hz}}{5000} = 1 \text{ Hz}
]
物理含义：在这个频谱图中，每个“柱子”（频点）代表一个 1 Hz 宽的频带。这意味着系统能够区分相差 1 Hz 以上的两个频率成分。

频谱结果解读

• 在 ( f = +1000 \text{ Hz} ) 和 ( f = -1000 \text{ Hz} ) 处，各出现一根清晰的谱线。
• 谱线幅度为0.5（因为 ( \sin(2\pi f t) ) 的双边谱幅度是单边幅度 1 的一半）。
• 其余频率处幅度几乎为 0（忽略数值计算中的微小浮点误差），证明信号非常纯净。

✂️ 第三部分：非线性失真——硬削波（对应代码 Q3）

任务描述

创建一个幅度为 1.0、频率为2.25 kHz的正弦信号（仍然使用 5 kHz 采样率）。
将这个信号“硬削波”至最大绝对值0.3：

• 所有大于 0.3 的样本值设为 0.3；
• 所有小于 -0.3 的样本值设为 -0.3。

绘制削波后的时域信号，并放大前 1 ms 以观察削波细节。

代码实现（逐行解析）

frequency=2250;% 2.25 kHzsampling_rate=5000;% 保持不变duration=1;max_amplitude=0.3;% 削波阈值num_samples=duration*sampling_rate;time=(0:num_samples-1)/sampling_rate;% 生成原始正弦波（幅度 1.0）signal=sin(2*pi*frequency*time);% === 硬削波操作（核心） ===% 利用 MATLAB 逻辑索引，将超出阈值的部分直接截断signal(signal>max_amplitude)=max_amplitude;signal(signal<-max_amplitude)=-max_amplitude;% 绘图：完整波形 + 局部放大figure(3);subplot(2,1,1);plot(time,signal);xlabel('Time (s)');ylabel('Amplitude');title('2.25kHz 正弦波，硬削波至 ±0.3 (Clipped)');grid on;subplot(2,1,2);plot(time,signal);axis([0,0.001,-0.4,0.4]);% 聚焦前 1 ms（约 2.25 个周期）xlabel('Time (s)');ylabel('Amplitude');title('局部放大 (Zoomed) - 可见削平顶部');grid on;

时域现象解读

• 原始 2.25 kHz 正弦波的峰值本来应达到 ±1.0，但现在被“削”成了平坦的直线 ±0.3。
• 这种“削平”动作在时域上产生了尖锐的转折角。转折角越尖锐，对应频域的高频成分越丰富。
• 关键前提检查：基波 2.25 kHz 本身满足 ( 2 \times 2250 = 4500 < 5000 )，所以基波自身不发生混叠。但削波产生的谐波会带来麻烦。

🎯 第四部分：削波信号的频谱与混叠定量分析（对应代码 Q4）—— 本文核心

任务描述

首先在纸上画出连续时域削波信号的预期频谱，再画出采样后预期看到的频谱。然后对削波信号做 FFT，对比实际与预期。

理论预期（连续时域）

在连续时间下，对称硬削波一个正弦波，其傅里叶级数只包含奇次谐波：

• 基波：( f_1 = 2.25 \text{ kHz} )
• 3次谐波：( f_3 = 6.75 \text{ kHz} )
• 5次谐波：( f_5 = 11.25 \text{ kHz} )
• 7次谐波：( f_7 = 15.75 \text{ kHz} )
• 9次谐波：( f_9 = 20.25 \text{ kHz} )
• ……（无穷多项）

谐波幅度随阶数增加而递减，但削波越严重（阈值越低），高次谐波衰减越慢。

采样后的混叠定量计算（关键）

当这个连续信号被 ( f_s = 5 \text{ kHz} ) 采样后，频谱会以 5 kHz 为周期进行无限复制。我们只关心主值区间 ( [-2.5 \text{ kHz}, +2.5 \text{ kHz}] )（奈奎斯特区间）。
应用混叠公式 ( f_{alias} = | f_{harm} - k \cdot f_s | )，其中 ( k ) 为整数，使得 ( f_{alias} \le 2.5 \text{ kHz} )：

关键结论：所有高次谐波无一例外地“折叠”进了低频区域。它们与基波（2.25 kHz）混合在一起，在频谱上制造出一系列“虚假”的峰值。

代码实现

figure(4);fft_signal=abs(fft(signal));% 取模frequency_axis=linspace(-sampling_rate/2,sampling_rate/2,num_samples);plot(frequency_axis,fftshift(fft_signal)/num_samples);xlabel('Frequency (Hz)');ylabel('Magnitude');title('削波信号 (2.25kHz) 的实际 FFT 频谱');xlim([-3000,3000]);% 聚焦在基波与混叠成分密集区grid on;

实际频谱解读（对比预期）

运行上述代码后，您将看到：

1. 最大的峰值出现在 ( \pm 2.25 \text{ kHz} )（基波）。
2. 显著的峰值出现在 ( \pm 1.75 \text{ kHz} )（这正是 3 次谐波混叠进来的位置）。
3. 较小的峰值出现在 ( \pm 1.25 \text{ kHz} )（5 次谐波混叠）。
4. 更微弱的峰值出现在 ( \pm 0.75 \text{ kHz} ) 和 ( \pm 0.25 \text{ kHz} )。
5. 实际频谱与我们的理论计算表完美吻合。这直观地证明了：采样率不足时，非线性失真产生的谐波会严重污染带内信号。

🚀 第五部分：过采样的神奇效果——混叠的“解药”（对应代码 Q6）

任务描述

将采样率从 5 kHz 大幅提升至100 kHz，再次对 2.25 kHz 正弦波进行削波（±0.3），观察时域波形并思考频谱变化。

代码实现

amplitude=1.0;frequency=2.25e3;sampling_rate=100e3;% 提高到 100 kHz！duration=1;t=linspace(0,duration,duration*sampling_rate);signal=amplitude*sin(2*pi*frequency*t);signal(signal>0.3)=0.3;signal(signal<-0.3)=-0.3;figure(6);plot(t(1:500),signal(1:500));% 只取前 5ms 展示（因为点太多）axis([0,0.005,-0.4,0.4]);xlabel('Time (s)');ylabel('Amplitude');title('过采样 (100kHz) 削波信号，时域分辨率大幅提升');

为什么过采样能拯救频谱？

• 此时奈奎斯特频率为 ( f_s/2 = 50 \text{ kHz} )。
• 我们计算一下各次谐波是否还在带内（远低于 50 kHz）：
  • 3次：6.75 kHz
  • 5次：11.25 kHz
  • 7次：15.75 kHz
  • 9次：20.25 kHz
  • 11次：24.75 kHz
  • 13次：29.25 kHz
  • 15次：33.75 kHz
  • 17次：38.25 kHz
  • 19次：42.75 kHz
  • 21次：47.25 kHz（仍小于 50 kHz！）
• 一直到 21 次谐波都不会发生混叠，因为全部低于奈奎斯特频率。
• 频谱表现为在 2.25k、6.75k、11.25k……处的一系列孤立谱线，没有折叠成分。
• 过采样为后续的数字低通滤波提供了巨大的“保护间隔”，使得我们可以轻松滤除高次谐波，从而恢复出纯净的基波。

🔲 第六部分：辅助验证——方波频谱（对应代码 Q7）

任务描述

生成一个 1.1 kHz 的方波（采样率 5 kHz），观察其频谱。方波天生包含奇次谐波，可以用来验证混叠的普遍性。

代码要点

f=1100;% 1.1 kHzsampling_rate=5e3;t=linspace(0,1,1*sampling_rate);square_wave=square(2*pi*f*t);% MATLAB 自带方波函数N=length(square_wave);freq_axis=(-N/2:N/2-1)*(sampling_rate/N);spectrum=abs(fftshift(fft(square_wave)));figure(7);plot(freq_axis,spectrum);xlabel('Frequency (Hz)');ylabel('Magnitude');title('1.1kHz 方波的频谱 (可见混叠)');xlim([-3000,3000]);

观察结果：方波的 3 次谐波（3.3 kHz）、5 次谐波（5.5 kHz）等，由于超过 2.5 kHz 的奈奎斯特频率，会折叠回低频区（如 0.7 kHz、0.5 kHz 等），进一步验证了混叠是采样系统的普遍现象，而非削波特例。

🔬 第七部分：削波阈值的影响（对应代码 Q8）

任务描述

将削波阈值改为更严苛的±0.1，采样率设为10 kHz，观察频谱变化。

现象分析

• 阈值越低（0.1），削波越“狠”，高次谐波的能量衰减越慢，频谱中的谐波幅度相对更大。
• 虽然采样率提高到 10 kHz（奈奎斯特 5 kHz），但 3 次谐波（6.75 kHz）仍然会混叠到 ( 6.75 - 5.00 = 1.75 \text{ kHz} )。
• 这组对比实验强调了：削波程度和采样率共同决定了频谱污染的程度。即使提高采样率，只要谐波超出新的奈奎斯特频率，混叠就依然存在。

💎 总结与工程启示

这份超详细的指南从数学公式、代码实现到频谱混叠的定量计算，进行了全方位的拆解。希望它能彻底解决您在采样、削波和混叠方面的所有疑惑！
This ultra-detailed guide breaks down everything from mathematical formulas and code implementation to quantitative aliasing calculations. We hope it thoroughly resolves all your questions regarding sampling, clipping, and aliasing!😊