求序列的功率与能量问题

lgf402

我有一个序列xn={976，975，976，978，...}；共307个这样的元素；采样频率为1HZ。
我的目的是：想求这个序列当频率在0.003HZ在0.4HZ区间内的能量。
我参考的资料中有两个matlab写的程序，一个是求能量的程序1，另一个是求功率谱的程序2。我为了满足我能取到0.003HZ对应点的能量，也就是说频率分辨率需要达到0.001HZ；所以我把N设置为1024 （我想这里应该设置为1000，但是暂时取最接近2的整数次幂1024测试）
   

程序1：
T=1;N=1024;
X=fft(xn,N);
k=3:400;
E=T/N*sum(abs(X(k+1)).^2)


程序2：
Fs=1; %采样频率
window=boxcar(length(xn)); %矩形窗
nfft=1024;
[Pxx,f]=periodogram(xn,window,nfft,Fs) %直接法

我想问的是：(1)求0.003HZ在0.4HZ区间能量是否可以用程序1？
            (2)求序列xn的功率谱是否可以用程序2？
            (3)如果程序1和程序2是对的，那么我用程序2求出功率谱后，是否可以利用功率求出其能量。（我是这样认为的，如果有：功率=能量/序列长度，那么我求能量的话，则                有：能量=功率*序列长度。是否这样求可以？）
            (4)如果上面的(3)是成立的，为什么用(3)求出的能量与用程序1求出的能量完全不一样？数据为什么会相差好几个数量级？（我看到的现象是：程序1求出的能量是3.5e+008，而我按思路(3)中求出的话就是2.5*307=767.5）
            (5)求序列xn当频率在0.003HZ在0.4HZ区间内的能量还有别的方法吗？


希望哪位大虾指点偶一下。多谢。。。！！！！

yangzj

(1)如果T是采样频率的话,改成E=T/N*sum(abs(X(k+1)/N).^2)*2  就对了.
(2)程序2求的是功率谱密度,功率谱=功率谱密度*频率分辨率Fs/nfft;也可以:功率谱=程序1中(abs(X)/N).^2 (加矩形窗时);
(3)正确
(4).(5)按以上求法再试下

lgf402

多谢yangzj的指点:

我现在在做一个项目：用硬件以1000HZ采样频率采集信号一个5分钟长的真实人体的心电波波形图信号。再在这个得到的图中，计算出一个心跳间期序列（两个心跳之间的时间间隔的序列，比如连续两次心跳的时间间隔为976，接下来的时间间隔为975，这样子就得到一组307个元素的序列：976、977、975。。。）。我现在要对计算出的这样的一组时间序列计算其能量与功率。
    我现在有几个问题还需要请教。

（1）我对这个时间序列进行能量与功率计算时。我的采样频率该取多少？1000HZ还是1HZ？虽然我的心电图信号的采样频率为1000HZ，但是我的时间序列是在此波形图上计算出来的。两个心跳间隔时间平均为976ms，所以时间序列元素间隔的时间就平均为976ms接近1s，也就是说时间序列的元素之间的间隔接近1S，是不是时间序列元素的采样周期为1S，频率为1HZ。那么哪个频率该作为为对这个时间序列进行能量与功率分析的采样频率呢？

(2) 我的这个时间序列是有限序列（5分钟采集的信号，得到的时间序列大概就是300多个有限长的元素）。（我参考的一些资料中，好象只有无限长的序列在进行截断处理时，才会有能量泄漏，之后再加窗处理）。针对我这个序列分析能量与功率谱时，会有能量泄漏的问题吗？需要加窗处理吗？如果需要该如何处理呢？

（3）我这样进行能量与功率的分析：假设我把序列的采样频率取1HZ，为了得到0.001HZ的分辨率，我取1024个点（307个元素后补加717个0），去掉直流分量（减去平均值975）后进行FFT。
     再计算能量：sum(abs(X(k+1)/N).^2)*2，针对我的这个实际情况，这样处理还有哪些需要注意的地方或不妥的地方？
       谢谢，望指教。。。。。

yangzj

首先要清楚目的是什么再来确定怎么来做，不知道你到底是想要求什么，它的物理意义是什么？

lgf402

我的目的是想计算：
(1)这个序列当频率在0.003HZ在0.4HZ区间内的能量
(2)这个序列的功率谱.

对应的物理意义:
(1) 0.003-0.4HZ范围内的能量在0.003-0.04\0.04-0.15\0.15-0.4这三个子区间的分布对应了不同的生理现象.
(2)用功率谱来分析生理变化的趋势

yangzj

那"计算出一个心跳间期序列（两个心跳之间的时间间隔的序列，比如连续两次心跳的时间间隔为976，接下来的时间间隔为975，这样子就得到一组307个元素的序列：976、977、975。。。）"这一步是用来做什么的?

lgf402

我这句话只是用来描叙我的这个时间序列是怎么得到的(已经实现)。
也就是说这一步的作用就是得到一个时间序列,
我的目的正是想计算这样一个时间序列的：
(1)这个序列当频率在0.003HZ在0.4HZ区间内的能量
(2)这个序列的功率谱.

yangzj

这样做有根据吗?直接在采集的原始序列上进行计算才对吧

lgf402

有医学意义的。
思路是这样的，先通过采集器采集心电信号。把采集的心电信号画出波形图，在波形图上识别出峰值特征点，再把波形的峰值间隔求出来，这个间隔医学上称为RR间期（心跳间隔），得到的这个RR间期序列后，这个序列值的变化反映了心率变异性，医学上对心率变异性的分析得到一些生理现象。
我现在要做的就是对这个序列的频域特性进行分析。也就是我希望得到这个序列的能量与功率谱。

lgf402

希望大家能再指点我一下。。
  多谢。

lgf402

麻烦帮我再分析下吧。。。谢谢

吃书的老虎

你的数据我感觉蛮有意思的，首先你的（xn={976，975，976，978，...}；）时间间隔是不相等的；你取0.003～0.4HZ的信息，取样的时间间隔接近1S，这对取0.4HZ附近的信息不准确啊，因为工程一般采样率大于信号最大频率的5倍啊；我想你可以取原始信号的包络线，在对包络线分析可以吗？或把你的原始信号发到网上，帮你看看。

yangzj

估计对所谓的心率变异性进行频谱分析时对应的频率物理意义以经不同于传统的物理意义了.传统频谱分析分析的序列是对应时间的物理量,而你这个序列其实是第n个心跳间隔,取N个点进行频谱分析得到的频谱图第k条谱线的意义应该是:每N/k次心跳间隔按余弦波动的幅度.
所以如果把频率的意义扩展下的话,你可以把提取后的序列的采样间隔看成是1(单位其实不是秒,而是次),取N点进行频谱分析和频率分辨率为1/N(这时的频率单位Hz其实不是1/秒,而是1/次).

lgf402

根据你的提示， 我对这个序列我又进行了重采样，比如我用2HZ的采样率对这个时间序列进行重采样，这样他们的采样频率就是一致的，对应的幅度值我就在原来时间序列的基础上用插值的方法求出来，得到一个新的时间序列，这样就和传统的频谱分析分析的序列物理意义是一致的了。
    （1） 这样子我就可以用传统的频谱分析方法对重采样后的序列进行频谱分析了。这样子操作可以吗？
     （2）我的这个重采样后的时间序列是有限序列（5分钟采集的信号，得到的时间序列大概就是600多个有限长的元素）。（我参考的一些资料中，好象只有无限长的序列在进行截断处理时，才会有能量泄漏，之后再加窗处理）。针对我这个重采样后的序列分析能量与功率谱时，会有能量泄漏的问题吗？需要加窗处理吗？如果需要该如何处理呢？
      谢谢指点。。。

lgf402

另；我进行FFT的长度假设为1024个元素。（这样时间序列的600个元素就没有截断处理，后面补加424个0。）