简介《数字信号处理MATLAB仿真》是一份面向电子信息、通信工程等专业学生和工程技术人员的PDF学习文档聚焦如何借助MATLAB完成数字信号处理中连续信号与离散信号的建模、运算和仿真分析。文档以实验形式组织实验目的明确设备要求简单只需一台装有MATLAB的计算机即可上手。内容从单位冲击信号、单位阶跃函数、斜坡函数、实指数函数、正弦函数等连续信号讲起逐步覆盖单位冲激序列、任意序列、单位阶跃序列、斜坡序列、正弦序列、实指数序列、复指数序列和随机序列等离散信号并延伸到信号延迟、卷积运算、DFT/FFT频域分析以及IIR与FIR数字滤波器设计等核心内容。针对每条知识点文档都提供了可直接运行的MATLAB程序示例例如用stem函数绘制冲激序列、用stepfun生成阶跃函数、用rand/randn产生均匀分布与高斯随机序列并附有输出图形与简要说明便于读者边看边练从时域和频域两个角度理解信号特性。资源包为单个PDF文件大小约1.15MB内容紧凑实用适合作为实验指导书、期末复习资料或自学入门手册。目前已有452人学习浏览是快速掌握MATLAB数字信号处理仿真的可靠参考。1. 数字信号处理的 MATLAB 仿真从实验手册到工程基线做完互联网语音传输VoIP里的回声消除模块后回头看实验课上那套 stepfun、stem、fft、butter 的操作几乎就是数字信号处理全部核心操作的骨架。这份实验指导的特别之处在于它把连续信号建模、离散序列生成、序列运算、DFT/FFT、IIR/FIR 滤波器设计串成了一条完整链路每个环节都给出了可直接运行的 MATLAB 代码。对刚接触 DSP 的本科生它是入门路径对需要在 MATLAB 里快速验证算法的工程师它是拿来即用的参数手册。下面按信号建模、序列运算、频域分析、滤波器设计四个层次逐段拆解。2. MATLAB 信号建模stepfun、stem 与连续/离散序列生成2.1 连续信号的近似与 MATLAB 表达从单位冲击到正弦函数离散化是连续信号仿真的前提。以单位冲击信号为例理论定义要求脉冲宽度趋近于零、幅度趋近于无穷计算机无法直接表示所以工程上的做法是用宽度为 1/A、幅度为 A 的窄矩形脉冲来近似。文档例 1.1 的实现值得细读clear all; t1 -0.5:0.001:0; A 50; A1 1/A; n1 length(t1); u1 zeros(1, n1); t2 0:0.001:A1; t0 0; u2 A * stepfun(t2, t0); t3 A1:0.001:1; n3 length(t3); u3 zeros(1, n3); t [t1 t2 t3]; u [u1 u2 u3]; plot(t, u); axis([-0.5 1 0 A2]);代码将时间轴分成三段左半段 t1 幅度为 0中间 t2 通过 stepfun 在 t00 处产生阶跃并放大 A 倍右半段 t3 归零。stepfun(t2, t0) 返回与 t2 等长的 0/1 向量元素在 t t0 时为 1否则为 0这正是阶跃函数的离散近似。u2 直接乘以 A 后得到宽度 0.02、高度 50 的窄脉冲由 plot 绘制。axis 范围留出顶部余量避免脉冲被坐标轴裁切。当 A 越大脉冲越窄越接近理论上的单位冲击 δ(t)。在单位冲击基础上其他连续信号就容易理解了。单位阶跃直接调用 stepfun(t, t0)斜坡函数 g(t)3(t-1) 用阶跃序列逐点乘以 (t - t0) 得到实指数 f(t)3e^{0.5t} 直接用 Aexp(at) 对向量逐元素求幂。正弦信号的实现需要注意相位参数clear all; t -0.5:0.001:1; A 3; f 5; fai 1; u A * sin(2*pi*f*t fai); plot(t, u); axis([-0.5 1 -3.2 3.2]);这里 sin 接收的是弧度制参数2pif*t 把频率 f 从赫兹转换为角频率fai 是初始相位。文档例子里 f5 对应周期 0.2 秒在 -0.5 到 1 的区间内可以观察到约 7.5 个完整周期axis 的 y 轴范围略大于幅度 3避免波形触顶。需要指出的是MATLAB 向量化要求 exp、sin 这类函数直接作用在向量上得到等长结果这也是文档所有连续信号实现共同遵循的规则。下面给出连续信号与对应 MATLAB 表达式的对照表实际使用时可以直接把表格里的表达式复制到脚本中调整参数信号类型理论表达式MATLAB 实现关键参数单位冲击δ(t)stepfun 构造窄脉冲A 决定脉冲宽度与高度单位阶跃u(t)stepfun(t, t0)t0 为跳变时刻斜坡g(t)B(t-t0)Bu.(t-t0)B 为斜率实指数f(t)Ae^{at}Aexp(at)a 为增长/衰减因子正弦f(t)Acos(2πftφ)Asin(2piftfai)f 为频率fai 为相位提示stepfun 等价于离散化的 heaviside它返回的是与时间向量等长的 0/1 向量直接参与乘法运算时会被当作数值处理。2.2 离散序列的生成stem 绘图与常用序列离散序列的绘制固定使用 stem它以样本索引为横轴、样本值为纵轴画出垂直杆。文档从单位冲激序列开始xzeros(1,N); x(1)1; 就构造了 64 点的单位冲激。类似地单位阶跃序列用 ones(1,N)斜坡序列在阶跃基础上逐点乘以 (i-k)正弦序列则把频率用样本索引归一化xAsin(2pif(xn/N)fai)。均匀分布与高斯随机序列的对比是实验中容易忽略的点rand(1,N) 产生 [0,1] 上的均匀分布randn(1,N) 产生均值为 0、方差为 1 的高斯白噪声。两者的差异直接决定仿真信号的统计特性。下面的代码一次生成四种常用离散序列clear all; N 32; % 实指数序列 A 3; a 0.7; xn 0:N-1; x_exp A * a.^xn; % 复指数序列j 为虚数单位 w 314; x_cplx A * exp((a j*w) * xn); % 均匀分布与高斯随机序列 x_rand rand(1, N); x_randn randn(1, N); % 分别绘图 subplot(2,2,1); stem(xn, x_exp); title(实指数序列); subplot(2,2,2); stem(xn, abs(x_cplx)); title(复指数序列幅度); subplot(2,2,3); stem(xn, x_rand); title(均匀随机序列); subplot(2,2,4); stem(xn, x_randn); title(高斯随机序列);实指数序列中的 a.^xn 使用了逐元素幂运算符a 是标量、xn 是向量逐元素运算得到 x(n)A·a^n。复指数序列 exp((aj*w)*xn) 中 j 是 MATLAB 内置虚数单位w314 约等于 100π 弧度对应 50 Hz 离散角频率。stem 对复数输入取实部绘图所以这里用 abs 取幅度。rand 与 randn 每调用一次产生一组新序列若需要可复现的仿真数据可以在脚本开头加入 rng(固定种子)。3. 离散序列运算的 MATLAB 实现延迟、对齐补零与 sum/prod3.1 信号延迟与翻转时间轴的移动离散信号延迟的定义很简单y(n)x(n-k) 表示序列右移 k 个抽样周期即在原序列前补 k 个零。文档例 3.1 用 zeros(1,k) 做前补零再拼接原序列实现延迟代码如下clear all; N 32; w 100; k 3; xn 0:N-1; x2 sin(100 * xn); % 原始正弦序列 x [zeros(1, k), x2]; % 前补 k 个零实现右移 stem(0:Nk-1, x);前补零的物理意义是让序列整体在时间轴上向右移动 k 个样本补零不改变原序列幅度只改变起始位置。与之相对信号翻转 y(n)x(-n) 用 fliplr 实现它对行向量做左右镜像。这里有一个容易踩的坑fliplr 只翻转数据顺序不改变横轴索引所以绘图时如果沿用原 n看起来是值翻转而轴不变更像镜像而不是时间反转。正确做法是为翻转后的序列单独构造索引轴 -n。3.2 相加与相乘长度对齐是前提两个序列相加前必须满足两个条件长度相等、样本位置对应。文档例 3.2 中 x1 只有 4 个点而 x2 有 8 个点直接相加会报维度错误正确做法是用 zeros 补零让两者对齐clear all; n1 0:3; x1 [2 0.5 0.9 1]; n2 0:7; x2 [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7]; n 0:7; x1 [x1 zeros(1, 8 - length(n1))]; % x1 右侧补零对齐 x2 的结束位置 x2 [zeros(1, 8 - length(n2)), x2]; % x2 长度已是 8左补零长度为 0 x x1 x2;补零方向决定了样本位置的对应关系x1 右补零意味着它占据时间轴的 0~3 位置x2 不变占据 0~7二者在 0~3 区间相加、4~7 区间只有 x2 的值。相乘与相加完全一致先把两个序列对齐补零再用 .* 做逐样本点乘。文档例 3.3 把加法改成乘法后4~7 区间结果全部为 0因为 x1 在该区间补的零与 x2 相乘仍为零。注意补零对齐是序列运算里最容易出错的环节。如果两个序列起始索引不同比如一个从 n2 开始需要分别在左侧补与偏移量相等的零才能保证样本在时间轴上真正对齐。3.3 信号和与信号积sum 与 prod序列和 ysum(x) 对所有样本求和序列积 yprod(x) 对所有样本连乘。两个函数都沿第一个非单一维度操作对行向量直接返回标量。序列积有一个特性值得注意只要序列中有一个零样本连乘结果必为零因此 prod 通常用于检验序列是否包含零值或者计算增益链的总衰减。下面把离散序列的常用运算整理成对照表运算数学表达MATLAB 实现注意事项延迟y(n)x(n-k)[zeros(1,k) x]k 为右移样本数翻转y(n)x(-n)fliplr(x)横轴需单独处理相加x1(n)x2(n)x1x2先补零对齐相乘x1(n)*x2(n)x1.*x2必须用点乘信号和Σ x(n)sum(x)返回标量信号积Π x(n)prod(x)含零则为 0提示对于不等长序列也可以用 interp1 或 resample 先做重采样再用算术运算但补零是最轻量、最容易验证正确性的方案适合教学与快速原型验证。4. 从 DFT 到 FFT 的 MATLAB 实现旋转因子矩阵与 fft 参数4.1 直接 DFT 矩阵计算理解旋转因子 WN文档例 4 用矩阵方式实现离散傅里叶变换代码很短但信息密度很高N 32; n 0:N-1; xn cos(pi*n/6); % 输入序列 k 0:N-1; WN exp(-j*2*pi/N); % 旋转因子 nk n * k; % 构造 n×k 矩阵 WNnk WN.^nk; % 旋转因子矩阵 Xk xn * WNnk; % 矩阵乘法完成 DFTDFT 的定义是 X(k)Σ x(n)·W_N^{nk}其中 W_N e^{-j2π/N}。nk 生成了一个 N×N 的指数矩阵每一行对应一个时间样本 n每一列对应一个频点 kWN.^nk 逐元素求幂后xnWNnk 就是向量与矩阵的乘法本质上是一次性完成 N 个频点的累加求和。这种写法在教学中非常直观它把 DFT 的累加过程展开成了矩阵运算便于对照公式检查实现是否有误。但矩阵方式的时间复杂度是 O(N²)。当 N32 时看不出问题当 N 达到数千甚至上万时矩阵会占据大量内存。FFT 通过蝶形运算把复杂度降到 O(NlogN)这就是为什么工程实现里几乎总是用 fft 而不是手写 DFT。4.2 fft 的函数语义与频谱绘制MATLAB 的 fft 函数是机器语言实现的执行速度远快于用户脚本。常用格式有两种yfft(x) 使用与输入长度相等的变换点数yfft(x,N) 强制使用 N 点变换如果 x 长度小于 N 则自动补零大于 N 则截断。补零不增加频率分辨率只是对频谱做插值这个区别经常被混淆处理方式含义频率分辨率fft(x)变换点数 数据长度Fs/len(x)fft(x, N) 且 N len(x)尾部补零仍为 Fs/len(x)谱线更密fft(x, N) 且 N len(x)直接截断可能丢失频谱细节绘制频谱的标准流程是对变换结果取幅度 abs(Xk)以频率轴为横轴用 stem 或 plot 绘图。文档例 4 中直接绘制了 abs(Xk) 的离散幅度谱。实际工程中还需要做两件事一是用 fftshift 把零频移到中心二是对幅度做 1/N 归一化否则幅度值不是信号的真实幅度。下面给出一个带归一化的频谱绘制模板Fs 1000; N 64; n 0:N-1; xn cos(2*pi*100*n/Fs); % 100 Hz 余弦信号 Xk fft(xn, N); mag abs(Xk) / N; % 幅度归一化 f (0:N-1) * Fs / N; % 频率轴 stem(f, mag); xlabel(频率 (Hz)); ylabel(幅度);这里 Fs 是采样率频率轴 f(0:N-1)*Fs/N 每个频点的间隔是 Fs/N。归一化除以 N 后单频余弦在峰值频点处的幅度约为 0.5对应 cos 展开成两个指数分量各占一半幅度。若信号包含直流分量则 0 Hz 频点归一化后幅度即直流值。需要留意的是 fft 输出的前 N/21 个点对应 0 到 Fs/2 的正频率后一半是负频率镜像对实信号而言通常只画前一半即可。提示文档中的旋转因子矩阵写法用到了 exp(-j2pi/N)这里的负号对应 DFT 的正变换。如果符号取反得到的是逆变换画出的频谱会沿频率轴镜像。5. IIR 数字滤波器设计butter 与 cheby1/cheby2 的参数化实践5.1 巴特沃斯带通滤波器的参数语义文档例 5.1 用一行代码设计了 10 阶带通巴特沃斯滤波器N 10; Wn [100 200] / 500; % 归一化通带边界 [b, a] butter(N, Wn, bandpass); freqz(b, a, 128, 1000); % 128 点频响采样率 1000 Hz [y, t] impz(b, a, 101); % 冲激响应 stem(t, y);butter 的第二个参数 Wn 是归一化频率范围在 0 到 1 之间1 对应奈奎斯特频率采样率一半。这里的 500 就是 1000 Hz 采样率对应的奈奎斯特频率100 Hz 和 200 Hz 分别除以 500 得到 0.2 和 0.4。freqz 的第三个参数 128 表示计算 128 个频点第四个参数 1000 是绘图用的采样率只影响频率轴标注不影响滤波器系数本身。巴特沃斯滤波器的特点是通带内最大平坦但过渡带较宽。10 阶带通的带外衰减斜率约为每倍频程 60 dB10 阶 × 6 dB/oct这个指标是否够用取决于实际应用场景。如果过渡带要求更陡就需要引入切比雪夫滤波器。5.2 切比雪夫Ⅰ型与Ⅱ型用纹波换过渡带切比雪夫滤波器的核心思想是允许通带或阻带存在等纹波换取更陡的过渡带。文档例 5.2 给出了一个完整的低通设计流程Wp 100; Rp 3; Ws 200; Rs 30; Fs 1000; [N, Wn] cheb1ord(Wp/(Fs/2), Ws/(Fs/2), Rp, Rs); [b, a] cheby1(N, Rp, Wn); freqz(b, a, 512, 1000);cheb1ord 的前两个参数分别是通带边缘 Wp 和阻带边缘 Ws 的归一化频率后两个是通带最大纹波 RpdB和阻带最小衰减 RsdB。函数返回满足指标的最小阶数 N 和对应的 3 dB 截止频率 Wncheby1 再用阶数、纹波和 Wn 计算滤波器系数。切比雪夫Ⅰ型在通带有纹波、阻带单调切比雪夫Ⅱ型则相反通带单调、阻带有纹波。文档例 5.3 的带通切比雪夫Ⅱ型设计使用 cheb2ord参数结构与 cheb1ord 完全一致只需要把通带边界写成两元素向量 Wp[100,250]阻带边界 Ws[50,300]Wp [100 250]; Ws [50 300]; Rp 3; Rs 30; Fs 1000; [N, Wn] cheb2ord(Wp/(Fs/2), Ws/(Fs/2), Rp, Rs); [b, a] cheby2(N, Rs, Wn); freqz(b, a, 512, 1000);注意 cheby2 的第二个参数是阻带衰减 Rs 而不是通带纹波 Rp这是与 cheby1 最容易混淆的地方。文档原代码此处写的是 Rp实际设计时应传 Rs否则阻带衰减指标无法满足。把三种 IIR 设计方式整理成对照表方便按需求选择设计方法函数通带特性适用场景巴特沃斯butter最大平坦对通带纹波敏感过渡带要求不高切比雪夫Ⅰcheb1ord cheby1等纹波需要更陡过渡带能容忍通带纹波切比雪夫Ⅱcheb2ord cheby2阻带等纹波对阻带衰减要求严格设计完成后用 freqz 绘制的幅频曲线可以直接读出通带纹波和阻带衰减是否达标。如果阶数过高导致数值不稳定优先考虑改用椭圆滤波器 ellipord/ellip或者在设计后通过 zp2sos 把传递函数分解为二阶节级联。6. FIR 滤波器设计与频响验证窗函数、fir1 与 freqz 排错6.1 窗函数选择与 fir1 的参数语义FIR 设计的核心参数是阶数、截止频率、窗函数类型。MATLAB 提供多种窗函数生成函数boxcar(n) 矩形窗、triang(n) 三角窗、hanning(n) 汉宁窗、hamming(n) 海明窗、blackman(n) 布拉克曼窗、kaiser(n,beta) 凯泽窗。窗函数的选择直接决定阻带衰减矩形窗约 21 dB汉宁约 44 dB海明约 53 dB布拉克曼约 74 dB凯泽窗可通过 beta 参数连续调节。fir1 的调用格式为 fir1(n, Wn, ftype, Window)其中 n 是阶数Wn 是归一化截止频率标量对应低通/高通[W1 W2] 对应带通/带阻ftype 缺省为低通high 和 stop 分别对应高通与带阻。下面用布拉克曼窗设计一个带通滤波器Window blackman(16); % 窗长度 阶数 1 b fir1(15, [0.3 0.5], Window); % 通带 0.3π~0.5π freqz(b, 1, 512);6.2 用 freqz 与 impz 完成设计验证验证滤波器是否符合指标我一般会做三步先用 freqz 看幅相特性确认通带边界是否落在设计值、阻带衰减是否达到预期再用 impz 看冲激响应的长度和对称性fir1 设计的线性相位滤波器冲激响应应对称最后用 filter 把真实信号送入滤波器对比输入输出频谱。文档例 6.3 中的 chirp 实验正是第三种验证方式的完整示范——加载 chirp.mat 信号用 chebwin(35,30) 设计 34 阶高通滤波后用 pburg 画滤波前后频谱直观看到低频分量被抑制。注意fir1 中窗函数长度必须是阶数加 1比如阶数 34 要用 chebwin(35,30)写错长度 MATLAB 会报错或静默截断。设计完成后检查 freqz 曲线在通带边缘是否出现明显的过冲如果有优先增大阶数或改用凯泽窗。本文还有配套的精品资源点击获取