10-11:功率谱估计及频率分辨率

1.频率相关的概念

对于周期信号，信号的频率等于 $1/T$ ， $T$ 为周期，就是完成往复运动一次所需的时间。频率是单位时间内某事件重复发生的次数，频率 $f=1/T$ 赫兹。 $T$ 大，说明相同时间内事件重复发生的次数少，即频率 $f$ 小。那么在DFT中，信号的频率是什么呢？
假如我要开始对脉搏信号进行采集，我每隔0.02s（即时域采集间隔 $t_s$ ）采集一次，一共采集了 $N$ （采样点数）个点，那么我采集这 $N$ 一共花了 $T_0=N*t_s$ 长的时间，我们认为这一段信号的时域长度为 $T_0$ 。那么这一段信号的频率上限为 $f_s=1/t_s$ （即采样频率）。根据采样定理，采样频率要大于信号频率的两倍。

2.FFT频谱分析

一直对FFT理解不到位，学习相关材料后梳理一下。
一个模拟信号，经过ADC采样之后，就变成了数字信号。 $N$ 个采样点，经过FFT之后，就可以得到 $N$ 个点的FFT结果（ $N$ 常取2的整数次方）。假设采样频率为 $f_s$ ，信号频率 $F$ ，采样点数为 $N$ 。那么FFT之后结果就是一个为 $N$ 点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。
【结论】假设原始信号的峰值为 $A$ ，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是 $A$ 的 $N/2$ 倍。而第一个点就是直流分量，它的模值就是直流分量的 $N$ 倍。而每个点的相位呢，就是在该频率下的信号的相位。第一个点表示直流分量（即 $0Hz$ ），而最后一个点则表示采样频率 $f_s$ ，这中间被 $N-1$ 个点平均分成 $N$ 等份，每个点的频率依次增加，点 $n(n=0,1,...,N)$ 所表示的频率为： $F_n=n*f_s/N$ 。由此可见， $F_n$ 所能分辨到频率为 $f_s/N$ 。
【根据FFT结果反推原信号的幅值、频率和相位】假设FFT之后某点 $n$ 用复数 $a+bi$ 表示，那么这个复数的模就是 $A_n=\sqrt[]{a^2+b^2}$ ，相位就是 $P_n=arctan(b/a)$ 。根据以上的结果，就可以计算出 $n$ 点（ $n≠0$ ，且 $n<N/2$ ）对应的信号的表达式为： $A_n/(N/2)*cos(2*pi*F_n*t+P_n)$ 即 $2*A_n/N*cos(2*pi*F_n*t+P_n)$ 。对于 $n=1$ 点的信号，是直流分量，幅度即为 $A_1/N$ 。
由于FFT结果的对称性，通常我们只使用前半部分的结果，即小于采样频率一半的结果。
【动手实践】假设我们有一个信号，它含有2V的直流分量和两个交流分量，一个是频率为50Hz、相位为-30度、幅度为3V的交流信号，另一个是一个频率为75Hz、相位为90度、幅度为1.5V的交流信号。
$S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)$ 式中cos参数为弧度，所以-30度和90度要分别换算成弧度。

原始信号.png

采样频率设置为256HZ，采样点数为256，即

f_s=256,N=256

根据上述的结论

F_n=n*f_s/N

，两个点之间的间距就是1Hz，其即为频谱图的横坐标

F

(频率)。我们的信号有3个频率：0Hz、50Hz、75Hz，应该分别频谱图的第0个点、第50个点、第75个点上出现峰值，其它各点应该接近0。

FFT取模.png

上图横坐标是采样点，纵坐标是FFT结果取模。我们可以看到除第0个点外上图是对称的，且在第50个点、第75个点取到峰值，这是因为刚好点与点之间的间隔是1Hz的缘故。第0个点、第50个点和第75个点的幅度值分别为：512，384，192。其余地方的幅度值接近0。

FFT换算后幅频图.png

此时，经过转换，上图横坐标是频率，纵坐标是幅值。按照公式，可以计算出直流分量为：

512/N=512/256=2

；50Hz信号的幅度为：

384/(N/2)=384/(256/2)=3

；75Hz信号的幅度为

192/(N/2)=192/(256/2)=1.5

。与原先设置的信号幅值相同。

FFT换算后相频图.png

直流信号没有相位，第50个点和第75个点经过

arctan(b/a)*180/pi

转换为角度后分别为-30度和90度，与原先设定一致。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['font.sans-serif'] = ['SimHei']  
plt.rcParams['axes.unicode_minus'] = False  

N = 256     # 采样点数
fs = 256    # 采样频率
ts = 1/fs   # 采样间隔，意味着0.001s采集一次
T = ts*N    # 信号时间长度，也会等于N/fs
t = np.arange(0,T,ts)  # 采样时刻
Signal = 2+3*np.cos(2*np.pi*50*t-np.pi*30/180)+1.5*np.cos(2*np.pi*75*t+np.pi*90/180) 

plt.figure(figsize=(15,6))
plt.plot(Signal)
plt.title('原始信号')
plt.xlabel('采样点')
plt.ylabel('信号值')
plt.savefig('原始信号.png')
plt.show()

Y = np.fft.fft(Signal,N)
AY = np.abs(Y)       # 取模

plt.figure(figsize=(15,6))
plt.stem(AY)
plt.title('FFT取模')
plt.xlabel('采样点')
plt.ylabel('模值')
plt.savefig('FFT取模.png')
plt.show()

AY_true = AY/(N/2)      # 换算成实际的幅度
AY_true[0] = AY[0]/N    # 直流分量
F = np.arange(N)*fs/N   # 换算成实际的频率

plt.figure(figsize=(15,6))
plt.stem(F[:N//2],AY_true[:N//2])
plt.title('FFT换算后幅频图')
plt.xlabel('频率/Hz')
plt.ylabel('幅度值')
plt.savefig('FFT换算后幅频图.png')
plt.show()

PY = np.angle(Y)
PY_angle = PY*180/np.pi

plt.figure(figsize=(15,6))
plt.stem(F[:N//2],PY_angle[:N//2])
plt.title('FFT换算后相频图')
plt.xlabel('频率/Hz')
plt.ylabel('相位[-180,180]')
plt.savefig('FFT换算后相频图.png')
plt.show()

【总结】假设采样频率为 $f_s$ ，采样点数为 $N$ ，做FFT之后，某一点 $n$ （ $n$ 从0开始）表示的频率为： $F_n=n*f_s/N$ ；该点的模值除以 $N/2$ 就是对应该频率下的信号的幅度（对于直流信号是除以 $N$ ）；该点的相位即是对应该频率下的信号的相位。要精确到 $xHz$ ，则需要采样长度为 $1/x$ 秒的信号，并做FFT。要提高频率分辨率，就需要增加采样点数。比较简单的方法是采样比较短时间的信号，然后在后面补充一定数量的0，使其长度达到需要的点数，再做FFT，这在一定程度上能够提高频率分辨力。

3.频率分辨率

【第一种理解】频率分辨率的定义是DFT频域相邻刻度之间的实际频率之差。设 $f_0$ 表示频率分辨率， $f_0=f_s/N=1/(N*ts)=1/T_0$ ，就相当于把 $f_s$ 分成 $N$ 等分。
如果采样频率 $f_s$ 为1024Hz，采样点数为1024点，则可以分辨到 $1024/1024=1Hz$ 。1024Hz的采样率采样1024点，刚好是 $1024*(1/1024)=1（即N*t_s）$ 秒，也就是说，在这个情况下采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时间的信号并做FFT，则结果可以分析到0.5Hz。如果要提高频率分辨力，则必须增加采样点数，也即采样时间 $T_0$ 。在采样频率一定情况下，频率分辨率和采样时间是倒数关系即 $f_0=1/T_0$ 。在采样频率一定情况下，采样点数的多少与要求多大的频率分辨率有关，同时还要考虑频率畸形和信号截断而产生泄露的问题。
【第二种理解】频率分辨率是指所用的算法（如功率谱估计）能将信号中两个靠得很近的谱峰保持分开的能力，是用来比较和检验不同算法性能好坏的指标。

4.泄露与窗函数

（相关资料）每次FFT变换只能对有限长度的时域数据进，行变换，因此，需要对时域信号进行信号截断。信号截断有两种，一种是周期截断，一种是非周期截断，哪怕原始信号是周期信号。若周期截断，则FFT频谱为单一谱线，得到的频率成分为原始信号的真实频率，并且幅值与原始信号的幅值相等。
若为非周期截断，截断后的信号起始时刻和结束时刻的幅值不等，将这个信号再进行重构，在连接处信号的幅值不连续，出现跳跃；对截断后的信号做FFT，频谱出现拖尾，峰值处的频率与原始信号的频率相近，但并不相等。另一方面，峰值处的幅值已不再等于原始信号的幅值，为原始信号幅值的64%（矩形窗的影响）。而幅值的其他部分（36%幅值）则分布在整个频带的其他谱线上。拖尾现象这种非常严重的误差，称为泄漏，是数字信号处理所遭遇的最严重误差。现实世界中，在做FFT分析时，很难保证截断的信号为周期信号，因此，泄漏不可避免。为了将这个泄漏误差减少到最小程度（注意是减少，而不是消除），我们需要使用加权函数，也叫窗函数。加窗主要是为了使时域信号似乎更好地满足FFT处理的周期性要求，减少泄漏。非周期截断的信号与窗函数相乘得到的信号起始点与最末点达到相同（比如都为0），变成一个类似周期截断的信号。窗函数只能减少泄漏，不能消除泄漏。

5.功率谱估计

上节关于功率谱估计的部分有一点错误，我没有理解分段加窗函数的含义。分段之后，需要补0至和原先信号一样长，这就相当于每小段加矩形窗，矩形窗不仅有1值还有0值。在连续的世界里非常理所当然的事情，在离散的世界里就发生了变化。如果不补0到原来长度，频率分辨率就大大降低了。
这里我重新用经典的四种方法（直接法、间接法、Bartlett法和Welch法）来处理一个随机信号。信号表示如下： $signal= cos(2*pi*1000*t)+3*cos(2*pi*2000*t)+noise$ 采样频率是3倍的最大频率，采样点数是1024点。四种功率谱估计的结果如下图所示：

四种功率谱估计.png

实验的代码如下，代码较长，没有对代码进行重构。可以看到在2000Hz的地方出现了峰。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False

def myAutocorrelation(SigA):
    '''
    @author:zengwei
    自相关函数,属于有偏估计，对应matlab的xcorr(SigA,'biased')
    https://ww2.mathworks.cn/help/matlab/ref/xcorr.html
    '''
    N = len(SigA)
    SigB = np.zeros(2*N-1)
    SigB[0:N] = SigA
    L = len(SigB)
    
    SigC = np.array(SigB.tolist()*(L+1))
    Matrix = np.zeros((L,L))
    start = N-1
    for i in np.arange(L):
        end = start + L               # 还可以用np.roll()函数
        Matrix[i,:] = SigC[start:end]
        start = end - 1
    return np.dot(Matrix,SigB)/N

def direct(signal):
    '''
    直接法（周期图法）求功率谱估计，并做归一化
    '''
    fft_signal = np.fft.fft(signal)
    P = (np.abs(fft_signal)**2)/len(signal)  # FFT取模的平方再除于N
    return P/np.max(P)

def indirect(signal):
    '''
    间接法（自相关法）求功率谱估计，并做归一化
    '''
    selfcorr = myAutocorrelation(signal) # 长度变成了2N-1
    fft_selfcorr = np.fft.fft(selfcorr)
    P = np.abs(fft_selfcorr)/len(selfcorr)
    return P/np.max(P)

def Bartlett(Signal,L):
    '''
    L:分成多少段
    '''
    N = len(Signal)
    M = N//L                # 每小段的长度
    signal = Signal.reshape(L,M)
    P = []
    for i in signal:
        everyParagraph = np.hstack((i,np.zeros(N-M)))  # 每一段要补0至原长度N
        selfcorr = np.fft.fft(everyParagraph)                    # 每一段做FFT,长度为N
        P_i = np.abs(selfcorr)**2                                # FFT取模的平方,长度为N
        P.append(P_i.tolist())
    P_ = np.sum(np.array(P),axis=0)/L                                   # L段对应元素相加;按列求和
    return P_/np.max(P_)

def Hamming(N):
    return np.array([0.54 - 0.46 * np.cos(2 * np.pi * n / (N - 1)) for n in range(N)])
def Hanning(N):
    return np.array([0.5 - 0.5 * np.cos(2 * np.pi * n / (N - 1)) for n in range(N)])
def Rect(N):
    return np.ones(N)

def Welch(Signal,M,s,window=None):
    '''
    M：每一段的长度
    s：可以重叠的长度
    '''
    N = len(Signal)
    L = int((N-s)//(M-s))
    signal = np.zeros((L,N))                      # 保证每一段跟原先一样长
    for m in np.arange(L):
        signal[m,:M] = Signal[int(m/2):int(m/2)+M]
        
    if window == 'Hanning':
        d2 = Hanning(M)
        U = np.sum(d2**2)/M
        d2 = np.hstack((d2,np.zeros(N-M)))
    elif window == 'Hamming':
        d2 = Hamming(M)
        U = np.sum(d2**2)/M
        d2 = np.hstack((d2,np.zeros(N-M)))
    else:
        d2 = 1
        U = 1
    
    P = []
    for i in signal:
        selfcor = np.fft.fft(i*d2)
        P_i = np.abs(selfcor)**2/(M*U)
        P.append(P_i.tolist())
    P_ = np.sum(np.array(P),axis=0)/L
    return P_/np.max(P_)

f0 = 1000                 # 信号频率1
f1 = 2000                 # 信号频率2

fs = 3*f1                # 采样频率，要满足采样定理
N = 1024                  # 采样点数
ts = 1/fs                # 采样间隔
T0 = N*ts                # 信号时间长度，采样时间
t = np.arange(0, T0, ts) # 采样时刻
f0 = fs/N                # 频率分辨率
F = np.arange(N)*f0      # 频率
"""
生成原始信号序列，在原始信号中加上噪声
"""
x = np.cos(2*np.pi*f0*t) + 3*np.cos(2*np.pi*f1*t) + np.random.randn(t.size)
#x = np.sin(2*np.pi*f0*t)+np.sin(2*np.pi*f1*t) + 10*(-SNR/20)*np.random.randn(N)

plt.figure(figsize=(20, 20))
ax=plt.subplot(321)
plt.plot(t,x)
ax.set_title('原始输入信号')
plt.xlabel('时间t')
plt.ylabel('信号值x')
plt.tight_layout()

"""
FFT变换
"""
fftx = np.fft.fft(x)
fft_x = np.abs(fftx)/(N/2)
fft_x[0] = fftx[0]/N

ax=plt.subplot(322)
plt.plot(F[:N//2],fft_x[:N//2])
ax.set_title('FFT幅频图')
plt.xlabel('频率/Hz')
plt.ylabel('幅值')
plt.tight_layout()

"""
直接法（周期图法）功率谱估计，并归一化
"""
P_direct = direct(x)
P_direct = 10*np.log10(P_direct)

ax=plt.subplot(323)
plt.plot(F[:N//2],P_direct[:N//2])
ax.set_title('直接法功率谱估计')
plt.xlabel('频率/Hz')
plt.ylabel('归一化功率谱P(k)/dB')
plt.tight_layout()

"""
间接法（自相关法）功率谱估计，并归一化
"""
P_indirect = indirect(x)
P_indirect = 10*np.log10(P_indirect)
F_indirect = np.arange(2*N-1)*fs/(2*N-1)

ax=plt.subplot(325)
#plt.plot(F_indirect[:N//2],P_indirect[:N//2])
plt.plot(F_indirect[:N],P_indirect[:N])
#plt.xlim((0, N//2))
ax.set_title('间接法功率谱估计')
plt.xlabel('频率/Hz')
plt.ylabel('归一化功率谱P(k)/dB')
plt.tight_layout()

"""
Bartlett法功率谱估计，并归一化
"""
L = 64    # 分成64段
P_Bartlett = Bartlett(x,L)
P_Bartlett = 10*np.log10(P_Bartlett)

ax=plt.subplot(324)
plt.plot(F[:N//2],P_Bartlett[:N//2])
ax.set_title('Bartlett法功率谱估计')
plt.xlabel('频率/Hz')
plt.ylabel('归一化功率谱P(k)/dB')
plt.tight_layout()

"""
Welch法功率谱估计，并归一化
"""
M = 32        # 每段长度
s = M//2      # 可以重叠的长度
P_Welch = Welch(x,M,s)
P_Welch = 10*np.log10(P_Welch)

ax=plt.subplot(326)
plt.plot(F[:N//2],P_Welch[:N//2])
ax.set_title('Welch法功率谱估计')
plt.xlabel('频率/Hz')
plt.ylabel('归一化功率谱P(k)/dB')
plt.tight_layout()
plt.savefig('四种功率谱估计.png')
plt.show()

同时，试了一下对Welch法换了一下汉明窗和汉宁窗看看效果。如下所示。可以看到比矩形窗要更平滑一些了。

两种窗的Welch功率谱估计.png

6.功率谱估计的分辨能力

如何能分辨两个很近的峰，如何能准确表征一个峰的频率？
频率分辨率 $f_0=f_s/N=1/(N*ts)=1/T_0$ ， $f_0$ 是频率之间间隔，越小分辨频率的能力越强。理论上在满足采样定理的情况下，减小采样频率，增大采样点数可以增大频率分辨能力。这些都做不到的情况下，还可以在后面补0类似于增大采样点数。
我实验测试之后发现，减小采样频率增大采样点数，确实比之前分辨能力更强了。在较大采样频率和较小采样点数的情况下，相近的两个频率混在一个峰上了。
【10月12日批注】今天请教老师后，原来频率分辨率只是跟信号时间长度有关。我在实验时同时减小采样频率增大采样点数才误以为减小采样频率也会有用。Amazing！