7.png (183.66 KB, 下载次数: 0)
2018-1-25 16:06 上传
图3 实验所用的随机信号
采样点数N分别取128、256、512和1024，周期图法matlab代码如下：
Fs=1000;
f1=50;
f2=125;
f3=135;
N=128;
Nfft=N;
n=0:N-1;
t=n/Fs;
%采用的时间序列
xn=cos(2*pi*f1*t)+1.5*cos(2*pi*f2*t)+cos(2*pi*f3*t)+1.5*randn(size(n));
figure;
plot(n,xn);
grid on;
title('时域信号');
P=10*log10(abs(fft(xn,Nfft).^2)/N);
%Fourier振幅谱平方的平均值，并转化为dB
f=(0:length(P)-1)*Fs/length(P);
%给出频率序列
figure;
plot(f(1:N/2),P(1:N/2));
grid on;
title('功率谱(dB图)');
ylabel('功率谱/dB');
xlabel('频率/Hz');
Pxx=abs(fft(xn,Nfft).^2)/N;
%Fourier振幅谱平方的平均值，并转化为dB
f=(0:length(Pxx)-1)*Fs/length(Pxx);
%给出频率序列
figure;
plot(f(1:N/2),Pxx(1:N/2));
%绘制功率谱曲线
xlabel('频率/Hz');
ylabel('功率谱');
title('周期图 N=128');
grid on;
std(Pxx)^2
N=256;
Nfft=N;
n=0:N-1;
t=n/Fs;
%采用的时间序列
xn=cos(2*pi*f1*t)+1.5*cos(2*pi*f2*t)+cos(2*pi*f3*t)+1.5*randn(size(n));
figure;
plot(n,xn);
grid on;
title('时域信号');
P=10*log10(abs(fft(xn,Nfft).^2)/N);
%Fourier振幅谱平方的平均值，并转化为dB
f=(0:length(P)-1)*Fs/length(P);
%给出频率序列
figure;
plot(f(1:N/2),P(1:N/2));
grid on;
title('功率谱(dB图)');
ylabel('功率谱/dB');
xlabel('频率/Hz');
Pxx=abs(fft(xn,Nfft).^2)/N;
%Fourier振幅谱平方的平均值，并转化为dB
f=(0:length(Pxx)-1)*Fs/length(Pxx);
%给出频率序列
figure;
plot(f(1:N/2),Pxx(1:N/2));
%绘制功率谱曲线
xlabel('频率/Hz');
ylabel('功率谱');
title('周期图 N=256');
grid on;
std(Pxx)^2
N=512;
Nfft=N;
n=0:N-1;
t=n/Fs;
%采用的时间序列
xn=cos(2*pi*f1*t)+1.5*cos(2*pi*f2*t)+cos(2*pi*f3*t)+1.5*randn(size(n));
figure;
plot(n,xn);
grid on;
title('时域信号');
P=10*log10(abs(fft(xn,Nfft).^2)/N);
%Fourier振幅谱平方的平均值，并转化为dB
f=(0:length(P)-1)*Fs/length(P);
%给出频率序列
figure;
plot(f(1:N/2),P(1:N/2));
grid on;
title('功率谱(dB图)');
ylabel('功率谱/dB');
xlabel('频率/Hz');
Pxx=abs(fft(xn,Nfft).^2)/N;
%Fourier振幅谱平方的平均值，并转化为dB
f=(0:length(Pxx)-1)*Fs/length(Pxx);
%给出频率序列
figure;
plot(f(1:N/2),Pxx(1:N/2));
%绘制功率谱曲线
xlabel('频率/Hz');
ylabel('功率谱');
title('周期图 N=512');
grid on;
std(Pxx)^2
N=1024;
Nfft=N;
n=0:N-1;
t=n/Fs;
%采用的时间序列
xn=cos(2*pi*f1*t)+1.5*cos(2*pi*f2*t)+cos(2*pi*f3*t)+1.5*randn(size(n));
figure;
plot(n(1:1000),xn(1:1000));
grid on;
title('时域信号');
P=10*log10(abs(fft(xn,Nfft).^2)/N);
%Fourier振幅谱平方的平均值，并转化为dB
f=(0:length(P)-1)*Fs/length(P);
%给出频率序列
figure;
plot(f(1:N/2),P(1:N/2));
grid on;
title('功率谱(dB图)');
ylabel('功率谱/dB');
xlabel('频率/Hz');
Pxx=abs(fft(xn,Nfft).^2)/N;
%Fourier振幅谱平方的平均值，并转化为dB
f=(0:length(Pxx)-1)*Fs/length(Pxx);
%给出频率序列
figure;
plot(f(1:N/2),Pxx(1:N/2));
%绘制功率谱曲线
xlabel('频率/Hz');
ylabel('功率谱');
title('周期图 N=1024');
grid on;
std(Pxx)^2
以上代码可得到的功率谱分别如图4、图5、图6和图7所示。分辨率能够直观的通过功率谱图形看出，方差的数值由表1给出。
8.png (62.96 KB, 下载次数: 0)
2018-1-25 16:06 上传
图4 N=128时周期图法得到的功率谱
9.png (64.84 KB, 下载次数: 0)
2018-1-25 16:06 上传
图5 N=256时周期图法得到的功率谱
10.png (77.42 KB, 下载次数: 0)
2018-1-25 16:06 上传
图6 N=512时周期图法得到的功率谱
11.png (72.13 KB, 下载次数: 0)
2018-1-25 16:07 上传
图7 N=1024时周期图法得到的功率谱
表1 不同N值得到功率谱的方差值
12.png (19.75 KB, 下载次数: 0)
2018-1-25 16:07 上传 通过上面实验结果的比较，我们很容易发现，周期图法得到的功率谱随着数据点数N的增大，分辨率变大、方差变也大。
二、平均周期图法 周期图法得到的功率谱与我们所期望的"分辨率大、方差小"是矛盾的。为了进一步降低方差，将N个观测样本数据点uN(n)分为L段，每段数据长度为M，分别对每段数据求周期图功率谱估计，然后求平均值，这种方法称平均周期图法。
那么这种方法会如何改善方差呢?下面给出证明：
13.png (10.34 KB, 下载次数: 0)
2018-1-25 16:08 上传 其中：
14.png (60.89 KB, 下载次数: 0)
2018-1-25 16:08 上传 由上我们可以看出，平均周期图法将原来的方差变为原来的1/L，L为分段数。
平均周期图法性能
取采样点数N为1024，分段数分别为8、4、2，平均周期图法的matlab代码如下：
clear;
Fs=1000;
f1=50;
f2=125;
f3=135;
n=0:1/Fs:1;
xn=cos(2*pi*f1*n)+1.5*cos(2*pi*f2*n)+cos(2*pi*f3*n)+1.5*randn(size(n));
N=1024;
Nsec=512;
%数据的长度和FFT所用的数据长度
Pxx1=abs(fft(xn(1:512),Nsec).^2)/Nsec;
%第一段功率谱
Pxx2=abs(fft(xn(513:1000),Nsec).^2)/Nsec;
%第二段功率谱
Pxx=10*log10((Pxx1+Pxx2)/2);
%Fourier振幅谱平方的平均值，并转化为dB
(std((Pxx1+Pxx2)/2))^2
f=(0:length(Pxx)-1)*Fs/length(Pxx);
%给出频率序列
figure;
plot(f(1:Nsec/2),Pxx(1:Nsec/2));
%绘制功率谱曲线
xlabel('频率/Hz');
ylabel('功率谱/dB');
title('N=2*512');
grid on;
N=1024;
Nsec=256;
%数据的长度和FFT所用的数据长度
Pxx1=abs(fft(xn(1:256),Nsec).^2)/Nsec;
%第一段功率谱
Pxx2=abs(fft(xn(257:512),Nsec).^2)/Nsec;
%第二段功率谱
Pxx3=abs(fft(xn(513:768),Nsec).^2)/Nsec;
%第三段功率谱
Pxx4=abs(fft(xn(769:1000),Nsec).^2)/Nsec;
%第四段功率谱
Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4)/4);
%Fourier振幅谱平方的平均值，并转化为dB
std((Pxx1+Pxx2+Pxx3+Pxx4)/4)^2
f=(0:length(Pxx)-1)*Fs/length(Pxx);
%给出频率序列
figure;
plot(f(1:Nsec/2),Pxx(1:Nsec/2));
%绘制功率谱曲线
xlabel('频率/Hz');
ylabel('功率谱/dB');
title('N=4*256');
grid on;
N=1024;
Nsec=128;
%数据的长度和FFT所用的数据长度
Pxx1=abs(fft(xn(1:128),Nsec).^2)/Nsec;
%第一段功率谱
Pxx2=abs(fft(xn(129:256),Nsec).^2)/Nsec;
%第二段功率谱
Pxx3=abs(fft(xn(257:384),Nsec).^2)/Nsec;
%第三段功率谱
Pxx4=abs(fft(xn(385:512),Nsec).^2)/Nsec;
%第四段功率谱
Pxx5=abs(fft(xn(513:640),Nsec).^2)/Nsec;
%第五段功率谱
Pxx6=abs(fft(xn(641:768),Nsec).^2)/Nsec;
%第六段功率谱
Pxx7=abs(fft(xn(769:896),Nsec).^2)/Nsec;
%第七段功率谱
Pxx8=abs(fft(xn(897:1000),Nsec).^2)/Nsec;
%第八段功率谱
Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7+Pxx8)/8);
%Fourier振幅谱平方的平均值，并转化为dB
std((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7+Pxx8)/8)^2
f=(0:length(Pxx)-1)*Fs/length(Pxx);
%给出频率序列
figure;
plot(f(1:Nsec/2),Pxx(1:Nsec/2));
%绘制功率谱曲线
xlabel('频率/Hz');
ylabel('功率谱/dB');
title('N=8*128');
grid on;
以上代码可得到的功率谱分别如图8、图9、图10所示。分辨率能够直观的通过功率谱图形看出，方差的数值由表2给出。
15.png (55.28 KB, 下载次数: 0)
2018-1-25 16:09 上传
图8 L=8时平均周期图法得到的功率谱
16.png (62.03 KB, 下载次数: 0)
2018-1-25 16:09 上传
图9 L=4时平均周期图法得到的功率谱
17.png (73.05 KB, 下载次数: 0)
2018-1-25 16:09 上传
图10 L=2时平均周期图法得到的功率谱
表2 不同L值得到功率谱的方差值
18.png (16.58 KB, 下载次数: 0)
2018-1-25 16:09 上传 L=1时，平均周期图法退化为周期图法。通过上面实验结果的比较，我们很容易发现，平均周期图法得到的功率谱随着分段数L变大，方差变小，但分辨率变小。
当观测样本序列数据个数N固定时，要降低方差需要增加分段数L。当N不大时分段长度M取值较小，则功率谱分辨率降低到较低的水平。若分段数L固定时，增加分辨率需要分段长度M，则需要采集到更长的检测数据序列。实际中恰恰是检测样本序列长度不足。
三、修正的平均周期图法 由于实际检测中样本序列长度是有限的。对现有数据长度N，如果能获得更多的段数分割，将会得到更小的方差。允许数据段间有重叠部分，来得到更多的段数。对段间重叠长度的选取，最简单是取为段长度M的一半。由平均周期图法可知更多的段数可以进一步降低方差。
数据截断的过程中相当于数据加矩形窗，矩形窗幅度较大的旁瓣会造成"频谱泄漏"。我们分段时采取的窗函数更为多样(三角窗,海明窗等), 以减小截断数据(加矩形窗)窗函数带来的影响[2]。
修正平均周期图法的matlab代码如下：
火狐截图_2018-01-25T08-14-50.942Z.png (132.47 KB, 下载次数: 0)
2018-1-25 16:15 上传
分别采用矩形窗、Blackman窗和Hamming窗，上述代码可得到的功率谱如图11所示。
19.png (95.65 KB, 下载次数: 1)
2018-1-25 16:10 上传
图11 不同窗函数的修正平均周期图法得到的功率谱
由上可以发现，矩形窗的分辨率最高，但是方差也最大，这是由于矩形窗频谱主瓣最窄，分辨率因此最高，旁瓣也高，导致频谱泄漏最严重，方差最大。
综而言之，周期图法获得的功率谱随着样本点数越多，分辨率越大、方差越大;平均周期图法以牺牲分辨率来进一步改善方差;修正的平均周期图法允许段的重叠来进一步增大分段数、或者分段数相同，每段样本点数变多。无论是哪种方法都没有彻底结局方差与分辨率之间的矛盾。
四、相关功率谱估计法-BT法 如前所述，要提高功率谱估计的分辨率，必须增加数据序列的长度N，但是较长的数据序列，由噪声引起的随机性得到更为充分的体现-较大的方差。事实上，当N无穷大时，方差为一非零常数，即周期图法无法实现功率谱的一致估计。下面我们来介绍相关功率谱估计法(BT法)，它是一致估计。
维纳辛钦定理指出，随机信号的相关函数与它的功率谱是一对傅里叶变换对。BT法就是基于这个原理。先由观测数据估计出自相关函数，然后求自相关函数的傅立叶变换，以此变换作为对功率谱的估计，也称为间接法。BT法要求信号长度N以外的信号为零，这是BT法的局限性。BT法可以表述为：
自相关函数：
20.png (11.48 KB, 下载次数: 0)
2018-1-25 16:11 上传 功率谱：
21.png (1.47 KB, 下载次数: 0)
2018-1-25 16:11 上传 取采样点数N分别为128、256、512和1024的BT法代码如下：
QQ截图20180125161349.png (110.19 KB, 下载次数: 1)
2018-1-25 16:13 上传
以上代码可得到的功率谱如图12、图13、图14和图15所示。
22.png (53.27 KB, 下载次数: 0)
2018-1-25 16:11 上传
图2-10 N=128时，BT法得到的功率谱
23.png (53.7 KB, 下载次数: 0)
2018-1-25 16:11 上传
图2-11 N=256时，BT法得到的功率谱
24.png (70.33 KB, 下载次数: 0)
2018-1-25 16:11 上传
图2-12 N=512时，BT法得到的功率谱
25.png (65.07 KB, 下载次数: 0)
2018-1-25 16:11 上传
图2-13 N=1024时，BT法得到的功率谱
由上面实验可以发现，M随着N的增大而增大时，分辨率提高，方差变大。BT法仍然没有解决分辨率与方差之间的矛盾，但是BT法得到的功率谱当N为无穷大时，方差会趋向于零，即为一致估计[2]。
周期图法与BT法的关系 相关函数可以写成卷积形式：
26.png (7.97 KB, 下载次数: 0)
2018-1-25 16:12 上传 设序列uN(n)的傅立叶变换为UN(ω)，则当M=N-1时，功率谱的估计可表示为：
27.png (6.04 KB, 下载次数: 0)
2018-1-25 16:12 上传 上式可以看出周期图法可以看作BT法在取M=N-1时的特例。
结 论 本文通过Matlab仿真，以一个具体的随机信号为例，简单介绍了周期图法、平均周期图法、修正的平均周期图法以及BT法的基本原理，并对这些方法的性能进行分析。可以看出，无论是周期图法及其改进算法还是BT法都没有从根本上解决分辨率与方差的矛盾。
经典功率谱估计是利用傅里叶变换估计功率谱，而我们之前分析随机信号不满足傅里叶变换的条件，所以经典功率谱估计方法不得不从无限长数据点截取有限长数据点，加入限制条件(周期图法实际上假定N点外数据周期重复、BT法假定N点外数据为零)来"强制"作傅里叶变换，这也是造成它局限性的原因。
参考资料
[1] 朱哲,钟宏伟. 非平稳随机信号分析处理方法研究[J] 安徽电子信息技术学院学报 2008.6:28-28
[2] 皇甫堪. 现代数字信号处理[M]. 电子工业出版社
本文根据博客园lulujianjie的博文整理
原文：http://www.cnblogs.com/jacklu/p/5140913.html
扫一扫
827
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
