再次将学习到的知识整理下来，来源有网络和课件。
1.平稳信号和非平稳信号
- 自然界中几乎所有信号都是非平稳信号，比如我们的语音信号就是典型的非平稳信号。
- 平稳信号在不同时间得到的采样值的统计特性(比如期望、方差等)是相同的，非平稳信号则与之相反，其特性会随时间变化。在信号处理中，这个特性通常指频率。
- 通常傅里叶变换只适合处理平稳信号，对于非平稳信号，由于频率特性会随时间变化，为了捕获这一时变特性，我们需要对信号进行时频分析，就包括短时傅里叶变换、小波变换、希尔伯特变换、希尔伯特黄变换这几种变换。
2.FT的局限
- 频谱分析，发现一些时域中看不到的信息，比如信号的频率分量组成、信号能量在频域上的分布等。
- 傅里叶变换是一种全局的变换，时域信号经过傅里叶变换后，就变成了频域信号，从频域无法看到时域信息，反变换同理。
- 【FT对非平稳信号和在分辨率上的局限性】无限长时间可以获得准确的频率定位（频率分辨率）；反之无限长的频率可以获得准确的时间定位（时间分辨力）；但是我们无法同时获得准确的时间定位和频率定位！
- 如果信号的频率特性在任何时间都不发生改变(即该信号是平稳信号)的话，使用傅里叶变换是没有问题的，然而如果该信号是非平稳信号，这时候时域信息就相当重要了。会出现不同频率变化但幅度谱一样的情况，我们希望能知道频率是怎么随时间变化的。
3.STFT
定义：
其中，为输入信号，
为窗函数，
是时间
和频率
的二维函数，它将信号的时域和频域联系起来，我们可以据此对信号进行时频分析。其中，
就是语音信号所谓的语谱图(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
可以看到，得到基本一致的结果。区别在于我取的每一段长度更短，频率分辨率没那么高。
网友评论