bt法 matlab,功率谱估计(一)— BT法与周期图法(附Mtalab代码)

本文同步发布在我的个人博客宅到没朋友,欢迎来玩。

1.前言

经典功率谱估计基于傅里叶变换的思想,典型代表为BT法和周期图法。

2.自相关函数

理论上求一个随机信号的自相关函数应该使用下面这个公式:R(s,t)=E[X(s)x(t)]R(s,t) = E[X(s)x(t)]\quadR(s,t)=E[X(s)x(t)]

但在实际应用中,我们只能得到一个随机信号有限长度的样本函数。

如果一个随机信号是均方遍历的,我们就可以用样本函数的时间自相关代替该随机信号的自相关函数,如下式:

168cb8c00510cf6b83b82c91836bb61a.png (式1)

理论上只要样本函数无限长,两者就完全相等,但实际应用中我们只能得到有限长样本函数,所以用样本函数的时间自相关去近似随机信号的自相关函数是有误差的,但只要误差在容许范围内,这就是一个很好的近似。

如果m和N都比较大,上式运算量很大,如何减小运算量,实现自相关函数的快速计算?

这就是我之前文章讲过的,用FFTFFTFFT,具体请移步还在按部就班的算信号自相关?FFT让你体验飞一般的感觉!

3.BT法进行功率谱估计

对上面所求得的r^(m)\hat{r} \left( m \right)r^(m)做傅里叶变换,即:SBT^(w)=∑−MMr^(m)ejwm∣m∣≤N−1\hat{S_{BT}} \left( w \right) = \sum_{-M}^{M} \hat{r} \left( m \right )e^{jwm}\quad \left | m \right|\leq N-1SBT​^​(w)=−M∑M​r^(m)ejwm∣m∣≤N−1

以上式为结果作为对真实功率谱S(w)S \left( w \right)S(w)的估计,BT法也称为间接法。

注:一般情况下,M

4.周期图法

周期图法又称为直接法,对式1两边直接求傅里叶变换,1式右边是两个信号卷积,时域相卷,频域相乘,所以我们可以得到周期图法表达式:SPRE^(w)=1N∣UN(w)∣2\hat{S_{PRE}} \left( w \right) = \frac{1}{N}\left | U_N \left( w \right ) \right |^{2}SPRE​^​(w)=N1​∣UN​(w)∣2

其中UN(w)=∑n=0N−1uN(n)e−jwnU_N \left( w \right ) = \sum_{n=0}^{N-1}u_N \left( n \right )e^{-jwn}UN​(w)=∑n=0N−1​uN​(n)e−jwn,这也是周期图法被叫做直接法的原因,它是直接观察数据的傅里叶变换求得的。

5.两者的联系与区别

BT法实质是在周期图法法的基础上加了一个矩形窗,即BT法是对周期图法法平滑,平滑使得BT法的方差小于周期图法,但分辨率下降。

6.MatlabMatlabMatlab仿真结果

图2 周期图法与BT法估计出的信号功率谱 (N=256,M=64)

219906ebd19e9680c294f5d0399bfaac.png

可以看到周期图法发分辨率要明显好于BT法,但同时起伏也更大。

7.MatlabaMatlabaMatlaba代码

下载代码请移步BT法和周期图法。

标签:Mtalab,周期,函数,谱估计,图法,BT,right,left

来源: https://blog.csdn.net/I_am_mengxinxin/article/details/105876548

  • 1
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
3种经典功率谱估计MATLABA代码-功率谱代码.doc 3种MATLAB的经典谱估计方 希望对大家有用~ 件所有代码: 直接: 直接又称周期图法,它是把随机序列x的N个观测数据视为一能量有限的序列,直接计算x的离散傅立叶变换,得X,然后再取其幅值的平方,并除以N,作为序列x真实功率谱的估计。 Matlab代码: clear; Fs=1000; %采样频率 n=0:1/Fs:1; %产生含有噪声的序列 xn=cos 3*cos randn); window=boxcar); %矩形窗 nfft=1024; [Pxx,f]=periodogram; %直接 plot); 改进的直接: 对于直接功率谱估计,当数据长度N太大时,谱曲线起伏加剧,若N太小,谱的分辨率又不好,因此需要改进。 1. Bartlett Bartlett平均周期图的方是将N点的有限长序列x分段求周期图再平均。 Matlab代码: clear; Fs=1000; n=0:1/Fs:1; xn=cos 3*cos randn); nfft=1024; window=boxcar); %矩形窗 noverlap=0; %数据无重叠 p=0.9; %置信概率 [Pxx,Pxxc]=psd; index=0:round; k=index*Fs/nfft; plot_Pxx=10*log10); plot_Pxxc=10*log10); figure plot; pause; figure plot; 2. Welch Welch对Bartlett进行了两方面的修正,一是选择适当的窗函数w,并再周期图计算前直接加进去,加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时,可使各段之间有重叠,这样会使方差减小。 Matlab代码: clear; Fs=1000; n=0:1/Fs:1; xn=cos 3*cos randn); nfft=1024; window=boxcar; %矩形窗 window1=hamming; %汉明窗 window2=blackman; %blackman窗 noverlap=20; %数据无重叠 range='half'; %频率间隔为[0 Fs/2],只计算一半的频率 [Pxx,f]=pwelch; [Pxx1,f]=pwelch; [Pxx2,f]=pwelch; plot_Pxx=10*log10; plot_Pxx1=10*log10; plot_Pxx2=10*log10; figure plot; pause; figure plot; pause; figure plot; 复制代码
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值