MATLAB数字信号处理实战：从源码包到滤波器设计

简介万永革《数字信号处理的MATLAB实现》高清扫描版是一本以大量实例和可运行程序贯穿始终的数字信号处理实践教材主要面向数理基础相对薄弱、重视动手能力的读者也适合本科生、研究生与工程技术人员作为案头参考。压缩包内共476个文件整体约54.24MB以364个.m脚本为主体附带.dat/.mat数据文件、.txt/.doc说明文档、.tex源文件及高清PDF文件类型覆盖程序、数据、文档三类配合完整目录结构方便按章节定位和调用。书内全部例题的MATLAB实现及其配套数据均已收录理解算法思路后只需修改参数即可复用至实际工程场景能有效降低上手门槛。目前已有530人学习下载适合希望在MATLAB环境下系统掌握数字信号处理应用方法的读者。1. 数字信号处理的MATLAB实现先跑通源码包再谈理解如果你想找的是数字信号处理的MATLAB实现而不是又一篇讲傅里叶变换定义的文章这份资源值得下。它在做的事是把离散信号生成、差分方程求解、FFT频谱分析、FIR/IIR滤波器设计、功率谱估计拆成一段段能直接运行的脚本让理论公式第一次变成可观察的曲线。它适合三类人期末复习时想验证概念的本科生、课设前需要快速搭框架的同学、以及做信号处理Demo时不想重写基础函数的从业者。核心价值不是提供答案而是告诉你MATLAB里那些函数真正要什么参数。我的结论比较直接跑通它比把书抄一遍更接近数字信号处理。2. 拆资源包按“时域→频域→滤波器”链路读脚本资源包里的脚本数量看上去不少但真正需要精读的不会超过二十个。按信号处理链路拆开你会发现它无非三块时域系统、频域变换、滤波器设计。时域脚本负责让你把离散信号和系统响应“看得见”频域脚本负责把谐波成分找出来滤波器脚本负责把不要的频率压下去。下面逐个拆每一步都给出能直接复现的写法。打开资源包后我建议不要从第一行顺序读到结尾而是按“时域→频域→滤波器”的顺序挑脚本看。先读时域基础因为filter的用法会贯穿整个滤波器章节再读FFT因为频谱图是验证所有滤波效果的依据最后读滤波器这时候你已经知道横轴怎么换算、幅值怎么取单边不会被图上奇怪的刻度带偏。2.1 时域篇信号生成与差分方程脚本怎么读时域基础部分的脚本一般负责生成单位脉冲、阶跃、指数衰减序列再用filter求解差分方程。很多同学一开始只盯着结果图我却建议先看y filter(b, a, x)这一行因为它是后面所有滤波例程的地基。filter函数并不神秘它按系统函数的分子分母系数逐点递归计算输出长度默认和输入一样。下面这段简化脚本就是这类例程最常见的骨架。n 0:15; % 离散时间轴从0到15 x (0.8).^n; % 指数衰减输入序列 b 1; % 系统函数分子系数 a [1, -0.8]; % 分母系数y[n] - 0.8*y[n-1] x[n] y filter(b, a, x); % 求解差分方程得到响应序列 subplot(2,1,1); stem(n, x, filled); title(输入 x[n]); subplot(2,1,2); stem(n, y, filled); title(系统响应 y[n]);逻辑说明filter不是狭义“滤波”函数而是按差分方程递归求解的通用工具。b和a来自传递函数H(z)B(z)/A(z)这里b1表示输入只直接作用于当前输出a[1,-0.8]表示当前输出还依赖上一时刻输出乘以0.8。参数说明subplot把两幅图纵向排开stem画离散序列如果这里改成plot会把点连成线掩盖离散信号的跳变特征。读这类脚本时我一般会把系数改成[1,-0.5]、[1,0.8]重新运行看波形变化。改分母系数是在验证系统是否稳定极点接近单位圆时序列衰减变慢极点一旦超出单位圆输出会指数发散。资源包里的例程通常只给稳定参数但课设答辩时老师最爱问“系数变了会怎样”所以动手改一改比背结论有用。2.2 频域篇FFT 脚本的取数与横轴换算频谱分析脚本是资源包里最能拉开差距的部分。fft函数本身只有一行但后续的取数方式直接决定图能不能看。常见错误是把fft结果直接plot(abs(X))横轴停留在“点序号”没法标Hz。另一个常见问题是忘了实信号频谱的镜像关系导致正弦波幅值显示成真实值的一半。下面这段脚本是频域例程的标准起手式。N 512; % FFT点数决定频率分辨率 fs 1000; % 采样率单位Hz t (0:N-1)/fs; % 时间轴单位秒 x sin(2*pi*50*t) 0.3*sin(2*pi*380*t); X fft(x); % 复数频谱长度等于N f (0:N-1)*fs/N; % 把点序号换算成物理频率 plot(f, abs(X), LineWidth, 0.8); xlabel(频率 (Hz)); ylabel(|X(f)|);这段代码里的频率换算必须记住fft返回的复数数组每个下标对应的数字频率是(idx-1)*fs/N所以频率轴写成(0:N-1)*fs/N。画实信号频谱时正负频率对称一般只取前N/21个点做单边谱如果直接显示全部50Hz正弦的峰值是0.5而不是1。修整方式很简单X(1)是直流X(N/21)是奈奎斯特频率这两个点不用乘其余部分乘2。如果要观察负频率比如调制信号可以在fft之后加Xfftshift(fft(x))频率轴改成f(0:N-1)*fs/N - fs/2。fftshift把零频移到中间正负频率对称结构直接可见。这个函数在资源包的对称性示例里不一定单独出现但读“信号调制”类脚本时总会遇到提前掌握能少卡壳一次。关于补零也要多说两句fft(x,2048)是把序列补零到2048点做变换曲线更平滑但两个频率靠得近时不一定会分开因为真实频率分辨率受原始信号长度限制。资源包如果用了补零你需要在报告里写清“补零提高的是插值密度不是真实分辨率”。2.3 滤波器篇设计函数与filtfilt的边界滤波器例程是资源包的结尾部分也是课设最常直接取用的部分。设计阶段一般调用fir1或butter使用阶段常见filter讲究一点会用filtfilt。设计时最需要算准的是归一化截止频率Wn它是物理截止频率除以奈奎斯特频率之后得到的无量纲数。下面这段脚本代表FIR低通的标准流程。fs 1000; % 采样率 fc 100; % 截止频率单位Hz Wn fc/(fs/2); % 归一化截止频率 b fir1(32, Wn, low); % 32阶FIR低通 x randn(1, 1000); % 白噪声信号 y filter(b, 1, x); % 普通滤波这里Wn0.2fir1内部在设计数字频率0.2π处的低通滤波器。阶数32表示滤波器有33个抽头阶数越高过渡带越窄但输出延迟和计算量同步上升。使用filter时注意起始瞬态因为滤波器初始状态为零输出开头几十个点受暂态影响展示波形时最好截掉前一段。如果资源包里使用了filtfilt它的作用是先正向滤波再反向滤波从而抵消相位偏移让输出看起来没有延迟。代价是非因果特性不能用在实际实时系统中。做课设演示时我会把filter和filtfilt两种结果都跑一遍报告里注明差异这样的细节比单纯贴代码更容易得分。2.4 随机信号与功率谱pwelch 为什么不能乱用窗功率谱估计脚本在资源包最后一段经常出现用的通常是pwelch。它把长序列分成多段每段加窗、做FFT、再对幅值平方做平均。这个方法比直接取FFT平滑因为平均降低了方差但代价是频率分辨率下降。常见的神逻辑是把pwelch输出和fft输出画在同一张图里然后比较幅值大小实际上两者量纲都不同一个是功率谱密度一个是线性谱。下面是pwelch的最小调用。fs 1000; N 8192; rng(0); x randn(N,1) 0.5*sin(2*pi*100*(0:N-1)/fs); [pxx, f] pwelch(x, hamming(512), 256, 1024, fs); plot(f, 10*log10(pxx)); xlabel(频率 (Hz)); ylabel(功率谱密度 (dB/Hz));参数很好理解512是窗长窗长越长频率分辨率越高但方差越大256是相邻窗的重叠点数重叠越大段间相关性越高1024是FFT点数fs是采样率。改窗长和重叠能看到曲线平滑度的明显变化。窗函数选择也有讲究矩形窗频率分辨率高但旁瓣泄漏大汉明窗旁瓣低但主瓣稍宽。资源包里常用汉明窗因为语音和噪声场景更看重旁瓣抑制。如果你把窗长从512改成128谱曲线会剧烈抖动这不是代码坏了而是平均段数增加导致的方差减小和分辨率下降同时发生。先多跑几个参数再下结论。3. 跑起来环境配置、参数设定与第一个可复现实验源码包是否能跑七成取决于运行环境三成才是脚本本身。这里先把“能不能跑”的检查做完再谈怎么把脚本改成自己的实验函数。很多第一次下源码包的同学会直接双击脚本文件弹出一堆红色报错后就关掉了其实问题往往出在工具箱缺失和路径冲突两件事上。3.1 先查工具箱与路径别急着双击拿到源码包后第一件事不是读代码而是检查运行环境。许多脚本会调用Signal Processing Toolbox里的函数缺了它会直接报“未定义函数”。在命令窗口执行以下几行5秒就能判断环境是否合格。ver(signal) % 查看信号处理工具箱版本 which fir1 % 找fir1所在路径 license(test,Signal_Toolbox) % 返回1表示有授权如果which fir1返回的不是工具箱内置路径而是某个以当前文件夹开头的路径说明资源包里有一个同名文件把系统函数盖住了。这个问题比工具箱缺失更隐蔽代码表面上能跑但结果可能不是官方实现。在MATLAB课程设计中最常见的人为冲突是自写fft函数一旦出现频谱图会变得很怪。版本兼容性也需要看一眼。老一些的脚本用wavread读音频新版MATLAB推荐audioread。如果运行提示函数已更名一般不影响结果但我会顺手改成audioread避免以后换机器再报“考虑使用audioread”的警告。这个警告虽然不致命但对答辩演示有点干扰。3.2 建项目文件夹处理路径与函数名冲突资源包解压后最好放进一个纯英文路径例如D:\DSP_Lab不要带中文或空格。路径含中文时某些版本的MATLAB读取数据文件会出玄学问题比如load(data.mat)明明存在却报文件不存在。这属于环境问题不是资源包的代码问题。把文件夹加进搜索路径后再检查一遍函数冲突。cd(D:\DSP_Lab); addpath(signal_scripts); which fft -allwhich fft -all会列出所有名为fft的路径如果排第一的不是toolbox内置路径就要改名。自写的函数一律改成dsp_fft这种不冲突的命名。这一点在把脚本改造成函数时尤其重要。运行脚本前还要注意很多旧例程第一行是close all; clear; clc;clear会清空工作区变量。如果你的实验脚本先加载了音频又调用了这一行音频变量会消失后面再引用就会报错。我一般会把这类开头的clear去掉改成close all; clc;保留工作区方便调试。3.3 把脚本改造成函数参数外置是复用关键源码包的脚本参数通常写死在代码里比如fs1000、N512。课设中你需要反复换采样率、换点数、换截止频率直接改脚本效率低还容易改乱。我习惯拿到一段能跑的脚本后先把它包成一个函数把主要参数提到函数签名里。function [X, f] my_spectrum(x, fs, N) % 输入x 时域序列fs 采样率N FFT点数 % 输出X 频谱复数f 频率轴 if nargin 3 N length(x); % 默认使用信号长度 end validateattributes(x, {double}, {vector}); % 确保输入是double向量 X fft(x, N); f (0:N-1)*fs/N; end函数化有一个好处主脚本不用关心计算过程只需要调用它。比如fs8000; [X,f]my_spectrum(x, fs, 1024);然后专心画图。nargin的作用是允许省略N时自动使用信号长度这在信号长度不固定时很实用。validateattributes做参数守卫输入不是double向量时直接报错避免把单精度或字符串误传进来。参数N的选择也有讲究N越大频率分辨率越细但计算量变大N取2的幂时FFT最快。如果信号只有8000点你取N65536结果只是曲线更光滑并没有真正提高分辨率。这一点在课设报告里写出来是能体现理解深度的。4. 避坑手册跑这套DSP代码最容易翻车的五个位置下面五条按“现象、原因、解决”的顺序列出来都是我在帮别人看代码时反复遇到的同类问题。不一定每个资源包都有但只要你用MATLAB做DSP迟早会碰到其中之一。4.1 频谱图的横轴显示的是点数而不是频率现象plot(abs(fft(x)))后横轴是0、100、200这类整数标不上Hz。原因fft输出的下标对应数字频率索引必须经过换算才能变成物理频率。解决先算频率轴f(0:N-1)*fs/N再plot(f, abs(X))如果只是xlabel(频率/Hz)而不设置横轴数据图上的刻度不会变。另一种更省事的方式是直接用freqz或periodogram它们可以接受采样率参数自动标Hz。4.2 用plot画离散序列把脉冲画成了锯齿现象单位脉冲序列用plot画出来是一条从0到1再回0的斜线完全看不出离散样点。原因plot把相邻样点用直线连接生成的是连续折线对离散序列是不准确的表示。解决画离散序列用stem如果你要显示包络可以用plot(n,x,o-)这种点线方式但正式报告里仍建议stem。这个坑在资源包的时域例程里几乎不会出现因为你直接复制代码就能看到正确画法但自己动手时容易顺手写plot。4.3 filter输出开头一段明显异常现象滤波后的信号前几十个点和稳态幅度差别很大看起来像毛刺或者零值。原因filter是因果递归滤波初始状态为零开头需要若干样本把内部状态建立起来这段暂态响应不是目标稳态。解决展示稳态波形时截掉开头一段或者用filtfilt做零相位滤波。如果是实时数据处理不能使用filtfilt需要设计合理的初始状态比如用filtic函数根据过去数据计算状态初值。课设里如果没做这个处理报告上的滤波波形会很丑。4.4 单边频谱幅度比预期小一半现象一个幅度为1的正弦信号fft后在50Hz处峰值显示0.5。原因实信号的频谱在正负频率上对称分布plot(abs(fft(x)))显示的是双边谱能量被对半分掉了。解决绘制单边谱时取前N/21个点乘以2直流和奈奎斯特频率两个点不乘。注意这个“乘2”只适用于实正弦信号如果信号本身是复指数频谱只有单边不需要乘2。这个细节很多教材不会强调但答辩时很容易被追问。4.5 butter滤波阶数过高导致输出异常现象用了butter(20, 0.4, low)后输出出现明显毛刺甚至MATLAB提示“Matrix is close to singular”。原因IIR滤波器阶数过高会让极点靠近单位圆量化误差和数值舍入误差被放大滤波器接近不稳定状态。解决IIR阶数控制在8以内需要更陡峭的过渡带时改用FIR或使用等波纹设计法。这也是资源包里FIR例程比IIR多的原因FIR没有反馈项天然稳定课设验收不容易翻车。如果必须用IIR先调用freqz检查幅频响应再决定最终阶数。5. 从源码包到课设成品做一个语音加噪与滤波实验源码包拆完下一步是把例程拼成一个能验收的完整项目。这里用一个最常见的语音滤波场景来串读入语音叠加固定频率干扰设计低通滤波器证明有效。这个思路你换到ECG信号、振动信号、音频降噪上都能复用只是信号源和截止频率不同。5.1 方案拼装把四段例程组合成一个主脚本典型的源码包会把“读文件、生成噪声、设计滤波器、播放对比”分散在多个脚本里课设报告不需要这么多文件。我习惯新建一个main_filter_demo.m把要点按流程合到一起。[y0, fs] audioread(speech.wav); % 读取语音 y0 y0(:, 1); % 取单声道 t (0:length(y0)-1)./fs; % 列时间轴 noise 0.05 * sin(2*pi*1800*t); % 1.8k正弦干扰 x y0 noise; % 加噪信号 [b, a] butter(4, 1200/(fs/2), low); % 4阶IIR低通 yf filtfilt(b, a, x); % 零相位滤波这个脚本看起来短但每一行背后都有对应例程支撑。audioread是读取模块sin是信号生成模块butter和filtfilt是滤波器模块。组合时最需要改的参数是截止频率这里的1200Hz是假设语音主要能量在1.2kHz以下而1.8k单频是噪声。如果你的语音采样率是44100Hz1200/(fs/2)只有0.0544低通会非常窄可能把语音也削掉。更稳妥的做法是先运行频谱分析脚本看噪声到底在哪个频点再定截止频率。5.2 验证用频域图与信噪比证明滤波器有效滤波器有没有效不能只看“听起来安静了”这一句主观判断。课设验收需要数字指标和谱图。常见做法是算滤波前后的SNR。snr_in 10*log10( sum(y0.^2) / sum((x - y0).^2) ); snr_out 10*log10( sum(y0.^2) / sum((yf - y0).^2) ); fprintf(SNR in %.2f dB, out %.2f dB\n, snr_in, snr_out);这里的y0是干净参考信号x是加噪信号yf是滤波输出。为什么用y0做参考因为SNR的定义是“信号功率/噪声功率”在没有加噪前录到的是参考信号源。如果资源包里没有干净参考信号可以用一段静音段估计噪声功率或者用频谱对比代替。注意如果用的是filter而不是filtfilt输出存在延迟直接做差算SNR会产生虚假误差需要先把yf对齐到y0再计算。5.3 参数标定的权衡滤波器设计不是越大越好也不是越陡越好。给一个直观的对比表参数影响建议值FIR阶数阶数越高过渡带越窄延迟和计算量越大20到100之间根据fs调整IIR阶数阶数高到8以上开始有数值稳定风险4到8截止频率太高滤不掉噪声太低削掉有效信号先用FFT看噪声频点再定采样率决定奈奎斯特频率与归一化折算使用录音原始fs不要随意resample这个表在报告里很有用老师问“为什么选4阶”时可以直接答兼顾衰减斜率和稳定性。如果源码包里的脚本直接给的是Wn而非Hz写报告时把你的转换公式Wnfc/(fs/2)放上去这属于加分项。5.4 验收时容易被问到的三个问题第一个问题是“为什么不用带通”如果1.8kHz噪声和语音频谱明显分离低通足够如果噪声落在语音频带中间低通会把语音一起削掉这时要改用自适应滤波或谱减法而不是一味堆阶数。第二个问题是“SNR提升为什么不明显”因为正弦噪声能量集中在单频低通能去掉它但语音本身也有1.8kHz以上的成分滤波后这些成分也被去掉。以干净语音y0为参考计算SNR时yf与y0的差异反而可能变大。所以SNR只是指标之一还要看时域波形和主观听感。第三个问题是“能不能改成实时处理”资源包里的filtfilt不能实时因为它是非因果的。要改成实时需要改用filter并接受相位延迟或者用FIR滤波器配合延迟补偿。在报告里把这个边界写清楚面试时非常加分。6. 把参数验证写成固定流程一个滤波器自检函数做课设和拆源码包最大的区别是拆包只要读懂课设要对自己写的东西负责。我后来给自己定了一条规矩任何滤波器交付之前先跑一遍自检函数。function verify_filter(b, a, fs) % 检查滤波器稳定性并画幅频和相频响应 if length(a) 1 ~all(abs(roots(a)) 1) error(滤波器不稳定存在极点超出单位圆); end [H, f] freqz(b, a, 1024, fs); subplot(2,1,1); plot(f, abs(H)); grid on; xlabel(频率 (Hz)); ylabel(|H(f)|); subplot(2,1,2); plot(f, angle(H)*180/pi); grid on; xlabel(频率 (Hz)); ylabel(相位 (度)); endroots(a)计算的是IIR分母的极点只要有一个极点落在单位圆外这个滤波器就不适合实际使用。freqz返回频率响应1024是频响计算点数fs用于把横轴换算成Hz。这个自检函数把稳定性和可视化一次做完省去每次手工敲freqz的麻烦。我的习惯是在加噪声脚本里加一行调用verify_filter(b, a, fs);如果报错说明参数不合适直接回第3.3节改阶数或截止频率。以前我总在阶数设得过高时看到输出毛刺原因往往是极点越界却没有快速定位的手段。从那以后每次设计滤波器我都强制先跑一遍自检再继续后面的处理效果确实稳定很多。希望帮到你。本文还有配套的精品资源点击获取