1、先导知识

1. 高等数学（微积分、线代、概率论及数理统计）
2. 信号与系统、数字信号处理
3. 随机信号分析

2、功率谱以及谱估计

• 功率谱的基本概念
• 谱估计的概念
• 自相关序列的估计

PAM的功率谱
$$P_s(f)=\frac{{\sigma }_a^2}{T_s}|G_T(f){|}^2+\frac{m_a^2}{T_s^2}\sum _{k=-\infty }^{+\infty }|G_T(k{R}_s){|}^2\delta (f-k{R}_s)$$

（1）功率谱的基本概念

我们首先有这样一个常识，对于一个确定的信号，它唯一对应了一个确定的频谱，这是信号与系统课程中我们学到的看世界的两种角度——时域、频域。

现在我们面对的是随机的信号，不知道信号在下一时刻的取值是1还是0，但是我们知道它出现1的概率，以及它出现0的概率。也就是说，不同时间观察到的信号是不一样的，也就是意味着观察到的频谱是不一样的。现在我们试图在变化中抓住本质，抓住那个不变的东西。

从随机信号分析这门课程里，我们学习到了自相关函数。对于一个广义平稳随机信号而言，它的均值为一个常数，自相关函数只与观察时差有关，而与观察时刻无关。与观察时刻无关就令人十分满意了，因为这意味着就算我明天来观察这个信号，今天所计算出来的自相关函数仍然是有用的。

但是自相关函数仍然不好看，因为它终归是从时域角度观察的。由维纳-辛钦定理我们知道，自相关函数的傅里叶变换就是这个信号的功率谱密度函数，常常简称其为功率谱。因此，我们得出个结论，广义平稳随机信号的功率谱是不变的，是一个能体现随机信号统计特性的统计量，其物理意义是单位频率的平均功率（单位是W/Hz）。

写成数学公式：
$$\begin{gathered} R_X(\tau )=E[X(t+\tau )X(t)] \\ {S}_X(j\omega )={\int }_{-\infty }^{+\infty }{R}_X(\tau )e^{-j\omega \tau }d\tau \end{gathered}$$

（2）谱估计的概念

从上述分析我们知道了功率谱是个好东西，但是现实是残酷的，因为我们不可能观察一个信号无穷时长，自然就得不到无穷宽的时差，也就得不到绝对精准的功率谱。

我们现在有的仅仅是有限时长的一段离散信号，这个时候就可以称其为序列了。但是我们试图在仅有的条件下，尽可能地反应现实世界。

在随机信号分析中，我们知道了各态历经这个概念，当观察时间长到足以把随机信号能经历的所有“状态”都经历一遍，那么咱就不观察了。我们就用此时得到的足够长的离散序列来求自相关序列，并通过快速傅里叶变换得到功率谱。

（3）自相关序列的估计

这部分主要涉及到自相关序列的无偏估计和有偏估计，关于其无偏性和有效性的数学推导较为复杂，但学过先导知识里面的课程之后，还是能够顺着参考书里面的步骤推导出来的。这里只列出结果。
$$\begin{gathered} u(0),u(1),...,u(N-1)\text{是观察到的样本序列。} \\ {r}_u(m)是自相关序列 \\ {\hat{r}}_u(m)是r_u(m)的有偏估计，\hat{r}{{}^{\prime }}_u(m)是r_u(m)的无偏估计 \end{gathered}$$

• 有偏估计

r ^ u ( m ) = 1 N ∑ n = 0 N − 1 − ∣ m ∣ u ( n ) u ∗ ( n + m ) , ∣ m ∣ ≤ N − 1 E [ r ^ u ( m ) ] = N − ∣ m ∣ N r u ( m ) V a r [ r ^ u ( m ) ] = ( N − ∣ m ∣ N ) 2 { ( 1 N − ∣ m ∣ ) 2 ∑ r = − ( N − 1 − ∣ m ∣ ) N − 1 − ∣ m ∣ ( N − ∣ m ∣ − ∣ r ∣ ) [ r u 2 ( r ) + r u ( r − m ) r u ( r + m ) ] } \hat r_u(m) = \frac{1}{N}\sum_{n=0}^{N-1-|m|}u(n)u^*(n+m),|m|\le N-1 \\ E[\hat r_u(m)] = \frac{N-|m|}{N} r_u(m) \\ Var[\hat r_u(m)] = (\frac{N-|m|}{N})^2 \{ (\frac{1}{N-|m|})^2 \sum_{r = -(N-1-|m|)}^{N-1-|m|}(N-|m|-|r|) [r^2_u(r)+r_u(r-m)r_u(r+m)] \} r^u​(m)=N1​n=0∑N−1−∣m∣​u(n)u∗(n+m),∣m∣≤N−1E[r^u​(m)]=NN−∣m∣​ru​(m)Var[r^u​(m)]=(NN−∣m∣​)2{(N−∣m∣1​)2r=−(N−1−∣m∣)∑N−1−∣m∣​(N−∣m∣−∣r∣)[ru2​(r)+ru​(r−m)ru​(r+m)]}

• 无偏估计
  $$\begin{gathered} {\hat{r}}_u^{\prime }(m)=\frac{1}{N-|m|}\sum _{N=0}^{N-1-|m|}u(n)u^{\ast }(n+m),|m|\le N-1 \\ E[{\hat{r}}_u^{\prime }(m)]=r_u(m) \\ Var[{\hat{r}}_u^{\prime }(m)]=(\frac{1}{N-|m|}{)}^2\sum _{r=-(N-1-|m|)}^{N-1-|m|}(N-|m|-|r|)[r_u^2(r)+r_u(r-m)r_u(r+m)] \end{gathered}$$

• 结论

有偏估计的方差和均方误差比无偏估计的小，并且不会带来负的功率谱估计值，因此通常使用有偏估计式来估计自相关序列。

3、谱估计的经典方法及matlab实现

所谓经典方法，就是通过有偏估计式${\hat{r}}_u(m)$进行傅里叶变换来计算出功率谱。

• 间接法
  • 平滑周期图法（Blackman-Turkey法）：加窗，减少频谱泄露
• 直接法
  • 周期图法：矩形窗，频谱分辨率低
  • 平均周期图法（Bartlett法）：分段，减少方差，降低分辨率
  • 平均修正周期图法（Welch法）：重叠分段，加归一化窗

（1）间接法

所谓间接法就是先求出自相关序列的有偏估计式，然后对有偏估计式进行傅里叶变换，由维纳辛钦定理得到功率谱的估计，这种方法最早由Blackman和Turkey两个人提出的，又叫BT法。

clf; f= 10;
N= 1000; Fs= 500;                   %数据长度和采样频率
n=0: N-1; 
t= n/Fs;                            %时间序列
Lag= 100;                           %延迟样本点数
x=sin(2* pi * f * t)+0.6* randn(1, length(t));   %原始信号
[c, lags]= xcorr(x, Lag, 'unbiased' );           %对原始信号进行无偏自相关估计
subplot(311),plot(t,x);                          %绘原始信号x
xlabel('t/s'); ylabel('x(t)'); grid;
legend('含噪声的信号x(t)');
subplot(312); plot(lags/Fs, c);                  %绘x信号自相关，lags/Fs为时间序列
xlabel('t/s'); ylabel('Rxx(t)');
legend('信号的自相关Rxx');grid;

Pxx= fft(c, length(lags));                       %利用FFT变换计算信号的功率谱
fp= (0:length(Pxx)- 1)'* Fs/length(Pxx);         %求功率谱的横坐标f
Pxmag= abs(Pxx) ;%求幅值
subplot(313);
plot(fp(1:(length(Pxx)-1)/2),Pxmag(1:(length(Pxx)-1)/2)); %绘制功率谐曲线
xlabel('f/Hz'); ylabel('|Pxx|');
grid;
legend('信号的功率谱Pxx');

这种方法通常需要对估计式加窗进行修正（这个时候其实也可以认为是在平滑周期图），减小自相关序列截断的影响（从无限到有限的截断）。加窗的本质就是在减少频谱泄露。

同时注意窗型和窗长：窗长和窗型同时影响分辨率，窗型主要影响谱估计的方差。

（2）直接法

直接法是指，先求出序列的离散傅里叶变换，然后求模平方，得到能量谱，再除以观察时间就得到功率谱，这种方法又称周期图法。为什么称为周期图（periodogram）法呢，这个和离散傅里叶变换的原理有关，离散傅里叶变换就是认为现在观察的时间已经足够长了，采样频率也足够高了，那么再去观察的话，不过是对现有观察序列的重复，那么也就是说现有序列是这个本该是无限长序列的一个周期。

显然这个假设在大多数时候是不符合实际的，而且谱分辨率也是低的，但是由于可以采用FFT算法，它的计算效率很高，在对谱分辨率要求不高的地方，这种方法还是很常用的。不过，后面还有对周期图法的改进方案。

%% 直接法
clf; f1= 60;f2 = 120;
N= 1000; Fs= 500;                   %数据长度和采样频率
x=sin(2* pi * f1 * t)+ sin(2* pi * f2 * t)+ 0.6* randn(1, length(t));   %原始信号
window = boxcar(length(x));
[Pxx,fx]=periodogram(x,window,N,Fs);
plot(fx,10*log10(Pxx));
% Pxx1 = abs(fft(x)).^2 /(N*Fs); %直接计算，跟peridogram的结果近似

（3）Bartlett法

这种方法的思路很简单，就是将原本的序列切割成两段，每段各自按照periodogram的方法求一个功率谱出来，然后把两个功率谱加起来除以二，从而将缩小原始的周期图法的方差。显然，可以切割成多段。但是这种切割是要付出代价的，代价就是谱分辨率降低。为什么会降低呢？因为原序列N个点，切割成两段，每段就是N/2个点，在采样率相同的情况下，N点FFT当然比N/2点FFT的分辨率高了（再深入的话就偏题了）。

%% Bartlett法
clf; f1= 60;f2 = 120;
N= 1000; Fs= 500;                   %数据长度和采样频率
n=0: N-1; 
x=sin(2* pi * f1 * t)+ sin(2* pi * f2 * t)+ 0.6* randn(1, length(t));   %原始信号
nfft = N/2;%切成两截
window = boxcar(length(x));
[Pxx2,fx]=periodogram(x,window,nfft,Fs);
plot(fx,10*log10(Pxx2));


% Pxx_2 = abs(fft(x(1:nfft),nfft)).^2 /(nfft*Fs)+abs(fft(x(nfft:2*nfft),nfft)).^2 /(nfft*Fs);%直接计算，结果相近
% plot(fx,10*log10(Pxx2),'r-',fx,10*log10(Pxx_2(1:nfft/2+1)),'b--');

（4）Welch法

Welch法是结合了Bartlett法和Blackman-Turkey法，并加以改进的一种方法。第一个改进的地方是切割成重叠的分段，用房地产行业的术语来说就是公摊面积。原本N个点分成两截，每截N/2，现在一人拿出N/4来公摊，那就是每截0.75N了，不过有0.5N是重复的，也就是帧长0.75N，帧移只有0.25N，重叠了0.5N。这样在一定程度上解决了谱分辨率和谱估计方差的矛盾，但有重叠，那么互相关性就会增强嘛，也就是说谱估计方差没有Bartlett法那么理想。

同时，也采用了BT法的思路，给每段序列加窗，改进的地方在于加了一个归一化因子——窗函数的能量，这样使得任何窗函数都可以得到非负的谱估计。

clf; f1= 60;f2 = 120;
N= 1000; Fs= 500;                   %数据长度和采样频率
n=0: N-1; 
x=sin(2* pi * f1 * t)+ sin(2* pi * f2 * t)+ 0.6* randn(1, length(t));   %原始信号
nfft = N/2;%切成两截，它同时也是帧长
overlap = nfft*0.5;
sa = nfft-overlap;%帧移
window = hamming(nfft);
U = sum(window.^2)/length(window);
[Pxx3,fx]=pwelch(x,window,overlap,nfft,Fs);
plot(fx,10*log10(Pxx3));
xlabel("f/Hz")
ylabel("dB")

% Pxx_3 = 1/U*(abs(fft(window' .* x(1:nfft),nfft)).^2 /(nfft*Fs)...
%     +abs(fft(window' .* x((1:nfft)+sa),nfft)).^2 /(nfft*Fs)...
%     +abs(fft(window' .* x((1:nfft)+2*sa),nfft)).^2 /(nfft*Fs));
% 
% plot(fx,10*log10(Pxx3),'r-',fx,10*log10(Pxx_3(1:nfft/2+1)),'b--');

4、参考书目

1、《通信原理（第二版）》李晓峰等著。

2、《谱估计与自适应信号处理教程》赵晓晖编著。

3、《MATLAB辅助现代工程数字信号处理》李益华主编。