虚拟商城网站,电商运营工作很难做吗,美食网站,个人网站开发软件经典功率谱估计的原理及MATLAB仿真#xff08;自相关函数BT法、周期图法、bartlett法、welch法#xff09; 文章目录 前言一、BT法二、周期图法三、Bartlett法四、welch法五、MATLAB仿真六、MATLAB详细代码总结 前言
经典功率谱估计方法包括BT法#xff08;对自相关函数求傅… 经典功率谱估计的原理及MATLAB仿真自相关函数BT法、周期图法、bartlett法、welch法 文章目录 前言一、BT法二、周期图法三、Bartlett法四、welch法五、MATLAB仿真六、MATLAB详细代码总结 前言
经典功率谱估计方法包括BT法对自相关函数求傅里叶变换求功率谱、周期图法、Bartlett法分段求平均、welch法有重合分段求平均。本文在总结各种方法的原理后将在MATLAB平台上进行仿真完成对功率谱密度的估计。 提示以下是本篇文章正文内容转载请附上链接
一、BT法
设 u N ( 0 ) , u N ( 1 ) , ⋅ ⋅ ⋅ , u N ( N − 1 ) u_N(0),u_N(1),\cdotp\cdotp\cdotp,u_N(N-1) uN(0),uN(1),⋅⋅⋅,uN(N−1) 为广义平稳随机过程 u ( n ) u(n) u(n) 的 N N N 个观测值为描述方便本章观测样本的序号为从0到 N − 1 N-1 N−1的整数)且设 u N ( n ) u_N(n) uN(n)其他时刻的值为零则 u N ( n ) u_N(n) uN(n)可表示为 u N ( n ) { u ( n ) , 0 ⩽ n ⩽ N − 1 0 , 其他 (1) u_N(n)\begin{cases}u(n),0\leqslant n\leqslant N-1\\\\0,\text{其他}\end{cases} \tag{1} uN(n)⎩ ⎨ ⎧u(n),0,0⩽n⩽N−1其他(1) u ( n ) u(n) u(n) 的自相关函数 r ( m ) r(m) r(m) 可用时间平均进行估计即 r ^ ( m ) 1 N ∑ n 0 N − 1 u N ( n ) u N ∗ ( n − m ) , ∣ m ∣ ⩽ N − 1 (2) \hat{r}(m)\frac{1}{N}\sum_{n0}^{N-1}u_N(n)u_N^*(n-m),\quad\mid m\mid\leqslant N-1 \tag{2} r^(m)N1n0∑N−1uN(n)uN∗(n−m),∣m∣⩽N−1(2) 根据维纳-辛钦定理对由式(2) 估计得到的自相关函数 r ^ ( m ) \hat{r}(m) r^(m)求傅里叶变换可得功率谱的估计为 S ^ B T ( ω ) ∑ m − ∞ ∞ r ^ ( m ) e − j ω m ∑ m − N 1 N − 1 r ^ ( m ) e − j ω m (3) \hat{S}_{\mathrm{BT}}\left(\omega\right)\sum_{m-\infty}^{\infty}\hat{r}\left(m\right)\mathrm{e}^{-\mathrm{j}\omega m}\sum_{m-N1}^{N-1}\hat{r}\left(m\right)\mathrm{e}^{-\mathrm{j}\omega m} \tag{3} S^BT(ω)m−∞∑∞r^(m)e−jωmm−N1∑N−1r^(m)e−jωm(3)
考虑到自相关函数 r ^ ( m ) \hat{r}(m) r^(m) 在 ∣ m ∣ N − 1 \mid m\midN-1 ∣m∣N−1 时为零且 r ^ ( m ) \hat{r}(m) r^(m) 在 m m m 接近 N − 1 N-1 N−1 时性能较差 式(3)经常表示为 S ^ B T ( ω ) ∑ m − M M r ^ ( m ) e − j ω m , 0 ⩽ M ⩽ N − 1 (4) \hat{S}_{\mathrm{BT}}\left(\omega\right)\sum_{m-M}^M\hat{r}\left(m\right)\mathrm{e}^{-\mathrm{j}\omega m}\:,\quad0\leqslant M\leqslant N-1 \tag{4} S^BT(ω)m−M∑Mr^(m)e−jωm,0⩽M⩽N−1(4)
以此结果作为对理论功率谱 S ( ω ) S(\omega) S(ω) 的估计因为这种方法估计出的功率谱是通过自相关函数间接得到的所以此方法又称为间接法。
当时延m≥0时r(m)的均值可以表示为 E { r ^ ( m ) } E { 1 N ∑ n 0 N − 1 u N ( n ) u N ∗ ( n − m ) } 1 N ∑ n m N − 1 E { u N ( n ) u N ∗ ( n − m ) } 1 N ∑ n m N − 1 r ( m ) N − m N r ( m ) (5) \begin{aligned} \operatorname{E}\{\hat{r}(m)\} \text{ E}\left\{\frac1N\sum_{n0}^{N-1}u_N\left(n\right)u_N^*\left(n-m\right)\right\} \\ \frac1N\sum_{nm}^{N-1}E\{u_N(n)u_N^*(n-m)\} \\ \frac1N\sum_{nm}^{N-1}r\left(m\right) \\ \frac{N-m}Nr\left(m\right)\tag{5} \end{aligned} E{r^(m)} E{N1n0∑N−1uN(n)uN∗(n−m)}N1nm∑N−1E{uN(n)uN∗(n−m)}N1nm∑N−1r(m)NN−mr(m)(5)
考虑到 m m m 的取值可正可负所以有 E { r ^ ( m ) } N − ∣ m ∣ N r ( m ) (6) \operatorname{E}\{\hat{r}(m)\}\frac{N-\mid m\mid}{N}r(m) \tag{6} E{r^(m)}NN−∣m∣r(m)(6)
BT法结论 r ^ ( m ) \hat{r}(m) r^(m) 是 r ( m ) r(m) r(m) 的渐近无偏估计。 r ^ ( m ) \hat{r}(m) r^(m) 是 r ( m ) r(m) r(m) 的渐近一致估计。
二、周期图法
由于随机过程 u ( n ) u(n) u(n) 的 N N N个观测值 u N ( n ) u_N(n) uN(n) 是确定信号对其进行傅里叶变换得 U N ( ω ) ∑ n 0 N − 1 u N ( n ) e − j ω n (7) U_N(\omega)\sum_{n0}^{N-1}u_N(n)\:\mathrm{e}^{-\mathrm{j}\omega n}\tag{7} UN(ω)n0∑N−1uN(n)e−jωn(7)
根据帕塞瓦尔定理上式的模的平方是确定信号 u N ( n ) u_N(n) uN(n)的能量谱对能量谱除以持续时间 N N N, 其结果应是 u N ( n ) u_N(n) uN(n) 的功率谱估计将其作为随机过程 u ( n ) u(n) u(n)的功率谱的估计表示为 S ^ P E R ( ω ) 1 N ∣ U N ( ω ) ∣ 2 (8) \hat{S}_{\mathrm{PER}}(\omega)\frac1N\mid U_N(\omega)\mid^2 \tag{8} S^PER(ω)N1∣UN(ω)∣2(8)
该方法称为功率谱估计的周期图法periodogram)。因为这种功率谱估计方法是直接通过观测数据的傅里叶变换求得的称之为直接法。
周期图法结论 r ^ ( m ) \hat{r}(m) r^(m) 是 r ( m ) r(m) r(m) 的渐近无偏估计。 r ^ ( m ) \hat{r}(m) r^(m) 不是 r ( m ) r(m) r(m) 的一致估计估计方差趋近于一个有限值。
三、Bartlett法
Bartlett法的基本步骤是将 N N N 点的观测数据 u N ( n ) u_N(n) uN(n) 分为 L L L 段每段数据的长度为 M M M, 整数 L L L 和 M M M 满足 L N M (9) L\frac NM \tag{9} LMN(9)
第 i i i 段数据加矩形窗后有 u N i ( n ) u ( n ( i − 1 ) M ) w M ( R ) ( n ) , 0 ⩽ n ⩽ M − 1 , 1 ⩽ i ⩽ L (10) u_N^i(n)u(n(i-1)M)w_M^{(R)}(n)\:,\quad0\leqslant n\leqslant M-1,1\leqslant i\leqslant L \tag{10} uNi(n)u(n(i−1)M)wM(R)(n),0⩽n⩽M−1,1⩽i⩽L(10) 其中 , w M ( R ) ( n ) ,w_{M}^{(R)}(n) ,wM(R)(n) 是长度为 M M M 的矩形窗。
对于每段数据 u N i ( n ) u_{\mathbb{N}}^{i}\left(n\right) uNi(n), 先利用周期图法求得其功率谱 S ^ P E R i ( ω ) \hat{S}_{\mathrm{PER}}^i(\omega) S^PERi(ω), 限 S ^ P E R i ( ω ) 1 M ∣ ∑ n 0 M − 1 u N i ( n ) e − j ω n ∣ 2 , 1 ⩽ i ⩽ L (11) \hat{S}_{\mathrm{PER}}^i(\omega)\frac{1}{M}\Big|\sum_{n0}^{M-1}u_N^i(n)\:\mathrm{e}^{-\mathrm{j}\omega n}\Big|^2,\quad1\leqslant i\leqslant L\tag{11} S^PERi(ω)M1 n0∑M−1uNi(n)e−jωn 2,1⩽i⩽L(11)
然后计算各段功率谱估计的平均得到平均周期图 S ‾ P E R ( ω ) \overline{\mathcal{S}}_{\mathrm{PER}}\left(\omega\right) SPER(ω), 即 S ‾ P E R ( ω ) 1 L ∑ i 1 L S ^ P E R i ( ω ) 1 L M ∑ i 1 L ∣ ∑ n 0 M − 1 u N i ( n ) e − j ω n ∣ 2 (12) \overline{S}_{\mathrm{PER}}\left(\omega\right)\frac{1}{L}\sum_{i1}^{L}\hat{S}_{\mathrm{PER}}^{i}\left(\omega\right)\frac{1}{LM}\sum_{i1}^{L}\Big|\sum_{n0}^{M-1}u_{N}^{i}\left(n\right)\mathrm{e}^{-\mathrm{j}\omega n}\Big|^{2}\tag{12} SPER(ω)L1i1∑LS^PERi(ω)LM1i1∑L n0∑M−1uNi(n)e−jωn 2(12)
Bartlett法结论 Bartlett功率谱估计使得数据点数从 N N N 减小为 M N / L MN/L MN/L, 于是窗函数的频谱宽度增大为周期图法的 L L L 倍因此频率分辨率下降为原来的 1 / L 1/L 1/L 。 Bartlett功率谱估计的方差减小为周期图法的 1 / L 1/L 1/L, 因此Bartlett功率谱估计较周期图法的结果更为平滑。
四、welch法
Welch法也对 N N N点的信号 u N ( n ) u_N(n) uN(n) 进行分段但允许分段时每段信号样本重叠。若每段的样本长度为 M M M, 信号被分为 L L L 段取相邻两段的信号样本重叠 50 % 50\% 50% 则 L L L满足 L N − M / 2 M / 2 (13) L\frac{N-M/2}{M/2}\tag{13} LM/2N−M/2(13)
将每段信号 u N i ( n ) u_N^i(n) uNi(n) 和窗函数 w ( n ) w(n) w(n)可以采用矩形窗、三角窗、汉宁窗或海明窗等相乘然后按式(11) 得到每段信号的功率谱估计为 S ^ P E R i ( ω ) 1 M ∣ ∑ n 0 M − 1 u N i ( n ) w ( n ) e − j ω n ∣ 2 (14) \hat{S}_{\mathrm{PER}}^i(\omega)\frac{1}{M}\Big|\sum_{n0}^{M-1}u_N^i(n)w(n)\:\mathrm{e}^{-\mathrm{j}\omega n}\Big|^2\tag{14} S^PERi(ω)M1 n0∑M−1uNi(n)w(n)e−jωn 2(14)
由此得到修正的周期图 S ‾ P E R ( ω ) \overline{\mathcal{S}}_{\mathrm{PER}}\left(\omega\right) SPER(ω) 为 S ‾ P E R ( ω ) 1 L ∑ i 1 L S ^ P E R i ( ω ) 1 L M ∑ i 1 L ∣ ∑ n 0 M − 1 u N i ( n ) w ( n ) e − j ω n ∣ 2 (15) \overline{\mathcal{S}}_{\mathrm{PER}}\left(\omega\right)\frac{1}{L}\sum_{i1}^{L}\hat{S}_{\mathrm{PER}}^i(\omega)\frac{1}{LM}\sum_{i1}^{L}\:\Big|\:\sum_{n0}^{M-1}u_{N}^{i}(n)w(n)\:\mathrm{e}^{-\mathrm{j}\omega n}\:\Big|^{2}\tag{15} SPER(ω)L1i1∑LS^PERi(ω)LM1i1∑L n0∑M−1uNi(n)w(n)e−jωn 2(15)
Welch法结论 Welch方法允许分段数据样本的重叠于是可以得到更多的周期图估计从而进一步减小估计的功率谱密度的方差。而与使用矩形窗实现数据分段的Bartlett方法相比Welch方法可以使用多种窗函数。通过窗函数加权可以减小每个分段起始和结束位置附近的样本在估计中的作用从而减小了相邻样本段之间的相关性。所以Welch 方法可以更好地控制功率谱密度估计的方差特性。
五、MATLAB仿真
设随机过程 u N ( n ) u_N(n) uN(n)由3个实正弦信号加噪声构成为 u N ( n ) ∑ k 1 3 s k ( n ) v ( n ) (16) u_N(n)\sum_{k1}^3s_k(n)v(n)\tag{16} uN(n)k1∑3sk(n)v(n)(16) 其中 s k ( n ) A k cos ( 2 π f k n φ k ) , k 1 , 2 , 3 s_k( n) A_k\cos ( 2\pi f_kn \varphi _k) , k 1, 2, 3 sk(n)Akcos(2πfknφk),k1,2,3,其归一化频率分别是 f 1 0.1 , f 2 0.25 f_10.1,f_20.25 f10.1,f20.25和 f 3 0.27 , φ k f_30.27,\varphi_k f30.27,φk是相互独立并在 [ 0 , 2 π ] \begin{bmatrix}0,2\pi\end{bmatrix} [0,2π]上服从均匀分布的随机相位 v ( n ) v(n) v(n)是均值为0、方差 σ 2 1 \sigma^21 σ21的实高斯白噪声序列。3个信号的信噪比分别为 S N R 1 30 \mathrm{SNR}_130 SNR130dB , S N R 2 30 ,\mathrm{SNR}_230 ,SNR230dB , S N R 3 27 ,\mathrm{SNR}_327 ,SNR327dB。 参数设置如下 通过以上结果可以看出点数对信号分辨率的影响点数加大时可以提高分辨率我这里就不贴图了感兴趣的可以自己下载代码验证。
六、MATLAB详细代码
经典功率谱估计超详细的MATLAB代码BT法和周期图法和bartlett法和welch法
总结
例如以上就是今天要讲的内容本文介绍了BT法、周期图法、Bartlett法、welch法四种经典功率谱估计方法并在MATLAB上面完成了功率谱估计仿真。
