经典功率谱估计及其实现

又到周五了，仿真实现了一半，回头来把这篇文章写了吧，两周前我决定写这篇文章时，对功率谱理解是一知半解的，现在不断地仿真、看论文，理解的比以前深了一点吧，一切都会好起来的~
参考书籍：
《现代信号处理》安颖、崔东艳著
《现代信号处理教程》胡广书著
《数字信号处理原理及其Matlab实现》从玉良编著

一、信号处理引言

作为信号处理方向的学生，经历过本科生和研究生的教育，回头来看信号处理，其实感觉脉络还是挺清晰的，只是学习的时候“只缘身在此山中”，一直是雾里看花，现在参考书籍和自己的理解，首先来明确这些知识的结构脉络。
对于信号的处理，总的来说有两类方法：一类是时域处理，比如使用维纳-卡尔曼滤波、自适应滤波这些；另一类是频域处理，此时如果针对确定性信号，当这些信号满足绝对可积或者可和时，直接转换到频域处理。如果针对随机信号（比较大自然中这些信号居多），采用随机信号的处理方法（比如协方差、自相关矩阵等）对信号进行功率谱估计。
功率谱的估计主要包括经典谱估计和现代谱估计两种方法。经典谱估计由于分析的输入信号为有限长，不能解决时间分辨率和频率分辨率的矛盾（测不准原理）、方差性能差，于是引入了现代谱估计：参数模型法具体表现为AR、MA、ARMA等方法估计出随机信号的功率谱。非参数模型法主要包括Pisarenko谱波分解等方法。
当然，还有一些其他的知识，比如滤波器、S变换等其他知识点，不过一切知识都是作为工具为我们服务的，我们只需要有重点的学习即可。

二、经典功率谱估计方法

经典功率谱估计采用的是传统傅里叶变换分析方法（又称线性谱估计），主要分为自相关法（间接法）和周期图法（直接法）两种。自相关法在1985年提出，先估计自相关函数，再计算功率谱。周期图法直接对观测数据进行快速傅里叶变换，得到功率谱。优点是可以直接FFT快速计算，所以应用比较广泛。
经典谱估计优点是计算效率高，缺点是频率分辨率低，常用于频率分辨率要求不高的场合。
放图：

自相关法

自相关法是根据维纳-辛钦定理，通过相关函数计算功率谱的。

周期图法

Schuster于1899年提出，把N点观测数据看做能量有限的信号，直接对其进行傅里叶变换，然后取其模值的平法，并除以N，得到观测数据真实的功率谱的估计。
质量分析：
（1）周期图法是功率谱的有偏估计。产生偏移的原因一方面是局部平均中主瓣的模糊租用，模糊的结果使谱估计的分辨率下降；另一方面是由于旁瓣的泄露。
（2）频率分辨率低。原因是傅里叶变换域是无限长的，而观测数据是有限长的，相当于将信号在时域添加了矩形窗，在频域真正功率谱卷积一个抽样函数sinc。
总结：
周期图法很简单，直接使用谱估计性能很差，因此在使用时一般使用改进的周期图法。

四种常见的改进周期图法

周期图法相当于一个矩形窗对时域信号的加窗，然后分辨率大约于数据长度N成反比。这里提出四种常见的改进周期图法。

平均周期图法

对一个随机变量观测时，得到L组独立数据，每组数据长为M，对每一组求PSD，之后将L个均值加起来求平均。这样得到的均值，其方差是原来的1/L。
质量分析：
平均周期图法仍然是有偏估计，偏移和每组数据长度M有关。由于每段FFT的长度变为M，分辨率更低（Fs/M），因此，平均周期图法以分辨率的降低换区了估计方差的减少。

窗函数法

窗函数法是将长度为N的观测数据乘以同一长度的数据窗w，数据加窗后，谱估计值的数学期望等于真实谱值与窗谱函数的平方卷积，因而是有偏估计。
质量分析：
选择低旁瓣的数据窗可使得杂散响应减少，但旁瓣的降低必然使得主瓣加宽，从而降低分辨率。
数据加窗后，谱估计值的方差>=谱估计值的平方，且不随数据长度增大而减少。
数据加窗可以降低谱估计值的旁瓣，但是降低了谱估计的分辨率。

修正的周期图平均法

又叫Bartlett法，首先把长度为N的数据分成L段，每段数据长为M，则N=LM。然后把窗函数w加到每段数据上，求出每段的周期图，之后对每段周期图进行评价。
质量分析：
Bartlett法在计算周期图前，先对各数据段加窗，是平均周期图法的估计方差减少，但是分辨率降低。

加权交叠平均法

又称Welch法。对Bartlett法的改进。首先，分段时相邻两段可以重叠，其次，窗函数使用汉宁窗或汉明窗，通过改进，达到进一步减小方差的目的。

总结：

经典功率谱虽然实现简单，计算速度快，但是具有方差性能较差、分辨率较低、旁瓣泄露等缺点。

三、仿真比较

（1）利用加矩形窗周期图法(自带API)实现功率谱估计。
代码：

%利用加矩形窗周期图法(自带API)实现功率谱估计
Fs = 500;                       %采样频率
NFFT = 1024;
n = 0:1/Fs:1;
vx = randn(1,length(n));        %产生含有噪声的序列
x = 4 * sin(2*pi*100*n) - 2*sin(2*pi*10*n) + vx;
window = boxcar(length(n));     %矩形窗
[Pxx,f]=periodogram(x,window,NFFT,Fs);  %将结果转换为dB
figure();
plot(f,10*log10(Pxx));title('加矩形窗周期图法实现功率谱估计');

效果：

（2）利用加三角形窗周期图法(自带API)实现功率谱估计

%利用加三角形窗周期图法(自带API)实现功率谱估计
Fs = 500;                       %采样频率
NFFT = 1024;
n = 0:1/Fs:1;
vx = randn(1,length(n));        %产生含有噪声的序列
x = 4 * sin(2*pi*100*n) - 2*sin(2*pi*10*n) + vx;
window = bartlett(length(n));     %三角形窗
[Pxx,f]=periodogram(x,window,NFFT,Fs);  %将结果转换为dB
figure();
plot(f,10*log10(Pxx));title('加三角形窗周期图法实现功率谱估计');

效果：

（3）自相关函数法实现功率谱估计

% 自相关函数法实现功率谱估计
Fs = 500;                       %采样频率
NFFT = 1024;
n = 0:1/Fs:1;
vx = randn(1,length(n));        %产生含有噪声的序列
x = 4 * sin(2*pi*100*n) - 2*sin(2*pi*10*n) + vx;
Cx = xcorr(x,'unbiased');
Cxk = fft(Cx,NFFT);
Pxx = abs(Cxk);
t = 0:round(NFFT/2 -1);
k = t*Fs/NFFT;
P = 10*log(Pxx(t+1));
figure();
plot(k,P);title('自相关函数法实现功率谱估计');

效果：

（4）利用FFT算法实现功率谱估计

%利用FFT算法实现功率谱估计
Fs = 500;
NFFT = 1024;
n = 0:1/Fs:1;
vx = randn(1,length(n));
x = 4 * sin(2*pi*100*n) - 2*sin(2*pi*10*n) + vx;
X = fft(x,NFFT);
Pxx = abs(X).^2/length(n);
t = 0:round(NFFT/2 -1);
k = t*Fs/NFFT;
P = 10*log(Pxx(t+1));
figure();
plot(k,P);title('利用FFT算法实现功率谱估计');

效果：

（5）利用FFT算法实现平均周期图法

% 平均周期图法(数据不重叠分段)
Fs = 1000;
NFFT = 1024;
Nsec = 256;
n = 0:NFFT-1;
t = n/Fs;
vx = randn(1,NFFT);
x = 4 * sin(2*pi*50*t) - 2*sin(2*pi*120*t) + vx;Pxx1 = abs(fft(x(1:256),Nsec).^2)/Nsec;
Pxx2 = abs(fft(x(256:512),Nsec).^2)/Nsec;
Pxx3 = abs(fft(x(512:768),Nsec).^2)/Nsec;
Pxx4 = abs(fft(x(768:1024),Nsec).^2)/Nsec;
Pxx = 10*log((Pxx1+Pxx2+Pxx3+Pxx4)/4);
f = (0:length(Pxx)-1)*Fs/length(Pxx);
subplot(2,1,1);
plot(f(1:Nsec/2),Pxx(1:Nsec/2));
xlabel('Frequency/Hz');ylabel('Powerspectrum/dB');
title('averaged periodogram(non-overlapping)N=4*256');
grid on%平均周期图法（数据重叠分段)
Pxx1 = abs(fft(x(1:256),Nsec).^2)/Nsec;
Pxx2 = abs(fft(x(129:284),Nsec).^2)/Nsec;
Pxx3 = abs(fft(x(257:512),Nsec).^2)/Nsec;
Pxx4 = abs(fft(x(385:640),Nsec).^2)/Nsec;
Pxx5= abs(fft(x(513:768),Nsec).^2)/Nsec;
Pxx6 = abs(fft(x(641:896),Nsec).^2)/Nsec;
Pxx7 = abs(fft(x(769:1024),Nsec).^2)/Nsec;
Pxx = 10*log((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7)/7); %求平均并转换为dB
f = (0:length(Pxx)-1)*Fs/length(Pxx);
subplot(2,1,2);
plot(f(1:Nsec/2),Pxx(1:Nsec/2));
xlabel('Frequency/Hz');ylabel('Powerspectrum/dB');
title('averaged periodogram(overlapping)N=1024');
grid on

效果：

（5）利用pwelch实现Welch法实现平均周期图法

clear;clc;close all;
% 自己实现的平均周期图平均法实现功率谱估计
Fs = 1000;
nfft = 256;
n = 0:1/Fs:1;
vx = randn(1,length(n));
x = 4 * sin(2*pi*100*n) - 2*sin(2*pi*10*n) + vx;
w = hanning(nfft)';
Pxx = ( abs(fft(w.*x(nfft*0/2+1:nfft*2/2))).^2+...abs(fft(w.*x(nfft*1/2+1:nfft*3/2))).^2+...abs(fft(w.*x(nfft*2/2+1:nfft*4/2))).^2+...abs(fft(w.*x(nfft*3/2+1:nfft*5/2))).^2+...abs(fft(w.*x(nfft*4/2+1:nfft*6/2))).^2+...abs(fft(w.*x(nfft*5/2+1:nfft*7/2))).^2 )/(norm(w)^2*6);
f = (0:(nfft-1))/nfft*Fs;
PXX = 10*log10(Pxx);
figure();
plot(f,PXX);
xlabel('频率/Hz');
ylabel('功率谱/dB');
title('自己实现的平均周期图平均法实现功率谱估计');
grid;%利用pwelch实现Welch法平均周期图平均法实现功率谱估计
Fs = 500;
NFFT = 1024;
n = 0:1/Fs:1;
vx = randn(1,length(n));
x = 4 * sin(2*pi*100*n) - 2*sin(2*pi*10*n) + vx;
window1 = boxcar(100);
window2 = hamming(100);
window3 = blackman(100);
noverlap = 20;%数据重叠
[Pxx1,f1]=pwelch(x,window1,noverlap,NFFT,Fs);
[Pxx2,f2]=pwelch(x,window2,noverlap,NFFT,Fs);
[Pxx3,f3]=pwelch(x,window3,noverlap,NFFT,Fs);
PXX1= 10*log10(Pxx1);
PXX2= 10*log10(Pxx2);
PXX3= 10*log10(Pxx3);
figure()
title('----')
subplot(3,1,1);
plot(f1,PXX1);
title('矩形窗');
subplot(3,1,2);
plot(f2,PXX2);
title('汉明窗');
subplot(3,1,3);
plot(f3,PXX3);
title('布莱克曼窗');