因此，在FFT问世之前，相关法是最常用的谱估计方法。
当FFT问世后，情况有所变化。
因为截断后的）（nxN可视作能量信号，由相关卷积定理可得
这就将相关化为线性卷积，而线性卷积又可以用快速卷积来实现。
我们可对上式两边取（2N-1）点DFT，则有
于是将时域卷积变为频域乘积，用快速相关求自相关函数估值的完整方案如下：
1.对N长）（nxN的充（N-1）个零，成为（2N-1）长的。
2.求（2N-1）点的FFT，得
3、求
4、求（2N-1）点的IFFT
（二）运算简要框图
输出
矩形窗截断
相关法谱估计运算简要框图
图中快速相关的输出时从-（N-1）到（N-1）的2N-1点，加窗后截取的是-（M-1）到（M-1）的
，最后做（2M-1）点FFT，即可得到结果。
（三）程序示例
程序的主要思路就是按照运算框图一步一步进行计算，下面附程序并进行简要解释：
N=512,n=0:
N-1;%N是FFT的变换区间
xn=2*cos（2*pi*0.2*n）+cos（2*pi*0.4*n）+2*randn（size（n））;
%产生加有均值为0，方差为1的AWGN的信号
Xk=fft（xn,1024）;%进行2N-1点FFT,系统会自动补0
Sk=abs（Xk）.*（abs（Xk））./N;%取频谱幅度的平方，并除以N，以此作为对xn真实功率谱的估计
Rn=ifft（Sk）;
Sk1=fft（Rn,512）;
figure
（2）
subplot（2,1,1）;
plot（n/N,Sk1）;
ylabel（'Sk'）
title（'（b）相关法估计功率谱密度'）
Sk2=10*log（Sk1）;%对估计出的Sk取对数，使画出的图更加突出特点
subplot（2,1,2）;
plot（n/N,Sk2）;
ylabel（'10log（PSD）'）
（四）结果分析
下面是程序运行后的结果
从上图中我们可以较为明显的看到信号中有两个频率分量，一个在0.2处，一个在0.4处，与产生的信号相一致。
但是我们不难看出，估计出的功率谱谱线非常不平坦，有很多起伏。
三、周期图法谱估计
（一）算法原理简介
周期图法又称直接法。
它是从随机信号x（n）中截取N长的一段，把它视为能量有限x（n）真实功率谱
的估计的抽样。
其具体步骤如下：
第一步：
由获得的N点数据构成有限长序列直接求傅里叶变换，得频谱。
第二步：
取频谱幅度的平方，并除以N，以此作为对x（n）真实功率谱的估计。
事实上，周期图法谱估计与自相关法谱估计的差异只是估计自相关函数的方法不同。
（二）运算简要框图
矩形窗（长度N）截断
X（n）
图中用FFT来代替傅里叶变换
（三）程序示例
N=512,n=0:
N-1;
xn=2*cos（2*pi*0.2*n）+cos（2*pi*0.4*n）+2*randn（size（n））;%产生加有均值为0，方差为1的AWGN的信号
Xk1=fft（xn,512）;%进行N点FFT
Sk3=abs（Xk1）.*（abs（Xk1））./N;%取频谱幅度的平方，并除以N，以此作为对xn真实功率谱的估计
Sk4=10*log（Sk3）;%对估计出的Sk取对数，使画出的图更加突出特点
figure（3）
subplot（2,1,1）;
plot（n/N,Sk3）;
title（'（c）周期图法估计功率谱密度'）
ylabel（'Sk'）
subplot（2,1,2）;
plot（n/N,Sk4）;
ylabel（'10log（PSD）'）
（四）结果分析
下面是程序运行后的结果
从上图中我们同样可以看到信号中有两个频率分量，一个在0.2处，一个在0.4处，与产生的信号相一致。
但是我们不难看出，估计出的功率谱谱线与相关法功率谱估计一样非常不平坦，有很多起伏。
四、Bartlett法功率谱估计
（一）算法原理简介
当我们用相关法或者周期图法对信号的功率谱进行估计时，都不是对
的一致估计，主要原因是方差大。
于是就产生了周期图法的改进。
改进的主要途径是平滑和平均。
平滑是用一个适当的窗函数与计算的功率谱进行卷积，是谱线平滑。
这种方法的出的谱估计是无偏的，方差也小，但分辨率下降。
平均就是将截取的数据段再分成L个小段，分别计算功率谱后取功率谱的平均。
因为L个平均的方差比随机变量的单独方差小L倍，所以当L趋于无穷时，L个平均的方差趋于零，可以达到一致估计的目的。
（二）运算简要框图
矩形窗截断
X（n）分成L小段输出
（三）程序示例
%L=2时bartlett法
N=256,n=0:
255;
x1n=2*cos（2*pi*0.2*n）+cos（2*pi*0.4*n）+2*randn（size（n））;%产生加有均值为0，方差为1的AWGN的信号的前半段
Xk1=fft（x1n,N）;%进行N点FFT
Sk5=abs（Xk1）.^2./N;%取频谱幅度的平方，并除以N，以此作为对xn真实功率谱的估计
n=256:
511;
x2n=2*cos（2*pi*0.2*n）+cos（2*pi*0.4*n）+2*randn（size（n））;%产生加有均值为0，方差为1的AWGN的信号后半段
Xk2=fft（x2n,N）;%进行N点FFT
Sk6=abs（Xk2）.^2./N;%取频谱幅度的平方，并除以N，以此作为对xn真实功率谱的估计
Sk7=（Sk5+Sk6）/2;%相加求平均
Sk8=10*log（Sk7）;
n=0:
255;
figure（4）
subplot（2,1,1）;
plot（n/N,Sk7）;
title（'（d）Bartlett法估计功率谱密度L=2'）
ylabel（'Sk'）
subplot（2,1,2）;
plot（n/N,Sk8）;
ylabel（'10log（PSD）'）
%L=4时bartlett法
N=128,n=0:
127;
x1n=2*cos（2*pi*0.2*n）+cos（2*pi*0.4*n）+2*randn（size（n））;%产生加有均值为0，方差为1的AWGN的信号的1/4段
Xk1=fft（x1n,N）;%进行N点FFT
Sk1=abs（Xk1）.^2./N;%取频谱幅度的平方，并除以N，以此作为对xn真实功率谱的估计
n=128:
255;
x2n=2*cos（2*pi*0.2*n）+cos（2*pi*0.4*n）+2*randn（size（n））;%产生加有均值为0，方差为1的AWGN的信号的1/4段
Xk2=fft（x2n,N）;%进行N点FFT
Sk2=abs（Xk2）.*（abs（Xk2））./N;%取频谱幅度的平方，并除以N，以此作为对xn真实功率谱的估计
N=128,n=256:
383;
x3n=2*cos（2*pi*0.2*n）+cos（2*pi*0.4*n）+2*randn（size（n））;%产生加有均值为0，方差为1的AWGN的信号的1/4段
Xk3=fft（x3n,N）;%进行N点FFT
Sk3=abs（Xk3）.^2./N;%取频谱幅度的平方，并除以N，以此作为对xn真实功率谱的估计
n=384:
511;
x4n=2*cos（2*pi*0.2*n）+cos（2*pi*0.4*n）+2*randn（size（n））;%产生加有均值为0，方差为1的AWGN的信号的1/4段
Xk4=fft（x4n,N）;%进行N点FFT
Sk4=abs（Xk4）.*（abs（Xk4））./N;%取频谱幅度的平方，并除以N，以此作为对xn真实功率谱的估计
Sk5=（Sk1+Sk2+Sk3+Sk4）/4;%相加求平均
Sk6=10*log（Sk5）;
n=0:
127;
figure（5）
subplot（2,1,1）;
plot（n/N,Sk5）;
title（'（e）Bartlett法估计功率谱密度L=4'）
ylabel（'Sk'）
subplot（2,1,2）;
plot（n/N,Sk6）;
ylabel（'10log（PSD）'）
（四）结果分析
下面是程序运行后的结果
上图分别是L=2和L=4时用bartlett法进行信号功率谱估计的波形。
从上图中我们同样可以看到信号中有两个频率分量，一个在0.2处，一个在0.4处，与产生的信号相一致。
但是我们不难看出，估计出的功率谱谱线与之前相关法功率谱估计和周期图法功率谱估计相比，波形相对平坦了一些。
L=2和L=4时用bartlett法进行信号功率谱估计的波形相比我们可以很明显的看出L大的那个波形更加平坦，这与之前在算法原理中介绍的一样，L越大，平均后的方差就越小，越能达到一致估计的目的。
五、Welch法功率谱估计
（一）算法原理简介
现在比较常用的功率谱估计改进方法是Welch法，又叫加权交叠平均法。
这种方法以加窗（加权）求取平滑，以分段重叠求得平均，因此集平均与平滑的优点于一体，同时也不可避免的带有两者的缺点。
其主要步骤如下：
第一步：
将N长的数据段分成L个小段，每小段M点，相邻小段见交叠M/2点。
第二步：
对个小段加同样的品挂窗后求傅里叶变换
第三步：
求个小段功率谱的平均，得
这里
（二）运算简要框图
矩形窗截断
X（n）分成L小段输出
（三）程序示例
N=128;n=0:
127;
x1n=2*cos（2*pi*0.2*n）+cos（2*pi*0.4*n）+2*randn（size（n））;%产生加有均值为0，方差为1的AWGN的信号一部分
wn=hanning（128）;
L=7;U=61.7942;%U是128长窗函数的能量,那个函数不会写，这个数是自己加出来的。
x11n=x1n.*wn';
Xk1=fft（x11n,128）;
Sk1=abs（Xk1）.^2.*1/U;
N=128;n=64:
191;
x2n=2*cos（2*pi*0.2*n）+cos（2*pi*0.4*n）+2*randn（size（n））;
x22n=x2n.*wn';
Xk2=fft（x22n,128）;
Sk2=abs（Xk2）.^2.*（1/U）;
N=128,n=128:
255;
x3n=2*cos（2*pi*0.2*n）+cos（2*pi*0.4*n）+2*randn（size（n））;
x33n=x3n.*wn';
Xk3=fft（x33n,128）;
Sk3=abs（Xk3）.^2.*（1/U）;
N=128,n=192:
319;
x4n=2*cos（2*pi*0.2*n）+cos（2*pi*0.4*n）+2*randn（size（n））;
x44n=x4n.*wn';
Xk4=fft（x44n,128）;
Sk4=abs（Xk4）.^2.*（1/U）;
N=128,n=256:
383;
x5n=2*cos（2*pi*0.2*n）+cos（2*pi*0.4*n）+2*randn（size（n））;
x55n=x5n.*wn';
Xk5=fft（x55n,128）;
Sk5=abs（Xk5）.^2.*（1/U）;
N=128,n=320:
447;
x6n=2*cos（2*pi*0.2*n）+cos（2*pi*0.4*n）+2*randn（size（n））;
x66n=x6n.*wn';
Xk6=fft（x66n,128）;
Sk6=abs（Xk6）.^2.*（1/U）;
N=128,n=384:
511;
x7n=2*cos（2*pi*0.2*n）+cos（2*pi*0.4*n）+2*randn（size（n））;
x77n=x7n.*wn';
Xk7=fft（x77n,128）;
Sk7=abs（Xk7）.^2.*（1/U）;
n=0:
127;
Sk=（Sk1+Sk2+Sk3+Sk4+Sk5+Sk6+Sk7）./L;
figure（6）
title（'（f）Welch法（汉宁窗，L=7,64点交叠）'）
subplot（2,1,1）;
plot（n/N,Sk）;
ylabel（'Sk'）
Skk=10*log（Sk）;
subplot（2,1,2）;
plot（n/N,Skk）;
ylabel（'10log（PSD）'）
（四）结果分析
下面是程序运行后的结果
上图分别是L=7，64点交叠用welch法进行信号功率谱估计的波形。
从上图中我们同样可以看到信号中有两个频率分量，一个在0.2处，一个在0.4处，与产生的信号相一致。
并且我们不难看出，估计出的功率谱谱线与之前上午功率谱估计相比，波形是最平坦的。
可见这种改进方法是比较准确的。
六、各波形直观比较
由以上六个图可以很清晰的看出不同功率谱估计的方法的优劣。
升级会员 