简介这一Matlab时频分析程序包面向信号处理领域的研究者与电子通信专业学生提供可交互指定Cohen类时频分布的计算工具覆盖WVD、CWD、PWVD等典型分布形式适用于非平稳信号的瞬时频率提取、时频特征分析与雷达、语音等场景的信号解译。包体十分轻量共139个文件其中以134个m源文件为核心覆盖分布计算、结果展示与示例程序另含少量mat数据、asv自动备份和txt说明压缩包仅1.01MB下载后可快速在MATLAB环境中运行。程序内整合了tfdemo系列演示脚本与tfrview、tfrqview等可视化模块可直观对比不同核函数下的时频聚集度辅助理解参数设置对分布结果的影响。目前已有560人浏览学习适合作为课程设计、竞赛准备或科研预研阶段的参考代码库。1. 时频图不是只有STFTCohen类核函数切换才是信号分析的关键很多人在振动信号故障诊断、语音特征提取或雷达回波分析里先用短时傅里叶变换画一张时频图然后发现窗长怎么调都别扭窗短了频率模糊窗长了瞬态被抹平。STFT本质上是给信号加一个固定窗窗函数一旦选定时间分辨率和频率分辨率就被锁死。而Cohen类时频分布用双线性变换和核函数把WVD、CWD、PWVD统一到一个框架里核函数一变时频图的“摊平方式”就变了可以在交叉项抑制和能量聚集度之间主动做选择。这套MATLAB程序把WVD、CWD、PWVD封装成可切换的Cohen类分布计算工具适合做非平稳信号分析、滚动轴承故障特征提取以及给深度学习模型准备时频图像输入。无论你刚看完matlab下载安装教程还是已经在用matlab做图像处理只要把扩展名为.m的文件放进当前目录就能跑。2. Cohen类分布统一框架WVD、CWD、PWVD的核函数选型2.1 WVD为什么会有交叉项核函数等于1的代价Wigner-Ville分布是所有Cohen类分布的原点。对实信号 x(t)WVD定义在瞬时自相关 R_x(t,τ)x(tτ/2)x*(t-τ/2) 的傅里叶变换上W_x(t,f)∫_{-∞}^{∞} x(tτ/2)x^*(t-τ/2)e^{-j2πfτ}dτ由于是信号与自身延迟相乘再变换双线性结构会让两个分量在时频平面中间位置产生交叉项。交叉项是虚的干扰能量幅值往往能达到真实分量的两倍但它也是WVD高分辨率的代价。核函数为1时不做任何平滑时频聚集度最高交叉项也最明显。用这套程序时如果只看WVD结果看到中间出现一大堆条纹状振荡不要先怀疑代码有问题而是要知道这是核函数没有抑制交叉项的正常表现。2.2 CWD的核函数如何压交叉项控制σ即控制惩罚范围Choi-Williams分布的核函数是高斯指数形式φ(θ,τ)exp(-θ²τ²/σ)这里的θ是频率延迟τ是时间延迟。核函数对离原点远的(θ,τ)区域做指数衰减等效于在模糊域里过滤掉远离坐标轴的交叉项部分。σ是缩放因子σ越小核函数越集中于原点附近抑制交叉项越强但自项也会被平滑σ趋于无穷大时核函数趋向1退化成WVD。实际调试时σ在1附近开始尝试交叉项重就调小到0.1分辨率不够就调大到10。2.3 PWVD是在时间和频率两个方向分别平滑伪Wigner-Ville分布PWVD的做法是对瞬时自相关加时间窗h(τ)相当于只沿频率方向做平滑去掉瞬时自相关在时间维上的多普勒响应。所以说PWVD是“伪”的是因为它只做了一半工作交叉项只被部分抑制。想要两个方向都平滑还要再引入时间平滑窗g(s)得到平滑伪WVDSPWVD。这套程序中的tfrspaw.m如果从命名习惯看很可能是带平滑窗的伪WVD核函数实现内部通过type参数切到不同的核函数形态。下表列出三种分布在程序里对应核函数的关键差异方便你在修改type时心里有底分布类型核函数形式主要调节参数交叉项时频聚集度WVDφ1无严重最高CWDexp(-θ²τ²/σ)σ明显抑制中高PWVD时间窗h(τ)窗长部分抑制中等2.4 为什么要把核函数设计成可切换参数如果每种分布都写一遍双线性积分代码重复率高而且不同核函数之间数值误差不一致对比实验结果时很难说是核函数的差异还是积分实现的差异。把核函数抽成一个统一入口外部只传分布类型字符串和对应参数内部走同一套自相关计算、同一套离散化采样才能保证WVD、CWD、PWVD三者结果可正交对比。这也是这个项目里tfrspaw.m和tfrspbk.m分开的原因核心计算与核函数表分开后续换核函数不需要动主流程。3. 这套MATLAB程序的调用结构从demo到指定type参数3.1 文件清单里每一份是干什么的拿到压缩包后先不要急着运行把文件分一下类。tolowcase.asv是MATLAB自动保存的备份文件内容与某个函数相同可以删除不影响运行。contents.m是工具箱说明文件MATLAB在命令行执行help contents时会展示这个目录下各函数的作用。tfdemo2.m到tfdemo5.m是一组演示脚本分别展示不同信号模型下的时频分布结果。tfrview.m和tfrqview.m是可视化工具前者画二维时频图后者画时频剖面或量化指标。tfrspaw.m和tfrspbk.m是核心计算函数一个负责指定Cohen类分布的计算另一个可能是某种备用核函数实现。打开contents.m可以看到函数清单。如果某个函数名右侧的描述带有“Cohen class”“time-frequency distribution”字样那就是你要改的入口。常见做法是先运行一个demo确认路径没有问题cd /your/path/cohen_tfd which contents.m tfdemo2这段命令先切换目录再用which确认文件真的在搜索路径上最后启动演示脚本。很多matlab下载安装教程都强调过set path的重要性这里直接调用可以避免找不到函数的报错。3.2 用switch包装一个可指定type的入口函数项目描述里说“可指定Cohen类分布类型WVDCWDPWVD”实际使用时不建议每次都去改tfrspaw.m源码。我一般会写一个薄包装函数把分布类型字符串解析成不同分支转到对应的核函数调用。下面这段代码可以作为模板function [tfr, t, f] cohen_type(x, fs, type, varargin) % x: 单通道信号fs: 采样率type: WVD/CWD/PWVD N length(x); t (0:N-1) / fs; f linspace(-fs/2, fs/2, N); switch upper(type) case WVD [tfr, t, f] tfrspaw(x, t, N, 0, 0); case CWD sigma varargin{1}; % 例如 sigma 1 [tfr, t, f] tfrspaw(x, t, N, sigma, 0); case PWVD win_len varargin{1}; % 时间窗长度奇数 [tfr, t, f] tfrspaw(x, t, N, 0, win_len); otherwise error(Unknown distribution type: %s, type); end end这段代码的逻辑是先根据采样率和信号长度构造时间轴和频率轴然后用switch语句把分布类型映射到tfrspaw的第四个和第五个参数上。第四个参数这里借用为CWD的核函数缩放因子σ第五个参数借用为PWVD的时间窗长度。具体签名要以你打开的tfrspaw.m头部注释为准有的版本会直接用type,param这种参数顺序。运行时注意varargin顺序CWD必须传σPWVD必须传窗长缺参数时应该在switch分支内直接报错而不是让NaN传到傅里叶积分里。3.3 参数怎么传才不踩坑从文件名看这套程序可能有多个demo脚本但每个demo里的信号长度和采样率都不一样。传参数时最常踩的坑有三个第一个是频率轴没有做fftshift导致时频图中心在四角第二个是CWD的σ不是越大越好σ大到100以上时核函数近似为1实际得到的就是WVD第三个是PWVD的窗长必须远小于信号长度否则窗内截断效应会引入虚假振荡。下面是个更安全的调用方式x randn(1, 1024); fs 1024; sigma 1; win_len 63; [tfr_wvd, t, f] cohen_type(x, fs, WVD); [tfr_cwd, t, f] cohen_type(x, fs, CWD, sigma); [tfr_pwvd, t, f] cohen_type(x, fs, PWVD, win_len);运行后检查tfr_wvd和tfr_cwd矩阵尺寸是否一致如果不一致说明tfrspaw内部对不同type做了不同下采样。此时不要直接做矩阵差值先size看维度再决定是插值还是截断。很多matlab教程里直接imagesc(tfr)如果矩阵维度不对就会提示赋值维度不匹配。4. 实测多分量信号下比较WVD、CWD、PWVD4.1 生成带跳变调频的实验信号要对比三种分布不能用单频正弦否则交叉项显示不出来。我一般构造一个两个线性调频分量叠加的信号一个频率从50Hz扫到300Hz另一个从250Hz扫到80Hz中间有一段接近交叉。下面是生成信号的代码fs 1024; t 0:1/fs:2-1/fs; f1 50 (300-50) * t / 2; phi1 2 * pi * cumsum(f1) / fs; f2 250 (80-250) * t / 2; phi2 2 * pi * cumsum(f2) / fs; x cos(phi1) cos(phi2);这里用cumsum而不是直接对f1积分是因为瞬时频率本身就是相位的一阶导累计求和能保证相位连续。两个分量在时间中点附近频率接近会产生显著交叉项适合暴露WVD的缺点。4.2 运行三种分布并用图像对比交叉项把信号送入前面写的cohen_type函数然后分别用imagesc显示[tfr_wvd, t_ax, f_ax] cohen_type(x, fs, WVD); [tfr_cwd, ~, ~] cohen_type(x, fs, CWD, 1); [tfr_pwvd, ~, ~] cohen_type(x, fs, PWVD, 63); figure; subplot(1,3,1); imagesc(t_ax, f_ax, abs(tfr_wvd).^2); axis xy; title(WVD); subplot(1,3,2); imagesc(t_ax, f_ax, abs(tfr_cwd).^2); axis xy; title(CWD); subplot(1,3,3); imagesc(t_ax, f_ax, abs(tfr_pwvd).^2); axis xy; title(PWVD);axis xy这行很重要默认imagesc会把y轴按矩阵行方向倒置不执行的话频率从高频到低频排列容易把上升扫频看成下降扫频。运行后你会看到WVD图中两个分量中间出现一串快速振荡的干扰条纹CWD图中振荡幅度明显下降但真实分量边缘也变模糊PWVD图介于两者之间靠近时间轴方向的部分初始段有轻微拖尾。4.3 交叉项定量指标和常见误判只看图像容易被视觉骗过去。可以计算一个简单的交叉项占比对时频能量矩阵做阈值分割统计低于最大幅值10%的像素能量占比占比越高说明交叉项越重ratio_wvd cross_ratio(tfr_wvd); ratio_cwd cross_ratio(tfr_cwd); ratio_pwvd cross_ratio(tfr_pwvd); function r cross_ratio(tfr) E abs(tfr).^2; Emax max(E(:)); r sum(E(E 0.1 * Emax)) / sum(E(:)); end这个比例不是严格意义上交叉项能量因为真实分量的低幅值拖尾也会被算进去但它能快速反映能量是否被分散。对同样的信号通常ratio_wvd ratio_pwvd ratio_cwd。如果你算出来CWD的比值反而更高先检查σ是否过小或者信号本身幅值波动太大。另一种常见误判是拿analyitic信号去算三种分布但忘了对实数信号做Hilbert变换WVD对实数信号会产生负频率镜像看起来像“多了一个分量”实际上不是代码bug。使用这些.m文件时如果demo脚本里信号直接是复数就不要额外取实部。5. 用tfrview和tfrqview把核参数从肉眼调参变成定量复验运行完demo后不要满足于imagesc画出来的亮图。tfrview.m和tfrqview.m才是这个程序包里真正提升效率的部分。tfrview.m打开一个交互式时频图光标移动时会显示当前点对应的时间和频率值还能用鼠标拖拽放大局部区域适合观察交叉项振荡分布在哪一片频率区间。tfrqview.m则展示指定时间点的频率切片或指定频率的时间切片还能计算时频聚集度指标比如Rényi熵。我一般用下面这组步骤复验三种分布的参数先运行一次tfrview手动确认交叉项集中在哪个区域然后回到命令窗口用tfrqview对比不同σ下的Rényi熵。Rényi熵公式里的q阶信息量对时频能量的集中度非常敏感熵值越低说明能量越集中在少数时频点时频聚集度越高。CWD的σ从0.1到10步长0.2每个σ算一次时频矩阵和熵值收敛到最小熵的σ就是当前信号下的较优核参数。这个流程把“多试几个参数看哪张图亮”的玄学问题变成了可重复的数值选择。如果要在脚本里批量做这个操作可以直接复用计算得到的时频矩阵不需要再调用GUIsigma_list 0.1:0.2:10; renyi zeros(size(sigma_list)); for i 1:length(sigma_list) [tfr_cwd, ~, ~] cohen_type(x, fs, CWD, sigma_list(i)); p abs(tfr_cwd).^2; p p / sum(p(:)); renyi(i) -log2(sum(p.^3)); end [~, best_idx] min(renyi); best_sigma sigma_list(best_idx);这段代码计算三阶Rényi熵p.^3突出了能量高点的权重所以熵值下降意味着时频平面出现了更尖锐的峰。注意这里没有加任何平滑如果时频矩阵维度过大直接对每个σ做全矩阵计算会占用较多内存可以先把矩阵缩放到256×256再算。看到best_sigma落在0.5到4之间时说明程序运行正常如果最佳σ落在边界说明信号含噪声过重需要先做带通滤波而不是继续调核函数参数。本文还有配套的精品资源点击获取
MATLAB实现Cohen类时频分布：WVD/CWD/PWVD核函数切换指南
简介这一Matlab时频分析程序包面向信号处理领域的研究者与电子通信专业学生提供可交互指定Cohen类时频分布的计算工具覆盖WVD、CWD、PWVD等典型分布形式适用于非平稳信号的瞬时频率提取、时频特征分析与雷达、语音等场景的信号解译。包体十分轻量共139个文件其中以134个m源文件为核心覆盖分布计算、结果展示与示例程序另含少量mat数据、asv自动备份和txt说明压缩包仅1.01MB下载后可快速在MATLAB环境中运行。程序内整合了tfdemo系列演示脚本与tfrview、tfrqview等可视化模块可直观对比不同核函数下的时频聚集度辅助理解参数设置对分布结果的影响。目前已有560人浏览学习适合作为课程设计、竞赛准备或科研预研阶段的参考代码库。1. 时频图不是只有STFTCohen类核函数切换才是信号分析的关键很多人在振动信号故障诊断、语音特征提取或雷达回波分析里先用短时傅里叶变换画一张时频图然后发现窗长怎么调都别扭窗短了频率模糊窗长了瞬态被抹平。STFT本质上是给信号加一个固定窗窗函数一旦选定时间分辨率和频率分辨率就被锁死。而Cohen类时频分布用双线性变换和核函数把WVD、CWD、PWVD统一到一个框架里核函数一变时频图的“摊平方式”就变了可以在交叉项抑制和能量聚集度之间主动做选择。这套MATLAB程序把WVD、CWD、PWVD封装成可切换的Cohen类分布计算工具适合做非平稳信号分析、滚动轴承故障特征提取以及给深度学习模型准备时频图像输入。无论你刚看完matlab下载安装教程还是已经在用matlab做图像处理只要把扩展名为.m的文件放进当前目录就能跑。2. Cohen类分布统一框架WVD、CWD、PWVD的核函数选型2.1 WVD为什么会有交叉项核函数等于1的代价Wigner-Ville分布是所有Cohen类分布的原点。对实信号 x(t)WVD定义在瞬时自相关 R_x(t,τ)x(tτ/2)x*(t-τ/2) 的傅里叶变换上W_x(t,f)∫_{-∞}^{∞} x(tτ/2)x^*(t-τ/2)e^{-j2πfτ}dτ由于是信号与自身延迟相乘再变换双线性结构会让两个分量在时频平面中间位置产生交叉项。交叉项是虚的干扰能量幅值往往能达到真实分量的两倍但它也是WVD高分辨率的代价。核函数为1时不做任何平滑时频聚集度最高交叉项也最明显。用这套程序时如果只看WVD结果看到中间出现一大堆条纹状振荡不要先怀疑代码有问题而是要知道这是核函数没有抑制交叉项的正常表现。2.2 CWD的核函数如何压交叉项控制σ即控制惩罚范围Choi-Williams分布的核函数是高斯指数形式φ(θ,τ)exp(-θ²τ²/σ)这里的θ是频率延迟τ是时间延迟。核函数对离原点远的(θ,τ)区域做指数衰减等效于在模糊域里过滤掉远离坐标轴的交叉项部分。σ是缩放因子σ越小核函数越集中于原点附近抑制交叉项越强但自项也会被平滑σ趋于无穷大时核函数趋向1退化成WVD。实际调试时σ在1附近开始尝试交叉项重就调小到0.1分辨率不够就调大到10。2.3 PWVD是在时间和频率两个方向分别平滑伪Wigner-Ville分布PWVD的做法是对瞬时自相关加时间窗h(τ)相当于只沿频率方向做平滑去掉瞬时自相关在时间维上的多普勒响应。所以说PWVD是“伪”的是因为它只做了一半工作交叉项只被部分抑制。想要两个方向都平滑还要再引入时间平滑窗g(s)得到平滑伪WVDSPWVD。这套程序中的tfrspaw.m如果从命名习惯看很可能是带平滑窗的伪WVD核函数实现内部通过type参数切到不同的核函数形态。下表列出三种分布在程序里对应核函数的关键差异方便你在修改type时心里有底分布类型核函数形式主要调节参数交叉项时频聚集度WVDφ1无严重最高CWDexp(-θ²τ²/σ)σ明显抑制中高PWVD时间窗h(τ)窗长部分抑制中等2.4 为什么要把核函数设计成可切换参数如果每种分布都写一遍双线性积分代码重复率高而且不同核函数之间数值误差不一致对比实验结果时很难说是核函数的差异还是积分实现的差异。把核函数抽成一个统一入口外部只传分布类型字符串和对应参数内部走同一套自相关计算、同一套离散化采样才能保证WVD、CWD、PWVD三者结果可正交对比。这也是这个项目里tfrspaw.m和tfrspbk.m分开的原因核心计算与核函数表分开后续换核函数不需要动主流程。3. 这套MATLAB程序的调用结构从demo到指定type参数3.1 文件清单里每一份是干什么的拿到压缩包后先不要急着运行把文件分一下类。tolowcase.asv是MATLAB自动保存的备份文件内容与某个函数相同可以删除不影响运行。contents.m是工具箱说明文件MATLAB在命令行执行help contents时会展示这个目录下各函数的作用。tfdemo2.m到tfdemo5.m是一组演示脚本分别展示不同信号模型下的时频分布结果。tfrview.m和tfrqview.m是可视化工具前者画二维时频图后者画时频剖面或量化指标。tfrspaw.m和tfrspbk.m是核心计算函数一个负责指定Cohen类分布的计算另一个可能是某种备用核函数实现。打开contents.m可以看到函数清单。如果某个函数名右侧的描述带有“Cohen class”“time-frequency distribution”字样那就是你要改的入口。常见做法是先运行一个demo确认路径没有问题cd /your/path/cohen_tfd which contents.m tfdemo2这段命令先切换目录再用which确认文件真的在搜索路径上最后启动演示脚本。很多matlab下载安装教程都强调过set path的重要性这里直接调用可以避免找不到函数的报错。3.2 用switch包装一个可指定type的入口函数项目描述里说“可指定Cohen类分布类型WVDCWDPWVD”实际使用时不建议每次都去改tfrspaw.m源码。我一般会写一个薄包装函数把分布类型字符串解析成不同分支转到对应的核函数调用。下面这段代码可以作为模板function [tfr, t, f] cohen_type(x, fs, type, varargin) % x: 单通道信号fs: 采样率type: WVD/CWD/PWVD N length(x); t (0:N-1) / fs; f linspace(-fs/2, fs/2, N); switch upper(type) case WVD [tfr, t, f] tfrspaw(x, t, N, 0, 0); case CWD sigma varargin{1}; % 例如 sigma 1 [tfr, t, f] tfrspaw(x, t, N, sigma, 0); case PWVD win_len varargin{1}; % 时间窗长度奇数 [tfr, t, f] tfrspaw(x, t, N, 0, win_len); otherwise error(Unknown distribution type: %s, type); end end这段代码的逻辑是先根据采样率和信号长度构造时间轴和频率轴然后用switch语句把分布类型映射到tfrspaw的第四个和第五个参数上。第四个参数这里借用为CWD的核函数缩放因子σ第五个参数借用为PWVD的时间窗长度。具体签名要以你打开的tfrspaw.m头部注释为准有的版本会直接用type,param这种参数顺序。运行时注意varargin顺序CWD必须传σPWVD必须传窗长缺参数时应该在switch分支内直接报错而不是让NaN传到傅里叶积分里。3.3 参数怎么传才不踩坑从文件名看这套程序可能有多个demo脚本但每个demo里的信号长度和采样率都不一样。传参数时最常踩的坑有三个第一个是频率轴没有做fftshift导致时频图中心在四角第二个是CWD的σ不是越大越好σ大到100以上时核函数近似为1实际得到的就是WVD第三个是PWVD的窗长必须远小于信号长度否则窗内截断效应会引入虚假振荡。下面是个更安全的调用方式x randn(1, 1024); fs 1024; sigma 1; win_len 63; [tfr_wvd, t, f] cohen_type(x, fs, WVD); [tfr_cwd, t, f] cohen_type(x, fs, CWD, sigma); [tfr_pwvd, t, f] cohen_type(x, fs, PWVD, win_len);运行后检查tfr_wvd和tfr_cwd矩阵尺寸是否一致如果不一致说明tfrspaw内部对不同type做了不同下采样。此时不要直接做矩阵差值先size看维度再决定是插值还是截断。很多matlab教程里直接imagesc(tfr)如果矩阵维度不对就会提示赋值维度不匹配。4. 实测多分量信号下比较WVD、CWD、PWVD4.1 生成带跳变调频的实验信号要对比三种分布不能用单频正弦否则交叉项显示不出来。我一般构造一个两个线性调频分量叠加的信号一个频率从50Hz扫到300Hz另一个从250Hz扫到80Hz中间有一段接近交叉。下面是生成信号的代码fs 1024; t 0:1/fs:2-1/fs; f1 50 (300-50) * t / 2; phi1 2 * pi * cumsum(f1) / fs; f2 250 (80-250) * t / 2; phi2 2 * pi * cumsum(f2) / fs; x cos(phi1) cos(phi2);这里用cumsum而不是直接对f1积分是因为瞬时频率本身就是相位的一阶导累计求和能保证相位连续。两个分量在时间中点附近频率接近会产生显著交叉项适合暴露WVD的缺点。4.2 运行三种分布并用图像对比交叉项把信号送入前面写的cohen_type函数然后分别用imagesc显示[tfr_wvd, t_ax, f_ax] cohen_type(x, fs, WVD); [tfr_cwd, ~, ~] cohen_type(x, fs, CWD, 1); [tfr_pwvd, ~, ~] cohen_type(x, fs, PWVD, 63); figure; subplot(1,3,1); imagesc(t_ax, f_ax, abs(tfr_wvd).^2); axis xy; title(WVD); subplot(1,3,2); imagesc(t_ax, f_ax, abs(tfr_cwd).^2); axis xy; title(CWD); subplot(1,3,3); imagesc(t_ax, f_ax, abs(tfr_pwvd).^2); axis xy; title(PWVD);axis xy这行很重要默认imagesc会把y轴按矩阵行方向倒置不执行的话频率从高频到低频排列容易把上升扫频看成下降扫频。运行后你会看到WVD图中两个分量中间出现一串快速振荡的干扰条纹CWD图中振荡幅度明显下降但真实分量边缘也变模糊PWVD图介于两者之间靠近时间轴方向的部分初始段有轻微拖尾。4.3 交叉项定量指标和常见误判只看图像容易被视觉骗过去。可以计算一个简单的交叉项占比对时频能量矩阵做阈值分割统计低于最大幅值10%的像素能量占比占比越高说明交叉项越重ratio_wvd cross_ratio(tfr_wvd); ratio_cwd cross_ratio(tfr_cwd); ratio_pwvd cross_ratio(tfr_pwvd); function r cross_ratio(tfr) E abs(tfr).^2; Emax max(E(:)); r sum(E(E 0.1 * Emax)) / sum(E(:)); end这个比例不是严格意义上交叉项能量因为真实分量的低幅值拖尾也会被算进去但它能快速反映能量是否被分散。对同样的信号通常ratio_wvd ratio_pwvd ratio_cwd。如果你算出来CWD的比值反而更高先检查σ是否过小或者信号本身幅值波动太大。另一种常见误判是拿analyitic信号去算三种分布但忘了对实数信号做Hilbert变换WVD对实数信号会产生负频率镜像看起来像“多了一个分量”实际上不是代码bug。使用这些.m文件时如果demo脚本里信号直接是复数就不要额外取实部。5. 用tfrview和tfrqview把核参数从肉眼调参变成定量复验运行完demo后不要满足于imagesc画出来的亮图。tfrview.m和tfrqview.m才是这个程序包里真正提升效率的部分。tfrview.m打开一个交互式时频图光标移动时会显示当前点对应的时间和频率值还能用鼠标拖拽放大局部区域适合观察交叉项振荡分布在哪一片频率区间。tfrqview.m则展示指定时间点的频率切片或指定频率的时间切片还能计算时频聚集度指标比如Rényi熵。我一般用下面这组步骤复验三种分布的参数先运行一次tfrview手动确认交叉项集中在哪个区域然后回到命令窗口用tfrqview对比不同σ下的Rényi熵。Rényi熵公式里的q阶信息量对时频能量的集中度非常敏感熵值越低说明能量越集中在少数时频点时频聚集度越高。CWD的σ从0.1到10步长0.2每个σ算一次时频矩阵和熵值收敛到最小熵的σ就是当前信号下的较优核参数。这个流程把“多试几个参数看哪张图亮”的玄学问题变成了可重复的数值选择。如果要在脚本里批量做这个操作可以直接复用计算得到的时频矩阵不需要再调用GUIsigma_list 0.1:0.2:10; renyi zeros(size(sigma_list)); for i 1:length(sigma_list) [tfr_cwd, ~, ~] cohen_type(x, fs, CWD, sigma_list(i)); p abs(tfr_cwd).^2; p p / sum(p(:)); renyi(i) -log2(sum(p.^3)); end [~, best_idx] min(renyi); best_sigma sigma_list(best_idx);这段代码计算三阶Rényi熵p.^3突出了能量高点的权重所以熵值下降意味着时频平面出现了更尖锐的峰。注意这里没有加任何平滑如果时频矩阵维度过大直接对每个σ做全矩阵计算会占用较多内存可以先把矩阵缩放到256×256再算。看到best_sigma落在0.5到4之间时说明程序运行正常如果最佳σ落在边界说明信号含噪声过重需要先做带通滤波而不是继续调核函数参数。本文还有配套的精品资源点击获取
