用matlab实现功率谱仿真分析
- 格式：docx
- 大小：590.58 KB
- 文档页数：6
1
功率谱估计性能分析及其MATLAB实现
一、 经典功率谱估计分类简介
1． 间接法
根据维纳-辛钦定理，1958年Blackman和Turkey给出了这一方法的具体实现，即先由N个观察值xN(n)，估计出自相关函数rx(m)，求自相关函数傅里叶变换，以此变换结果作为对功率谱Px(w)的估计。
2． 直接法
直接法功率谱估计是间接法功率谱估计的一个特例，又称为周期图法，它是把随机信号的N个观察值xN(n)直接进行傅里叶变换，得到XN(ejw)，然后取其幅值的平方，再除以N，作为对功率谱Px(w)的估计。
3． 改进的周期图法
将N点的观察值分成L个数据段，每段的数据为M，然后计算L个数据段的周期图的平均P̅per(w)，作为功率谱的估计，以此来改善用N点观察数据直接计算的周期图P̂per(w)的方差特性。根据分段方法的不同，又可以分为Welch法和Bartlett法。
Welch法
所分的数据段可以互相重叠，选用的数据窗可以是任意窗。
Bartlett法
所分的数据段互不重叠，选用的数据窗是矩形窗。
二、 经典功率谱估计的性能比较
1． 仿真结果
为了比较经典功率谱估计的性能，本文采用的信号是高斯白噪声加两个正弦信号，采样率Fs=1000Hz，两个正弦信号的频率分别为f1=200Hz，f2=210Hz。所用数据长度N=400.
仿真结果如下：
2
(a)
(b)
(c)
(d)
(e)
(f)
Figure1 经典功率谱估计的仿真结果
Figure1(a)示出了待估计信号的时域波形；
Figure2(b)示出了用该数据段直接求出的周期图，所用的数据窗为矩形窗；
Figure2(c)是用BT法（间接法）求出的功率谱曲线，对自相关函数用的平滑窗为矩形窗，长度M=128，数据没有加窗；
Figure2(d)是用BT法（间接法）求出的功率谱曲线，对自相关函数用的平滑窗为Hamming窗，长度M=64，数据没有加窗；
Figure2(e)是用Welch平均法求出的功率谱曲线，每段数据的长度为64点，重叠32点，使用的Hamming窗；
Figure2(f)是用Welch平均法求出的功率谱曲线，每段数据的长度为100点，重叠48点，使用的Hamming窗； 3
2． 性能比较
1) 直接法得到的功率谱分辨率最高，但是方差性能最差，功率谱起伏剧烈，容易出现虚假谱峰；
2) 间接法由于使用了平滑窗对直接法估计的功率谱进行了平滑，因此方差性能比直接法好，功率谱比直接法估计的要平滑，但其分辨率比直接法低。
3) Welch平均周期图法是三种经典功率谱估计方法中方差性能最好的，估计的功率谱也最为平滑，但这是以分辨率的下降及偏差的增大为代价的。
3． 关于经典功率谱估计的总结
1) 功率谱估计，不论是直接法还是间接法都可以用FFT快速计算，且物理概念明确，因而仍是目前较常用的谱估计方法。
2) 谱的分辨率较低，它正比于2π/N，N是所使用的数据长度。
3) 方差性能不好，不是真实功率谱的一致估计，且N增大时，功率谱起伏加剧。
4) 周期图的平滑和平均是和窗函数的使用紧密关联的，平滑和平均主要是用来改善周期图的方差性能，但往往又减小了分辨率和增加了偏差，没有一个窗函数能使估计的功率谱在方差、偏差和分辨率各个方面都得到改善，因此使用窗函数只是改进估计质量的一个技巧问题，并不能从根本上解决问题。
三、 AR模型功率谱估计
1． AR模型功率谱估计简介
AR模型功率谱估计是现代谱估计中最常用的一种方法，这是因为AR模型参数的精确估计可以用解一组线性方程（Yule-Walker方程）的方法求得。其核心思想是：将信号看成是一个p阶AR过程，通过建立Yule-Walker方程求解AR模型的参数，从而得到功率谱的估计。
由于已知的仅仅是长度有限的观测数据，因此AR模型参数的求得，通常是首先通过某种算法求得自相关函数的估计值，进而求得AR模型参数的估计值。常用的几种AR模型参数提取方法有：
1) 自相关法
假定观测数据区间之外的数据为0，在均方误差意义下使得数据的前向预测误差最小。
由此估计的自相关矩阵式正定的，且具有Toeplitz性，可以用Levison-Durbin算法求解。
2) 协方差法
不作观测数据区间之外的数据为0的假设，在均方误差意义下使得数据的前向预测误差最小。由此估计的自相关矩阵式半正定的，且不具有Toeplitz性，得到的AR模型可能不稳定。
3) 修正的协方差法
不作观测数据区间之外的数据为0的假设，在均方误差意义下使得数据的前向预测误差与后向预测误差之和最小。由此估计的自相关矩阵式半正定的，且不具有Toeplitz性，得到的AR模型可能不稳定。但得到的一阶AR模型是稳定的。
4
4) Burg法
在约束AR模型的参数满足Levison-Durbin递归条件的前提下，在均方误差意义下使得数据的前向预测误差与后向预测误差之和最小。得到的AR模型是稳定的，但有时可能出现谱线分裂现象。
仍然用前面的仿真信号，取AR模型的阶数p=48，用上述各种AR模型参数提取方法估计的功率谱如figure2所示。
Figure2 AR模型功率谱估计的仿真结果
2． AR模型功率谱估计与经典功率谱估计性能比较
分别采用直接法和AR模型功率谱估计中的自相关法得到的上述信号的功率谱估计如figure3所示：
Figure3 AR模型功率谱估计与经典功率谱估计比较 5
小结：
1) 由于AR模型是一个有理分式，因而估计出的谱要比经典法的谱平滑。
2) AR谱估计的一些方法隐含着数据和自相关函数的外推，可获得高的分辨率。
3． 关于AR模型功率谱估计总结
Figure4给出了阶数对AR模型功率谱估计的影响，采用的是AR模型功率谱估计中的Burg法。
Figure4 不同阶数的AR模型功率谱估计
小结：
阶数越高，得到的谱的分辨率也越高，但方差也越大，将会产生更多的虚假谱峰；
四、 附录
本文所用来仿真的MATLAB程序如下：
% Research On Classic PSD Estimate Methods
% Author: Chen Feiqiang
% Date: 2010-12-14
%% Generate Signal
clear all; close all; clc;
randn('state',0);
Fs = 1000; % sample frequency
t = 0:1/Fs:.3;
sigma = 1; % noise variance
x = cos(2*pi*t*200) + cos(2*pi*t*210) + sigma*randn(size(t)); % generate signal:
% cosine + noise
figure;
plot(t,x), xlabel('t'), ylabel('x'); 6
title('Signal In Time Domain');grid on;
%% Estimate PSD using periodogram method
window = rectwin(length(x)); % window function used
xn = x'.*window;
index = nextpow2(length(x));
N = pow2(index);
Xw = fft(xn, N); % do N-points FFT
Pw = Xw.*conj(Xw)/N; % Calculate PSD
k = 0:floor(N/2 - 1);
figure;
plot(k*Fs/N, 10*log10(Pw(k+1)/max(Pw)));
title('Periodogram PSD Estimate');
xlabel('Frequency(Hz)');ylabel('Normalized PSD(dB)'),grid on;
hold on
%% Estimate PSD using BT method
window_a = rectwin(length(x)); % window function for data x(n)
xn = x'.*window_a;
Rx = xcorr(xn); % auto-correlation function estimate
N = length(Rx);
M = floor(N/4); % the length of smooth window
% M = 100;
window_v = rectwin(M); % smooth window for BT method
RxWin = Rx(1:M).*window_v; % smooth window multiply auto-correlation
function
Pw = abs(fft(RxWin, N));
k = 0:floor(N/2 - 1);
figure;
plot(k*Fs/N, 10*log10(Pw(k+1)/max(Pw)),'r');
title('BT Method PSD Estimate');
xlabel('Frequency(Hz)');ylabel('Normalized PSD(dB)'),grid on;
%% Estimate PSD using Welch method
window = 32; % the length of each segment
noverlap = 8; % overlap number for each segment
nfft = pow2(nextpow2(length(x))); % nfft-points FFT for each segment
[Pxx,f] = pwelch(x,window,noverlap,nfft,Fs); % call pwelch function
figure;
plot(f,10*log10(Pxx/max(Pxx)),'g');
title('Pwelch Method PSD Estimate');
xlabel('Frequency(Hz)');ylabel('Normalized PSD(dB)'),grid on;
