c语言编译功率谱密度函数,科学网—6、功率谱密度函数估计 - 柏世平的博文
功率谱密度函数估计,在随机信号处理中具有极其重要的意义。不管是为了目前增加对信号基本属性的了解,还是为了以后对信号作进一步分析处理,现在对敝人各生理信号作一下功率谱估计,都是很有必要的。
MATLAB信号处理工具sptool中,有8种已经很成熟的功率谱估计方法(FFT法,Welch法,Covariance法,Mod.Covariance法,MTM法,MUSIC法,Burg法,Yule AR法)可供调用,非常方便。
将体温信号序列TW461512_0导入sptool,并Create(加载)于spectra栏,即打开功率谱处理、观察窗口Spectrum Viewer。
(1)、先看看快速傅立叶变换FFT方法,也就是所谓的经典周期图法吧。取FFT执行点数Nfft=2^16=65536。在2的n次幂中,它比TW461512_0的长度55380长,且最接近55380。采样频率Fs取Fs=Nfft=65536(次/单位采样时间,我亦称之为“Hz”,将“1次/单位采样时间”的采样频率提高Nfft倍,是为了使频率轴都用正整数表示。
图6-1 TW461512_0快速傅立叶FFT法估计的功率谱图
此图中左标尺所标示的谱线位置,为频率f=5461Hz处。其波形周期T=Fs/f=65536Hz/5461Hz=12.0007单位采样时间,即1天。其小数点后面的误差为数字化误差。此谱线右边的谱线为其倍频谱线。
图6-2 FFT法功率谱最低频处。频率轴放大128倍,幅值轴放大4倍。
波峰狭窄尖锐,幅值相近,噪声谱与信号谱混杂在一起,很难分辨样本信号谱峰。
图6-3 FFT法功率谱最高频处。放大倍数同前。
谱峰形状与最低频处一样。
(2)、来看看使用Welch(经典的韦尔奇)方法估计功率谱。此方法需要选择的参数比较多。除了FFT执行点数Nfft外,还需要确定窗函数Window及其长度Nwind、分段重叠点数Overlap。对于没有足够的先验知识的人来说,要一次选对各个参数是很难的。先比较随意地选定一组参数:Window=hanning,Nwind=4096,Overlap=1024,作一功率谱图如下:
图6-4 TW461512_0经典welch法估计的功率谱图
图6-5 图6-4最低频处。横轴放大16倍,纵轴放大1.2倍。
图6-6 图6-4最高频处。横轴放大16倍,纵轴无放大。
直觉上感到此谱图各波形太过平缓,波峰幅值相近,难以确定信号波峰及其频率点。现在以窗函数宽度及其重叠的比例为变量,作搜索式计算。编程如下:
%welch法搜索参数
L=length(TW461512_0);%=55380;
m0=512;%窗函数最小宽度
c=floor(L/(2*m0));%信号序列最大分段数
for i=1:4
M=zeros(300,c);
for j=1:c
m=j*m0;%窗宽按每次增加一个最小窗宽数m0搜索
n=m*(i-1)/4;%数据分段重叠比例
Pxx=PWELCH(TW461512_0,hanning(m),n,65536,65536);
M(1:149,j)=Pxx(1:149);%数据太长,无法全部清晰作图显示,只能
%取最低频与最高频段各149点
M(150:151,j)=ones(2,1);%人工隔板
M(152:300,j)=Pxx(length(Pxx)-148:length(Pxx));
end
M=M';%横坐标表示频率,纵坐标表示窗宽
figure(i)
surf(log10(M))
end
运行,得:
图6-7 TW461512_0的Welch法估计功率谱。分段无重叠
图6-8 TW461512_0的Welch法估计功率谱。分段重叠比例为1/4
图6-9 TW461512_0的Welch法估计功率谱。分段重叠比例为2/4
图6-10 TW461512_0的Welch法估计功率谱。分段重叠比例为3/4
由于数据太长,无法在一张图片中显示整个结果,只能分别取最低频最高频段各149Hz合在一起显示。中间的“隔板”是另加的。
上面四张图片,分别是信号序列分段无重叠、重叠1/4、2/4、3/4时的谱图。其纵坐标表示窗函数的宽度(以512点为一个单位)。我认为谱估计效果最好的应该是最后一张图,即信号序列分段重叠3/4的时候。因为它在整个纵轴方向上,随着窗宽增加,波峰略为变高变窄,而位置(频率点)都很少改变,而不像前面三张图,各波形波峰高度及其频率点都随着窗函数宽度的改变而有所改变,尤其是窗函数宽度比较小的时候。
我认为一种好的估计方法,当需要选择的参数落在一个大范围后,其估计结果应该对参数的变化是“不敏感”的。原因很简单:真理只有一个。如果反之,估计结果对参数的变化是“敏感”的,参数稍微一变化,估计结果就跟着很快变化,那么我们最终到底要选哪一个参数?如果随便选一个,那岂不就象买彩票一样?科学要追求确定性的结果,而不是中彩。而这里的参数必定是由人来选定的。如果参数可以由一套算法来计算、确定,那就不用设参数,而设置一个函数好了。
根据上面的分析,本信号序列采用Welch方法估计功率谱时,窗宽重叠率须取1/2以上。至于窗宽,取一个中间偏大的值如512×(35~54)都可(第一行的窗宽为512点,最后一行窗宽为512×54=27648点)。谱图图片暂就不贴了,等到在同一窗口比较各方法谱估计曲线时一并贴出。
(3)Covariance(协方差)法。此法需选择参数,除了快速傅立叶变换执行点数Nfft,就是阶数Order了。取Nfft=65536,任取几个阶数Order运行,发现此法非常耗用计算机资源。当Order>500时即Out of memory(超出内存)了。
我的内存为2×1G。要增加内存,牵涉到整机升级的问题,不是换两根内存条就可以了的。据说可以设置虚拟内存,但我设置之后,也没见运算能力的明显提高。不知是何缘故。
自己编程,看看Order<=500时的谱估计情况:
%Covariance法搜索阶数的计算
M=zeros(300,10);
for i=1:10
order=50*i;
Pxx=pcov(TW461512_0,order,65536,65536);
M(:,i)=Pxx(1:300);
end
M=M';
surf(log10(M))
运行,得:
图6-11 Covariance法搜索阶数的计算
什么也没有。所以在个人电脑上是无法使用Covariance(协方差)法的。
(4)、Mod.Covariance(改进的协方差)法。在Spectrum Viewer窗口中
略一运行,便知此法比Covariance法更耗电脑资源。pass。
(5)、MTM(Multitaper多窗口)法。此方法需要选择的参数有时间-带宽乘积nw、fft执行点数Nfft、采样频率Fs、权重算法Weights。
取Nfft=Fs=65536,Weights='adapt'(Thomson的非线性自适应组合算法),各取nw=2、3、4,得到3条谱估计曲线。使3条曲线在同一窗口显示,见图6-12。
图6-12 时间-带宽乘积nw不同,其余参数相同时,TW461512_0的功率谱估计曲线。
其中绿线nw=2,红线nw=3,蓝线nw=4。
图6-13 图6-12最低频段处。横轴放大128倍,纵轴放大(为原来)2倍。
图6-14 图6-12最高频段处。横轴放大128倍,纵轴放大(为原来)2倍。
图中3条曲线在大的走势上虽然保持了一致,但在更细节的形状上差异也不小,而且很多波峰形状也不明显,幅值接近,使人很难判断真实的谱线。
此方法中,2*nw-1表示采用的窗口数。nw的典型取值有2,5/2,3,7/2,4。但我发现在上述取值范围内,谱估计函数Pxx=pmtm(x,nw,Nfft,Fs,Weights)对所有实数nw都是有意义的。仿照我前面搜索welch法参数的方法,作搜索nw的计算。程序如下:
%搜索MTM谱估计法中时间-带宽乘积参数nw
Nfft=65536;fft执行点数
Fs=65536;%采样频率
M=zeros(300,30);
for i=16:45
nw=0.1*i;%时间-带宽乘积
Pxx=pmtm(TW461512_0,nw,Nfft,Fs,'adapt');
M(1:149,i)=Pxx(1:149);%最低频段149HZ
M(150:151,i)=ones(2,1);%人造隔板
M(152:300,i)=Pxx(length(Pxx)-148:length(Pxx));%最高频段149HZ
end
M=M';
surf(log10(M))
运行,得:
图6-15 TW461512_0的pmtm法估计功率谱。权重算法为adapt
觉得还是很难确定nw到底取多大为好。改程序中权重算法因子为unity(一致权重的线性组合算法)、eigen(特征值权重的线性组合算法),运行,得:
图6-16 TW461512_0的pmtm法估计功率谱。权重算法为unity。
图6-17 TW461512_0的pmtm法估计功率谱。权重算法为eigen。
发现上面三张图中,曲面形状非常一致。注:上面三张图及图6-7~图6-10的纵坐标应乘10才是分贝值,因为程序中“surf(log10(M))”应为“surf(10*log10(M))”。懒得再截图了,说明一下。
回到Spectrum Viewer窗口,取定nw=3,Nfft也不变,分别取Weights='adapt'、'unity'、'eigen',得到3条曲线。令其在同一窗口显示,如下图:
图6-18 Weights不同,其余参数相同时,TW461512_0的功率谱估计曲线。最低频段处,横轴放大128倍,纵轴放大为2倍。其中绿线Weights='eigen',红线Weights='unity',蓝线Weights='adapt'
可以看出,最低频段处,3条曲线基本上重合了,后面的绿线、蓝线基本上都被前面的红线挡住。
图6-19 TW461512_0的功率谱估计曲线。最高频段处。其余情况同图6-18
最高频段处,绿线与红线曲线基本上重合了,蓝线在某些部位不与红线、绿线重合。
所以,MTM法中参数Weights的选择是不敏感的,选哪一个都差不多。
(6)、MUSIC(多信号分类)法。该法需要选择的参数比较多,有信号子空间数Signal Dim.、阈值、fft长度Nfft、窗函数、窗函数宽度Nwind、窗函数重叠点数Overlap。随选几组参数,在Spectrum Viewer窗口试运行,即知此方法非常占用电脑资源。信号子空间Signal Dim.参数达到4000、加窗宽度Nwind达到7680(=512*15)时便“Out of memory”了。
现在我觉得,Matlab中Sptool的Spectrum Viewer窗口功率谱估计,最适合的是那些做固定项目工程应用的人。他们对谱估计相关的参数选择都已经很熟悉,只是由于样本数据的变化,而在此窗口查看一下谱估计的结果。对于做科研工作的人来说,大多数是第一次接触到新项目,谱估计参数选择不熟悉。此窗口只能作辅助性的研究工具。
准备看看在Signal Dim.<=4000、Nwind<=7680时的功率谱估计情况。先固定窗函数及其宽度、窗函数重叠点数,对Signal Dim.作搜索计算。程序在傍晚时开始运行,结果运行到第二天临晨还没有结束,只好强行中断。虽然感觉再要不了多久就会结束,但也不是十分有把握。如果还需要运行几十、几百……个小时,那不麻烦了。第二天先试循环少量几次,对整个程序的运行时间有一个大概的评估后才开始正式的运行。所以对有循环语句的程序,这个工作是少不了的,除非非常有经验。
将原程序略作改动,分段运行:
pack%整理工作空间
Nfft=65536;fft执行点数
Fs=65536;%采样频率
d=100;%每段程序循环搜索次数
k=20;%信号子空间数每次增加数
Dim0=2000;
M_2000_20_4000=zeros(300,d);
m=5120;%窗函数宽度
n=m/2;%数据分段重叠比例
for i=1:d
Dim=Dim0+i*k;
Pxx=pmusic(TW461512_0,Dim,Nfft,Fs,hanning(m),n);
M_2000_20_4000(1:150,i)=Pxx(1:150);
M_2000_20_4000(151:300,i)=Pxx(length(Pxx)-149:length(Pxx));
end
M_2000_20_4000=M_2000_20_4000';
运行完毕后,将其中d=100; k=20;Dim0=2000;分别改成d=100、100、75;
k=10、7、4;Dim0=1000、300、0;M_2000_20_4000改成M_1000_10_2000、M_300_7_1000、M_0_4_300,分4段运行。4段程序大约运行了15、6个小时。
最后:
M=[M_0_4_300;M_300_7_1000;M_1000_10_2000;M_2000_20_4000];
surf(10*log10(M))
运行,得:
图6-20 TW461512 _ 0的pmusic法估计功率谱,对信号子空间Signal Dim.的搜索。hanning窗宽5120点,重叠2560点
可以看出,Signal Dim.太小的时候,基本上没有什么谱峰。即便在Signal Dim.最大值4000处,谱峰也很平缓。这是低频段的情况。在高频段则始终没有出现谱峰。
总而言之,本人现在的电脑配置是没法使用MUSIC方法估计本样本数据的功率谱的。最后再贴几张不同参数下的功率谱估计图吧:
图6-21 Dim=4000,Nwind=15*512=7680
图6-22 Dim=4000,Nwind=7*512=3584
图6-23 Dim=3200,Nwind=2*512=1024
图6-24 Dim=1000,Nwind=2*512=1024
(7)、Burg法。前面6种估计方法都是属于非参数估计方法,Burg法及最后的Yule AR法属于参数估计方法。在Spectrum Viewer窗口中,Burg法需要选择的参数,除了fft长度Nfft,就只有阶数Order。试选择几个Order数运行,知道其可运行范围很大(Order=20000都可以)。反正我对于Order的具体选择也没有经验,那就编程作搜索式运算吧。
%用pburg法估计TW461512_0功率谱,搜索阶数参数Order
Nfft=65536;
M=zeros(Nfft/2+1,240);
for i=1:240
Order=50*i;
Pxx=pburg(TW461512_0,Order,Nfft);
M(:,i)=Pxx;
end
P_TW4615_b50_50_12000=M;%保存
M(151:(Nfft/2+1)-150,:)=[];
M=M';
surf(10*log10(M))
运行,得:
图6-25 用pburg法估计TW461512_0功率谱,搜索阶数参数Order
上图数据太密集,图形太暗。
M(1:2:240,:)=[];%删去一部分数据
surf(10*log10(M))
图6-25(1) 用pburg法估计TW461512_0功率谱,搜索阶数参数Order
可以看出,低频段,阶数达到8、9000以上,高频段,阶数达到5、6000以上后,谱峰形状与位置都很稳定。这也符合我在用welch方法估计功率谱时所说的“我认为一种好的估计方法,当需要选择的参数落在一个大范围后,其估计结果应该对参数的变化是“不敏感”的”的特点。
贴几张Spectrum Viewer窗口谱图:
图6-26 用pburg法估计TW461512_0功率谱,阶数Order=10000
图6-27 图6-26最低频段处,横轴放大128倍,纵轴放大为2倍
图6-28 图6-26最高频段处,横轴放大128倍,纵轴放大为2倍
从波形上来看,总的感觉,参数估计法比非参数估计法效果就是要好得多。
(8)、Yule AR法。此法与Burg法极其相似,在Spectrum Viewer窗口中需选择的参数完全一样。搜索阶数Order的程序只需将其谱估计函数“pburg”换成“pyulear”就可以了。
%用pyulear法估计TW461512 _ 0功率谱,搜索阶数参数Order。
Nfft=65536;
M=zeros(Nfft/2+1,120);
for i=1:120
Order=100*i;
Pxx=pyulear(TW461512_0,Order,Nfft);
M(:,i)=Pxx;
end
P_TW4615_y100_100_12000=M;%保存
M(151:(Nfft/2+1)-150,:)=[];
M=M';
surf(10*log10(M))
运行,得:
图6-29 用pyulear法估计TW461512_0功率谱,搜索阶数参数Order
第一感觉,就是此图与图6-25也极其相似。再贴几张Spectrum Viewer窗口谱图:
图6-30 用pyulear法估计TW461512_0功率谱,阶数Order=10000
图6-31 图6-30最低频段处,横轴放大128倍,纵轴放大为2倍
图6-32 图6-30最高频段处,横轴放大128倍,纵轴放大为2倍
感觉上面三图与图6-26、图6-27及图6-28也极其相似。将图6-27与图6-31、图6-28与图6-32中的曲线分别放在同一窗口中进行比较。见下二图:
图6-33 最低频段处。蓝线是Burg法、红线是Yule AR法、绿线是Welch法
图6-34 最高频段处。其余情况同上图
绿线Welch法谱估计参数是hanning窗,窗宽Nwind=25600,窗函数重叠点数Overlap=19200。
可以看到,蓝线与红线重叠程度非常高。绿线作为一种非参数估计方法,其曲线形状变化步调与蓝线、红线能够保持这么高,已经很不简单了。
上面就目前最常用的8种谱估计方法,对体温信号序列TW461512_0作了初步的功率谱估计。对这些谱估计的结果作进一步的分析处理,就放到以后的篇目中再做。
c语言编译功率谱密度函数,科学网—6、功率谱密度函数估计 - 柏世平的博文相关推荐
- 计算机科学家的摇篮,科学网—科学家有多少摇篮? - 籍利平的博文
在一则题为<懒得考清华北大?16岁女生考上中国版"麻省理工">的微信公众号文章后面,我曾经这样留言(个别笔误已经更正): 各有所好吧,不必因此贬低北大清华.上少年班,的 ...
- vc2010解决方案项目编译顺序_科学网—VS2012 (2008,2010) 编译问题解决合集 - 冯博远的博文...
问题一: VS2012 (包括从VS2008,VS2010) 出现编译错误:LINK : fatal error LNK1104: cannot open file 'LIBC.lib' 的解决办法: ...
- R语言导出为html,科学网—[转载]R语言中数据的导入与导出(笔记) - 刘朋的博文...
!!!help文档!!! 1.导入数据语句为mydata 要分析的.csv数据(.xlsx另存为...)导入.数据导入后可以edit(mydata),R语言工作区就会弹出数据, 可以进行编辑和修改.还 ...
- r语言pls分析_科学网—R语言统计:偏最小二乘路径模型(plspm) - 涂波的博文...
R包"plspm" 作者:Gaston Sanchez 单位:Berkeley, California. 包使用说明文件:http://www.gastonsanchez.com/ ...
- c语言酶切算法,科学网—FitHiC V1算法解析(一) - 卢锐的博文
FitHiC V1主要用于识别中程顺式互作 处理Hi-C数据,最自然的分辨率划分方法是基于限制性内切酶切出来的酶切片段,即一个酶切片段为一个最小单位.但是,因为测序深度和基因组上感兴起的size不同, ...
- docker安装gamit_科学网—Ubuntu系统GAMIT/GLOBK程序安装 - 陈超的博文
最近开始学习GAMIT,网上资料还是蛮多的,但是感觉都是东拼西凑的,一点都不系统,一点不适合初学者.安装教程也是乱七八糟的,下面把我的安装过程分享一下:(我也是参考的网上一个教程,很久以前下载的,地址 ...
- python做社会网络分析_科学网-python 社会网络分析工具之igraph-郗强的博文
1.networkx 2.igraph 3.SNAP 2.igraph igraph是免费的复杂网络(graphs)处理包,可以处理百万级节点的网络(取决于机器内存).igraph提供了R和C语言程序 ...
- 信息与计算机,科学网—信息与计算机(1) - 姜咏江的博文
信息与计算机(1)信息的概念姜咏江 当今科学研究中使用最多的概念是什么?信息!然而信息却是那样地让人们困惑.有句话说:"把复杂的问题说简单了,那是学问:将简单的问题说复杂了,那是蒙人!&qu ...
- html哭脸字符,科学网—Unicode中的符号 - 丁祥欢的博文
可是,这句话基本上是正确的,只要你提到的符号不是组合符号或者是你自己临时创造的符号,你都可以在以下文件中找到,这个文件是就是Unicode 13.0的字符集. 提取码: vis3 当然,你也可以在线查 ...
- python 网络_科学网-python 社会网络分析工具之networkx-郗强的博文
1.networkx 2.igraph 3.SNAP 1.networkx NetworkX是一个用Python语言开发的图论与复杂网络建模工具,内置了常用的图与复杂网络分析算法,可以方便的进行复杂网 ...
最新文章
- nginx中的数组结构ngx_array_t
- 10 个优质的 Laravel 扩展推荐
- 新媒体中的MCN机构是什么意思
- java企业人事管理系统源码_企业人事管理系统完美版源代码 - 源码下载|行业应用软件|企业管理(财务/ERP/EIP等)|源代码 - 源码中国...
- 零基础小白10分钟用Python搭建小说网站!网友:我可以!
- Matlab--m代码转C与C++代码)1(简单示例涉及到函数调用)
- 传网易云音乐高管变动:市场副总裁李茵离职 CEO被降权
- tpc-c 服务器性能,IBM创英特尔8处理器服务器TPC-C性能记录
- 利用skipList(跳表)来实现排序(待补充)
- 自然语言处理----词干提取器
- L298N电机驱动模块详解
- 使用Git的Kdiff3解决合并冲突 显示乱码的问题
- 2021-07-19支付宝扫码点餐推广怎么做(干货来了)
- JVM-JConsole:Java监视与管理控制台(windows)
- 「津津乐道播客」#301 这是一期价值3000元的当代社畜科学点餐指南
- 锋利的jQuery读书笔记-第1章 认识jQuery
- choerodon-ui/pro入门 - dataset 的使用
- OEPNCV 轮廓提取函数findContours中所用的算法原理疑问。
- 在Mac电脑的状态栏右上角显示自己的名字
- 设计一个用于人事管理的“人员”类