写论文的人估计都有同感:小波去噪、压缩感知、信号稀疏表示这一类题目,最后都绕不开一张对比图——几个经典信号在各种方法下重构/去噪后的波形和SNR表。而这张图里出镜率最高的,就是Blocks、Bumps和Doppler这三个名字。
这三个信号是Donoho和Johnstone在1994年前后提出的标准测试信号,出自他们关于小波收缩的经典论文Ideal Spatial Adaptation by Wavelet Shrinkage。后来随着WaveLab工具箱的流传,几乎成了信号处理领域约定俗成的"基准套餐"。无论你是做小波阈值去噪、压缩感知重构、稀疏字典学习,还是设计新的自适应滤波算法,都需要一组能代表不同类型信号特征的测试数据来验证算法的性能。Blocks模拟的是边缘突变信号,Bumps模拟的是稀疏峰值信号,Doppler模拟的是频率随时间变化的非平稳信号。用这三个信号跑一遍,基本就能看出一个算法对突变、脉冲、时变频谱三类问题的适应能力。
这篇文章我会从数学原理讲到纯Matlab实现,把这三个信号的生成代码、参数含义、校验方法以及实际使用中的坑一次说清楚。就算你手头没有WaveLab工具箱,也不需要联网下载任何第三方包,直接复制代码就能在自己机器上生成文献同款信号。
1. 为什么文献里的去噪对比实验总绕不开这三个信号
1.1 Donoho-Johnstone测试集在信号处理中的地位
先说点背景。90年代初,Donoho和Johnstone提出小波收缩(wavelet shrinkage)方法时,面临一个很实际问题:怎么客观评价一个新去噪算法到底好不好?如果用真实语音、真实心电信号,每个应用场景的数据特征完全不同,噪声类型也千差万别,很难横向量化比较。如果纯用白噪声或者正弦波,又太简单,无法体现算法的真实能力。
于是他们设计了一组人工合成信号,把实际信号中常见的几种结构特征提炼出来。这套信号的初衷就是当"标准试纸",让不同算法、不同参数、不同实现方式能在同样的输入上直接对比输出。后来这套信号被封装进WaveLab(斯坦福和莱斯大学联合开发的一个Matlab/Octave工具箱),跟着工具箱一起在学术界传播开来。现在你在IEEE、Signal Processing、IEEE Transactions等期刊上随便翻一篇关于去噪或重构的文章,实验部分大概率会出现这三个名字。
使用这套信号还有一个隐性的优势:它们不是某个具体应用场景的数据,没有专利和版权问题,任何人可以在任何论文中自由复现和使用。对做对比实验的人来说,这是很大的便利。
1.2 三个信号分别模拟了哪类真实信号特征
这三个信号虽然都是合成信号,但它们各自抓住了真实世界信号的一类核心特征,这也是它们三十年都没被淘汰的原因。
Blocks是分段常值信号,每个区间内幅值恒定,突然在某几个时间点发生跳变。它模拟的是地质分层信号、图像中的物体边缘、基因表达中的拷贝数变化这类带有明显突变边界的数据。这类信号的特点是"局部平滑但整体含有剧烈跳变",在小波域中只有少数几个大系数对应跳变位置,非常适合检验算法对边缘的保护能力。
Bumps则是若干尖锐脉冲的叠加,每一个脉冲都有自己的位置、高度和宽度。它模拟的是光谱信号中的特征峰、核磁共振谱线、雷达回波中的目标脉冲这类稀疏脉冲型数据。这类信号的特点是在时域上大部分区域接近零,只有少数位置有非零脉冲,适合检验算法对孤立峰值和微弱脉冲的检出能力。
Doppler是一个频率随时间增加而升高的振荡信号,而且幅度服从一个sqrt(t(1-t))的包络。它模拟的是多普勒雷达回波、地震波、语音信号中带有频率调制特性的非平稳成分。这类信号的最大难点在于频率不是恒定的,传统傅里叶分析在全局上很难准确刻画,而小波分析天然具有时频局部化能力,用它来检验算法的时频分辨率很合适。
1.3 从稀疏性的角度理解为什么是这三个组合
从数学上看,这三个信号共同构成了对小波基"稀疏表示能力"的三重考验。小波变换能把能量集中在少数系数上,这是小波去噪和压缩感知方法理论正确的前提。但不同信号在小波域的稀疏程度差异很大。
Blocks在Haar小波基下的表示是极度稀疏的,因为Haar小波的形状天然匹配分段常数函数,只要确定跳变位置和幅值,几乎就能完美重建。Bumps需要更高阶的小波基才能较好地稀疏化,因为它的脉冲形状更平滑,用Haar基表示会产生较多系数,而用Daubechies或Symlet小波则表现更优。Doppler虽然是非平稳信号,但小波包的时频分解可以把它在若干个尺度上展开,同样存在稀疏表示的可能性,只是对基函数的选择更敏感。
用这三个信号一起测试,本质上就是在测试算法对"稀疏性假设是否成立"这个前提的鲁棒性。如果你提出的去噪算法只在Blocks上效果好,在Bumps和Doppler上效果差,说明算法依赖了某种特殊的结构先验,在一般场景下可能并不可靠。这就是为什么审稿人往往要求实验部分同时给出这三个信号的对比结果。
2. 三个信号的数学原理:参数不是随便定的
2.1 Blocks的分段常值构造原理
Blocks信号的构造公式如下:
[ f(t) = \sum_{k=1}^{K} h_k \cdot \frac{1 + \text{sgn}(t - t_k)}{2}, \quad t \in [0, 1] ]
其中(t_k)是跳变位置,(h_k)是对应位置的跳变幅度,(\text{sgn}(\cdot))是符号函数。当(t < t_k)时,(\frac{1+\text{sgn}(t-t_k)}{2}=0);当(t > t_k)时,它等于1。所以每一项都相当于在时间轴上的某个位置"开启"一个常量电平,把所有项叠加起来就得到一个阶梯状的分段常值信号。
标准参数中一共用了11个跳变点,位置分别位于0.10、0.13、0.15、0.23、0.25、0.40、0.44、0.65、0.76、0.78、0.81这些时间点,对应的幅度为4、-5、3、-4、5、-4.2、2.1、4.3、-3.1、2.1、-4.2。注意幅度有正有负,所以信号不是单调递增的阶梯,而是上下反复跳变的不规则轮廓。之所以选这些"看起来不怎么整齐"的位置和幅度,是为了避免信号本身包含某种周期性或对称性——那样会让算法占便宜,比如恰好利用对称结构来改善效果,导致测试结果失真。
跳变点一共11处,这对算法来说有明确的挑战:到底能不能准确定位每一个跳变位置?跳变幅度差的信号段能否被保留?
2.2 Bumps的尖峰叠加原理
Bumps信号的构造公式是:
[ f(t) = \sum_{k=1}^{K} h_k \cdot \left(1 + \left|\frac{t - t_k}{w_k}\right|\right)^{-4} ]
每一项都是以(t_k)为中心的钟形脉冲。分母中有4次方,这让脉冲的衰减速度非常快,远离中心的区域迅速归零,只在中心附近留下一个尖锐的"峰"。参数(h_k)控制峰高,(w_k)控制峰的宽度,(w_k)越小峰越窄越尖锐。
标准参数里,位置和Blocks共用了同一组(t_k),但高度全部取为正数(4、5、3、4、5、4.2、2.1、4.3、3.1、2.1、4.2),宽度参数为0.005、0.005、0.006、0.01、0.01、0.03、0.01、0.01、0.005、0.008、0.005。可以看到大部分宽度都在0.01以下,相比整个区间[0,1]来说,这些峰是非常窄的,整体看起来就是在一段几乎平坦的背景上冒出一串尖刺。其中0.40处的峰宽度为0.03,是其中较宽的一个,这故意制造了一个尺度差异,考察算法对不同宽度脉冲的适应能力。
Bumps的4次方衰减有一个特点:它介于高斯衰减(指数速度)和Lorentzian长尾衰减(平方速度)之间。4次方衰减不会像高斯那么快,从而保留了一定的"旁瓣"信息;也不会像Lorentzian那么慢,导致峰之间互相粘连。
2.3 Doppler的非平稳频率调制原理
Doppler信号的构造公式更简洁:
[ f(t) = \sqrt{t(1-t)} \cdot \sin\left(\frac{2\pi(1+\epsilon)}{t+\epsilon}\right), \quad \epsilon = 0.05 ]
这个式子从多普勒现象中获得启发:当一个辐射源相对观测者运动时,接收到的信号频率会随距离变化而变化。公式里的(\sin)参数是((2\pi(1+\epsilon))/(t+\epsilon)),当(t)从0向1增大时,分母从0.05增大到1.05,所以整个分式的值从约(2\pi\times21)下降到(2\pi\times1),也就是说瞬时频率在信号过程中从高到低变化,跨越了大约一个数量级。
前面的系数(\sqrt{t(1-t)})是一个抛物形包络:在(t=0)和(t=1)处取值为0,在(t=0.5)处取得最大值0.5。这个包络让信号在两端趋于零,避免了边界上出现突兀的跳变,同时也给信号赋予了随时间的幅度调制特性。
Doppler信号最具挑战性的一点在于,它具有明显的非平稳性:频率持续变化,且变化的跨度很大。任何全局性的变换手段在这种信号面前都容易"顾此失彼";而理想的去噪或重构算法应该能够在不丢失高频细节的情况下,同时保存低频部分的整体形态。
2.4 Wavelab标准参数:为什么我用这组数
你自己设计测试信号的时候,可能会疑惑:跳变位置为什么选0.13、0.23这些数而不是0.1、0.2这种整十数?这其实是"反规律性"设计。测试信号最怕的是隐藏某种均匀分布或对称规律,一旦信号结构有规律,某些善于利用特定结构先验的算法就会"作弊式"地取得好成绩,而真正通用的性能反而被掩盖了。所以标准参数刻意把跳变位置设计得不规则,让不同跳变之间没有固定间隔,这样任何算法都不容易利用位置上的先验信息。
这组参数另一个重要特点是尺度不变性。无论你生成多少采样点,这些位置始终在归一化的[0,1]区间内。采样点数n只是决定了每个区间内的点数密度,并不会改变跳变点的时间位置和幅值。所以在不同采样率下生成的信号之间可以横向比较,这也是它能成为标准测试集的原因之一。
3. 不依赖Wavelab工具箱:三个函数的完整Matlab实现
3.1 编写前的约定:向量方向与函数封装
在实际动手写代码之前,我先说明几个约定,这些都是实际使用中最容易出问题的地方。
第一,输出信号统一用列向量。Matlab里linspace默认生成行向量,但大多数信号处理函数(如wdenoise、dwtenc、dct)和滤波器设计函数对列向量更友好。所以我们用t = linspace(0, 1, n)';这种转置操作把t变成列向量,后续生成的信号自然就是列向量。
第二,函数输入参数给默认值。写成function y = genBlocks(n),然后用if nargin < 1, n = 1024; end做默认值处理。这样当你只调用genBlocks而不传参时,也能生成一个默认长度的信号,调试时方便。
第三,向量化计算。Matlab里sin、sqrt、sign等函数对向量输入自动按元素操作。我们构造信号时不要用循环去逐点计算;但对Blocks和Bumps来说,每个信号分量需要循环叠加,这个循环是在分量维度上的,循环次数极少,完全没有性能压力。
3.2 genBlocks.m完整代码
function y = genBlocks(n) % GENBLOCKS 生成Donoho-Johnstone Blocks标准测试信号 % y = genBlocks(n) 生成长度为n的Blocks信号 % 默认n = 1024 % % 参考: Donoho & Johnstone (1994), Biometrika 81(3):425-455 if nargin < 1 n = 1024; end t = linspace(0, 1, n)'; % 跳变位置(归一化时间坐标) pos = [0.10, 0.13, 0.15, 0.23, 0.25, 0.40, 0.44, 0.65, 0.76, 0.78, 0.81]; % 跳变幅度 hgt = [4, -5, 3, -4, 5, -4.2, 2.1, 4.3, -3.1, 2.1, -4.2]; % 用符号函数构造阶梯信号 y = zeros(n, 1); for k = 1:length(pos) y = y + hgt(k) * (1 + sign(t - pos(k))) / 2; end end这段代码的核心只有一行:y = y + hgt(k) * (1 + sign(t - pos(k))) / 2;。sign(t - pos(k))在t小于pos时返回-1,大于pos时返回1,所以(1+sign)/2就是逻辑上的"开关"函数。之所以用1+sign再做除法,而不是直接用t > pos(k),是因为sign函数同时处理了所有采样点,整个表达式是一个向量计算,后续维护起来很清晰。
3.3 genBumps.m完整代码
function y = genBumps(n) % GENBUMPS 生成Donoho-Johnstone Bumps标准测试信号 % y = genBumps(n) 生成长度为n的Bumps信号 % 默认n = 1024 % % 参考: Donoho & Johnstone (1994), Biometrika 81(3):425-455 if nargin < 1 n = 1024; end t = linspace(0, 1, n)'; % 脉冲中心位置 pos = [0.10, 0.13, 0.15, 0.23, 0.25, 0.40, 0.44, 0.65, 0.76, 0.78, 0.81]; % 脉冲高度 hgt = [4, 5, 3, 4, 5, 4.2, 2.1, 4.3, 3.1, 2.1, 4.2]; % 脉冲宽度 wth = [0.005, 0.005, 0.006, 0.01, 0.01, 0.03, 0.01, 0.01, 0.005, 0.008, 0.005]; y = zeros(n, 1); for k = 1:length(pos) y = y + hgt(k) ./ (1 + abs((t - pos(k)) ./ wth(k)).^4); end endBumps的计算中有一个容易写错的地方:分母是4次方而不是平方。有些博客和论坛上给的是平方版本,其实那已经偏离了原始定义。4次方的好处是峰腰部收窄、尾部衰减比高斯更慢,这样峰的形态更接近光谱类信号的实际特征。如果你做对比实验,一定要用4次方版本,否则生成的波形和文献中的标准图会不同,审稿人对图有印象,一眼就能发现。
3.4 genDoppler.m完整代码
function y = genDoppler(n) % GENDOPPLER 生成Donoho-Johnstone Doppler标准测试信号 % y = genDoppler(n) 生成长度为n的Doppler信号 % 默认n = 1024 % % 参考: Donoho & Johnstone (1994), Biometrika 81(3):425-455 if nargin < 1 n = 1024; end t = linspace(0, 1, n)'; epsilon = 0.05; y = sqrt(t .* (1 - t)) .* sin((2 * pi * (1 + epsilon)) ./ (t + epsilon)); endDoppler的实现最简洁,但要注意一点:sqrt(t .* (1 - t))中用的点乘,因为t是一个列向量,这里是要逐元素相乘,不是矩阵乘法。如果漏掉点号,轻则报维度错误,重则在某些自动广播机制下得到完全错误的矩阵,返回一个n×n的结果,后续绘图和计算都得乱。我总是建议写完一行涉及向量的式子后,先看一眼size(y)是否正确。
3.5 主脚本:一次性生成、绘图、保存
下面给一个可以直接跑通的主脚本,把三个信号生成后放在同一个figure里对比展示,并保存到.mat文件备用。
%% 生成三个标准测试信号并可视化 clear; close all; clc; n = 1024; t = linspace(0, 1, n)'; blockSig = genBlocks(n); bumpsSig = genBumps(n); dopplerSig = genDoppler(n); figure('Color', 'w', 'Position', [100 100 820 900]); subplot(3,1,1); plot(t, blockSig, 'LineWidth', 1.2); ylabel('幅值'); title('Blocks 标准测试信号'); grid on; xlim([0 1]); subplot(3,1,2); plot(t, bumpsSig, 'LineWidth', 1.2); ylabel('幅值'); title('Bumps 标准测试信号'); grid on; xlim([0 1]); subplot(3,1,3); plot(t, dopplerSig, 'LineWidth', 1.2); ylabel('幅值'); xlabel('时间 t'); title('Doppler 标准测试信号'); grid on; xlim([0 1]); %% 保存为MAT文件 save('test_signals_1024.mat', 't', 'blockSig', 'bumpsSig', 'dopplerSig'); disp('信号已生成并保存到 test_signals_1024.mat');如果你和我一样经常在论文里引用这些信号,我强烈建议你直接把上面三个函数保存成genBlocks.m、genBumps.m、genDoppler.m放到一个专门的目录里,加进Matlab路径。之后任何脚本里直接调用,不用每次复制粘贴一遍。
4. 结果校验、可视化与常见改动
4.1 如何确认生成的信号和文献一致
生成完后别急着跑算法,先校验一下。最容易做的校验是检查信号的关键统计量,因为文献和标准代码里这些信号的特征值是有共识的。
以n=1024为例,你可以算一下每个信号的最大值、最小值和标准差。Blocks的最大值大约在8.3左右(因为最大一段是从0.40跳到5后叠加了前方的2.1,再加上0.44处的-4.2跳变前是2.1+5=7.1,到0.65处又加4.3达到11.4附近,具体要按叠加关系算);Bumps的最大值大约在5左右(最尖的峰高度为5);Doppler的最大值约在0.45左右。不同机器上浮点运算的尾数可能有微小差异,但前三位有效数字应该一致。
还有一个更可靠的校验方法:用Matlab的wavedec对生成的Blocks做一级Haar小波分解,你会发现除了对应跳变位置的系数较大之外,大部分细节系数为零。这正是Blocks信号在小波域极度稀疏的直接体现。如果你分解后细节系数到处都是中等大小的值,那你的信号生成多半有问题。
最后,可以和自己手头的文献对照。翻开任何一篇提到Blocks或Bumps的论文,它们展示的信号图像基本都是同一组位置、同一组幅度、同一组宽度画出来的。把子图画出来和论文里的图比一比,形状大致吻合就说明代码没问题。
4.2 改变信号长度n时的注意事项
很多初学者会问:我改成n=2048、4096,或者n=1000,信号还准不准?答案是只要保证t = linspace(0,1,n),n可以任意取,跳变点和幅值都不受影响。这是因为位置参数已经归一化到[0,1]区间,n只影响每个跳变间隔内的采样密度。
但有一个细节要注意:当n比较小时(比如n=64或n=128),Blocks和Bumps的部分峰可能只被几个采样点覆盖。Bumps里最窄的峰宽度只有0.005,在n=128时,0.005×128 = 0.64个采样点,也就是说这个峰可能恰好落在两个采样点之间,导致峰形严重失真。如果你需要在低采样率下使用Bumps信号,建议考虑增加n,或至少确认最窄峰被至少3到5个采样点覆盖。
Doppler对采样率的要求更高。因为它的瞬时频率在靠近t=0时非常高,要准确刻画高频部分,n至少要取到1024或2048。如果你用n=128,你生成的Doppler其实已经严重欠采样,高频振荡部分几乎看不出来。这也是为什么标准实验里一般默认n=1024,如果要做多尺度分析,取2048或4096更稳妥。
4.3 添加噪声与归一化处理
标准测试信号经常和加性高斯白噪声配合使用。添加噪声的公式很简单:
rng(2024); % 固定随机种子,保证结果可复现 sigma = 0.5; % 噪声标准差 noisyBlock = blockSig + sigma * randn(n, 1);这里有一个新手经常忽略的点:randn每次运行结果不同。如果你要对比不同算法,或者在论文中给出某个具体的SNR结果,必须在调用randn之前固定随机种子,用rng(2024)或者rng(0)。否则同样的代码跑两次,输出的SNR值会略有差异,读者也无法复现你的实验。
关于噪声强度的设置,文献里常见的有两种方式。一种是直接指定sigma,比如sigma=0.1、0.5、1.0;另一种是先把信号归一化到单位标准差,再按目标SNR计算sigma。第二种方式更公平,因为不同信号的原始幅度差异很大,直接给固定sigma可能出现某个信号"噪声几乎淹没信号",另一个信号"噪声几乎看不见"的情况。常见的做法是:
x = blockSig / std(blockSig); % 单位标准差化 sigma = sqrt(1 / SNR_linear); % 根据目标SNR计算噪声标准差 noisyX = x + sigma * randn(n, 1);这里SNR_linear由SNR_dB换算得到:SNR_linear = 10^(SNR_dB / 10)。经过单位标准差化后,给定SNR=10dB就表示噪声标准差约为0.316,此时噪声的功率是信号的十分之一。这样处理可以让不同信号、不同采样长度之间的信噪比都有统一含义。
4.4 保存为mat文件的方法
保存数据用save函数,我推荐保存成-v7.3版本格式的.mat文件,如果数据量很大或者要在不同Matlab版本之间共享:
save('test_signals_1024.mat', 't', 'blockSig', 'bumpsSig', 'dopplerSig', '-v7.3');如果只是自己本地用,默认格式就行。注意变量名要带引号,写成一整个字符串或字符向量元胞数组都可以。读取时用load('test_signals_1024.mat'),读进来的变量会以保存时的名字出现在工作区。
如果你和别人协作,对方的Matlab里没有你的函数文件,那么直接把信号变量保存成.mat传过去是最省事的。当然,也可以只保存文本格式:
writematrix([t, blockSig, bumpsSig, dopplerSig], 'test_signals.csv');这样就算对方不用Matlab,也能用Excel、Python或Origin打开。我在多语言协作项目里通常两者都留一份。
5. 从生成到实战:在典型算法验证工作流中的应用
5.1 小波去噪的基准对比
生成这三个信号之后,最直接的用途就是小波去噪对比。以Blocks信号为例,下面是完整的去噪流程框架:
rng(2024); n = 1024; t = linspace(0, 1, n)'; x = genBlocks(n); x = x / std(x); % 单位标准差化 sigma = sqrt(1 / 10^(10/10)); % 目标SNR=10dB对应的噪声标准差 y = x + sigma * randn(n, 1); % 含噪信号 % 用默认软阈值小波去噪 xd = wdenoise(y, 'Wavelet', 'db4', 'DenoisingMethod', 'Bayes'); % 计算去噪前后的SNR SNR_in = 10 * log10(sum(x.^2) / sum((y - x).^2)); SNR_out = 10 * log10(sum(x.^2) / sum((xd - x).^2)); fprintf('输入SNR = %.2f dB, 输出SNR = %.2f dB\n', SNR_in, SNR_out);把同样的流程分别套到Bumps和Doppler上,就能得到一张算法在三个信号上的对比表。这也是论文实验部分最基础的呈现方式。为了公平,去噪参数(小波基、分解层数、阈值规则)在三个信号上应保持一致,除非你有明确理由说明某个算法需要针对信号类型调参——那样的话审稿人会要求你解释这种依赖性的合理性。
5.2 压缩感知重构实验
这三个信号也是压缩感知实验的常客。在压缩感知理论中,信号需要在某个变换域稀疏。Blocks在Haar小波基下稀疏,Bumps在Daubechies小波基下较稀疏,Doppler在小波包或DCT基下有一定稀疏性。你可以分别尝试用不同的稀疏基表示这三个信号,观察表示系数的衰减速度。
一个典型实验设计是:用DCT字典或小波字典作为稀疏基,对信号做随机欠采样(只保留一部分傅里叶系数或随机观测矩阵),然后用OMP或BPDN算法重构,画出重构误差随采样率变化的曲线。这时三个信号的差异会体现得很明显——Blocks在很低采样率下就能恢复得不错,Doppler则需要更多观测才能达到同等重构精度。
DCT时,需要注意信号长度最好选2的幂,方便小波变换的层数选择。比如n=1024正好是2的10次方,最多可以分解10层。如果你取n=1000,虽然也能做,但分解层数和边缘拼接上会有一些不必要的麻烦。
5.3 与深度学习信号处理任务的衔接
这几年深度学习方法大火,很多人也把这些标准测试信号搬进了深度学习实验里。如果你的研究涉及深度学习去噪或超分辨,你可以批量生成大量不同长度、不同噪声强度、不同信噪比的Blocks、Bumps和Doppler样本,作为训练集和测试集:
rng(42); numSamples = 1000; dataset = cell(numSamples, 1); labels = cell(numSamples, 1); n = 1024; sigmas = linspace(0.05, 0.5, numSamples); % 噪声强度从低到高变化 for idx = 1:numSamples % 随机选一种信号类型 switch randi(3) case 1 clean = genBlocks(n); case 2 clean = genBumps(n); case 3 clean = genDoppler(n); end clean = clean / std(clean); noisy = clean + sigmas(idx) * randn(n, 1); dataset{idx} = noisy(:); labels{idx} = clean(:); end save('trainset.mat', 'dataset', 'labels', '-v7.3');因为三个信号都有解析表达式,你可以无限生成样本而不需要做数据增广,这对训练深度学习模型是很大的优势。另一方面也要提醒,用合成信号训练出的模型一度被认为在真实数据上泛化能力欠佳,所以这类实验结果比较适合论证方法在受控条件下的有效性。如果要论证实际应用价值,还需要在真实数据集上补充验证。
6. 生成过程中容易踩的坑与排查思路
6.1 断点处sign函数的取值问题
sign(0)在Matlab里返回0,不是1。这意味着如果你让某个跳变位置恰好落在某个采样点上,比如t=0.10恰好等于某个采样值,那么这一项在断点处的贡献是hgt/2,而不是hgt或0。由于linspace(0,1,n)产生的采样点一般不会恰好等于0.10(浮点数的0.1本来就不是精确值),实际使用中很少遇到真正等于0的情况,但理论上你要知道这个行为。
如果你希望信号在断点处的值采用“断点左侧值”,建议把(1 + sign(t - pos(k))) / 2改写成逻辑索引:
mask = t > pos(k); y(mask) = y(mask) + hgt(k);这样断点处严格从0跳到hgt,更接近理想化阶梯信号。但这种写法会让代码可读性变差,如果你只是做常规实验,用sign版本完全够用。我自己的经验是,在标准测试信号上去噪重构时,断点处单点的取值对整体SNR影响通常小于0.1dB,不值得为这个微差牺牲代码的简洁性。
6.2 向量维数与绘图失败
最常见并且让人摸不着头脑的报错是“矩阵维度必须一致”。这通常由两种原因造成。一是你用冒号运算符生成时间轴:t = 0:0.001:1,此时t是行向量,而信号函数返回的是列向量,两者相加就出错了。解决办法是统一使用转置,或者统一用linspace,我建议统一用linspace(0,1,n)'生成列向量。
二是你在调用genBlocks(n)时传的n是向量而不是标量,比如不小心传了1:1024,导致linspace(0,1,1:1024)报错。这种问题一般出现在你把代码封装成循环、变量名冲突的场景里。排查时用dbstop if error打开错误模式,查看工作区中每个变量的size,一眼就能定位。
还有一种看似不报错但结果错误的场景:你把randn噪声写成了矩阵加法y = x + randn(n,n),Matlab在某种广播机制下可能不报错但生成一个二维数组,后续计算虽然不报错但结果完全是错的。所以我每次生成完信号都会习惯性检查size(blockSig),确认是[1024, 1]而不是[1024, 1024]。
6.3 除零、NaN与Doppler信号的端点行为
Doppler公式在t=0时有t + epsilon = 0.05,所以不会出现除零。但你会发现在t靠近0时,正弦函数的参数极大,振荡非常剧烈,这也意味着t=0附近是采样最困难的地方。如果在某些非常规条件下你让epsilon=0,那么t=0处的分母会变成0,直接产生NaN。所以建议理解这个公式时保持epsilon为0.05不变。
另一个容易忽略的点是,Matlab中sqrt(t .* (1 - t)),如果t包含0和1,那么包络在两端都是0,整个信号在端点的取值就是0,没有问题。但如果你的t不是严格的0到1,比如不小心用了t = linspace(0, 1, n+1)然后再去掉最后一项,导致t的区间变成[0,1)或(0,1],信号边界的包络值会稍有不同。虽然对结果影响不大,但如果你特别在意和文献严格一致,请务必确认t(1)=0且t(end)=1。
6.4 随机噪声可复现性问题
我在前面已经强调过rng的重要性,这里再展开说一个实战场景。假设你对比了两篇论文给出的去噪后SNR,却发现对不上。除了算法实现不同,一个可能性很大的原因就是两篇论文用了不同的随机噪声实现。A论文用了randn(n,1)生成噪声,B论文用了randn(n,1)但之前调用了别的随机数生成函数,或者A和B的Matlab版本不同,导致randn内部算法略有变化。
为了让读者能够精确复现你的结果,我建议在论文的实验描述中明确写出“使用Matlab R2021a,随机种子固定为2024,噪声标准差按单位标准差信号计算得到”。同时在代码中统一用rng(seed)且每个独立实验都用不同的已知seed。这样哪怕别人在别的版本上跑,大概率也能复现出近似的数值结果,至少趋势是稳定的。
如果做批量化实验,还有一个技巧是保存每次实验的随机状态:
seed = 2024; rng(seed); noiseState = rng; % 保存当前随机流状态 % 第一次实验... rng(noiseState); % 恢复,保证第二实验使用相同的噪声序列这样可以在多个实验之间使用完全相同的噪声序列,避免噪声本身的随机差异干扰实验结果对比。这个方法在我自己写对比实验时经常使用,比单纯固定seed更灵活。
生成这三个标准测试信号本身并不难,难的是确保生成结果和文献一致、在不同实验之间可复现、并且能把它们有效嵌入到算法验证流程里。把这三个函数存好,把rng的用法养成习惯,你之后做小波去噪、压缩感知、稀疏表示相关实验的时候,就能把更多时间花在核心算法本身,而不是反复折腾这些基础信号的构造细节上。希望这篇文章能帮你少走这些弯路。