10-21：声音的STFT分析

再次将学习到的知识整理下来，来源有网络和课件。

1.平稳信号和非平稳信号

自然界中几乎所有信号都是非平稳信号，比如我们的语音信号就是典型的非平稳信号。
平稳信号在不同时间得到的采样值的统计特性(比如期望、方差等)是相同的，非平稳信号则与之相反，其特性会随时间变化。在信号处理中，这个特性通常指频率。
通常傅里叶变换只适合处理平稳信号，对于非平稳信号，由于频率特性会随时间变化，为了捕获这一时变特性，我们需要对信号进行时频分析，就包括短时傅里叶变换、小波变换、希尔伯特变换、希尔伯特黄变换这几种变换。

2.FT的局限

频谱分析，发现一些时域中看不到的信息，比如信号的频率分量组成、信号能量在频域上的分布等。
傅里叶变换是一种全局的变换，时域信号经过傅里叶变换后，就变成了频域信号，从频域无法看到时域信息，反变换同理。
【FT对非平稳信号和在分辨率上的局限性】无限长时间可以获得准确的频率定位（频率分辨率）；反之无限长的频率可以获得准确的时间定位（时间分辨力）；但是我们无法同时获得准确的时间定位和频率定位！
如果信号的频率特性在任何时间都不发生改变(即该信号是平稳信号)的话，使用傅里叶变换是没有问题的，然而如果该信号是非平稳信号，这时候时域信息就相当重要了。会出现不同频率变化但幅度谱一样的情况，我们希望能知道频率是怎么随时间变化的。

3.STFT

定义：
$X(n,w)=\sum^{ \infty}_{m = -\infty}x(m)g(n-m)e^{-jwn}$
其中， $x(m)$ 为输入信号， $g(m)$ 为窗函数， $X(n,w)$ 是时间 $n$ 和频率 $w$ 的二维函数，它将信号的时域和频域联系起来，我们可以据此对信号进行时频分析。其中， $S(n,w)=|X(n,w)|^2$ 就是语音信号所谓的语谱图(Spectrogram)。通过STFT，我们可以很容易地得出非平稳信号的时变特性。
计算语谱时采用不同窗长度，可以得到两种语谱图：

窄带语谱图。长时窗（至少两个基音周期）常被用于计算窄带语谱图。窄带语谱图具有较高的频率分辨率和较低的时间分辨率，良好的频率分辨率可以让语音的每个谐波分量更容易被辨别，在语谱图上显示为水平条纹。
宽带语谱图。短窗用于计算宽带语谱图。宽带语谱图具有较高的时间分辨率和较低的频率分辨率，低频率分辨率只能得到谱包络，良好的时间分辨率适合用于分析和检验语音的发音。

接下来开始利用上节学到的知识（时域分析和STFT）对声音进行分析。

4.代码实践

def load_data(name):
    '''
    name:文件路劲，字符串类型
    '''
    mat_contents = sio.loadmat(name)
    sig = mat_contents['sig']   # 
    fs = mat_contents['fs']     # 48000
    sig = sig[:,0]              # 选择第一个通道
    # https://www.osgeo.cn/numpy/reference/generated/numpy.asscalar.html
    fs = np.asscalar(fs)        # 将大小为1的数组转换为其标量等效值
    sigLen = sig.size           # 信号点数48000
    t = np.arange(0,sigLen)/fs  # 采样的时刻
    return sig,fs,sigLen,t

def FenDuanloc_pow(segL,overlap,sig):
    #信号分段
    #计算需要多少段 
    N = len(sig)
    delta = segL-overlap
    
    # 这里算需要多少段，(N-overlap)/(M-overlap ),M表示段长
    segNum = np.int32(np.ceil((N-overlap)/delta));
    
    #扩展信号：看最后有没有多出来一点，补0处理
    padNum = segNum*delta+overlap-N
    if padNum==0:
        sigEx = sig
    elif padNum>0:
        sigEx = np.hstack((sig,np.zeros(padNum)))    
    
    #分段标签:其实就是找每一段的起始位置
    segIdx = np.arange(0,segNum)*delta
    #生成分段矩阵
    segMat = as_strided(sigEx,shape=(segNum,segL),strides=(sigEx.strides[0]*delta,sigEx.strides[0]))
    loc_pow = np.sum(segMat**2,axis=1)/segL
    loc_max = np.max(segMat,axis=1)
    
    return segIdx,segMat,loc_pow,loc_max

sig,fs,sigLen,t = load_data('sound6')
plt.figure(figsize=(12,4))
plt.plot(t,sig)

原始信号.png

1.分析有哪些音

sig,fs,sigLen,t = load_data('sound6')
ff,P,nfft = Myfft(sig,fs,48000)
# 画出0-1000Hz的FFT结果
pIdx = findPink(P[:1000],-70)
plt.figure(figsize=(12,4))
plt.scatter(ff[pIdx],P[:1000][pIdx],c='r')
plt.plot(ff[:1000],P[:1000])
plt.title('1s信号FFT在0-1000Hz的值')
plt.xlabel('Hz')
plt.ylabel('dB')
plt.show()
print('峰值的频率为：',ff[pIdx])

找峰值.png

峰值的频率为： [348. 389.]

从结果来看，有F4音和G4音。

2.分析这些音在什么时刻出现（STFT分析）

sig,fs,sigLen,t = load_data('sound6')

NFFT = 4096              # 每一段信号的长度
overlapSize = NFFT/3     # 每一段可以重叠的长度
plt.figure(figsize=(12,4))
spectrum,freqs,ts,fig = plt.specgram(sig,NFFT = NFFT,Fs =fs,window=np.hanning(M = NFFT),noverlap=overlapSize,mode='default',scale_by_freq=True,sides='default',scale='dB',xextent=None)#绘制频谱图
plt.show()

时频图.png

plt.figure(figsize=(12,4))
plt.contourf(ts, freqs[30:45], spectrum[30:45,:],200)
plt.xlabel('时间/s')
plt.ylabel('频率/Hz')
plt.show()

时频图.png

可以看到大约在1.1s出现了440Hz的音（A4）和在4.4s左右有一个392Hz的音（G4）。但是上面找到了F4音，扩大氛围查看:

plt.figure(figsize=(12,4))
plt.contourf(ts, freqs[25:45], spectrum[25:45,:],200)

时频图.png

可以看到，在刚开始左右有一个349.23Hz的音（F4）。

3.用自己写的STFT进行时频分析

sig,fs,sigLen,t = load_data('sound6')

segL = 1024
overlap = int(segL/3)
segIdx,segMat,loc_pow,loc_max = FenDuanloc_pow(segL,overlap,sig)

ts = t[segIdx]

ff = np.linspace(0,sigLen,segL)*(fs/sigLen)

plt.figure(figsize=(12,4))
pMat = np.fft.fft(segMat)
pMat = np.abs(pMat)**2
plt.contourf(ts, ff[5:15], pMat[:,5:15].T,50)
plt.xlabel('时间/s')
plt.ylabel('频率/Hz')
plt.show()

时频图.png

可以看到，得到基本一致的结果。区别在于我取的每一段长度更短，频率分辨率没那么高。