经典功率谱估计——BT法探究

正弦信号的自相关。

广义上的自相关，或者说最一般的自相关，在matlab 中可用卷积conv 和相关函数xcorr 来实现。

序列x(n)和y(n)的互相关公式为：
∑∞
-∞
=-=
n
xy m n y n x m )()()(φ
同理，自相关公式为：
∑∞
-∞
=-=
n
xx m n x n x m )()()(φ
有以下程序验证： close all clear,clc
Fs=1000; %采样率 f = 5; % f = .5;
% f1 = 2; %信号基频 % f2 = 2.5;
N = 1; % 周期数 % t = N/f;
t=N; % 信号时长 s
% t = N/f1; % N 个周期的时间,针对小频率信号，在此时间内，大频率信号周期更多 n=0:1/Fs:t -1/Fs; % 采样时间点，刚好采N 个周期
len = length(n); % 信号点数，也是FFT 变换点数，即采集多少点就做多少点的FFT
y = sin(2*pi*f*n); % 采集到的离散信号
% y = sin(2*pi*f1*n) + sin(2*pi*f2*n); % 采集到的离散信号 plot(y)
y1 = fliplr(y);
y_corr = conv(y,y1);
% [y_corr1,lags] = xcorr(y,100,'unbiased'); y_corr1 = xcorr(y,y);
figure;plot(y_corr)
% hold on
figure;plot(y_corr1)
相关的图像是一模一样的！
但是，在BT 法求功率谱中，信号的自相关却不是这样求的。

由于信号的功率谱与自相关函数互为傅里叶变换关系，因此，信号的功率谱估计可以先通过对自相关函数进行估计，再对其进行傅里叶变换即可，这种方法称为自相关函数法，由Blackman 和Tukey 于1958年提出，故也称为BT 法，又称为间接法。

设N 个样本序列{xn}的值为x(0),x(1),…,x(N -1),现需要用此N 个数据来估计自相关函数)(m xx φ.由于xn 只能观察到0<<n<<N -1的N 个值，而n<0和n>N -1时的xn 值是未知的，一般只能假定为0.根据自相关函数的定义得到：
∑-+=
10
)()(1
)(N xx m n
x n x N
m φ
由于x(n)只有N 个观测值，因此对于每个固定的延迟m ，可以利用的数据只有（N -|m|）个，且在[0，N -1]范围内，所以实际计算)(m xx φ为：
∑--+=
1
||0
)()(1
)(m N xx m n
x n x N
m φ
考虑乘积项的长度，自相关序列的估计为：
1||其中,)()(||1
)(1
||0
-≤+-=∑--N m m n x n x m N m m N xx φ
式中，m 取绝对值是因为)(m xx φ=)(m xx -φ，m 为负值时上式仍适用。

规定的求和上下限的原则是保证充分利用全部数据。

在matlab 中，上述这种相关方式成为无偏相关。

只需要在调用相关函数xcorr 时选中一个参数’unbiased’即可。

有以下程序： clf;f=10;
N=1000;Fs=500;
n=[0:N-1];t=n/Fs;
Lag=100;
randn(‘state’,0);
x=sin(2*pi*f*t)+0.6*randn(1,length(t));
[c,lags]=xcorr(x,Lag,’unbiased’);
subplot(311),plot(t,x);
xlabel(‘t/s’);ylabel(‘x(t)’);grid;
legend(‘含噪声的信号x(t)’);
subplot(312),plot(lags/Fs,c);
xlabel(‘t/s’);ylabel(‘Rxx(t)’);
legend(‘信号的自相关Rxx’);
grid;
Pxx=fft(c,length(lags));
fp=(0:lenth(Pxx)-1)’*Fs/length(Pxx);
Pxmag=abs(Pxx);
subplot(313)
plot(fp(1:length(Pxx)/2),Pxmag(1:length(Pxx)/2)); xlabel(‘f/Hz’);ylabel(‘|Pxx|’);grid;
legend(‘信号功率谱Pxx’)；
结果：
功率谱放大：
可见主频是在10Hz 。

上述是指定了自相关时信号点数的延迟区间为[-100,100],则自相关点数为200-1=199点。

也可以不设置，直接在信号全长度的区间依次延迟运算时平移的信号，即： [c,lags]=xcorr(x,’unbiased’); 得到结果：
00.20.40.60.8
1 1.
2 1.4 1.6 1.82-505
t/s
x (t )
-0.2
-0.15-0.1-0.05
00.050.10.150.2-101
t/s
R x x (t
)
50100
150200250
050100
f/Hz
|P x
x
|
功率谱放大：
可以看出，功率谱上主频依然是10.
有个问题：这两种计算得到的功率谱峰值大小不等，这该怎么解决呢？
00.20.40.60.8
1 1.
2 1.4 1.6 1.82-505
t/s
x (t )
-2
-1.5-1-0.5
00.51 1.52-101
t/s
R x x (t
)
50100
150200250
0500
f/Hz
|P x x
|。