2020-01-18 python实现stft并绘制时频谱

Python016

2020-01-18 python实现stft并绘制时频谱,第1张

官方文档中给出了非常详细的安装方法

http://librosa.github.io/librosa/install.html

函数声明:

librosa.core.stft(y, n_fft=2048, hop_length=None, win_length=None, window='hann', center=True, dtype=<class 'numpy.complex64'>, pad_mode='reflect')

常用参数说明:

y:输入的numpy数组,要求都是实数

n_fft:fft的长度,默认2048

hop_length:stft中窗函数每次步进的单位

win_length:窗函数的长度

window:窗函数的类型

return:一个1+n_fft/2*1+len(y)/hop_length的二维复数矩阵,其实就是时频谱

参考:

http://librosa.github.io/librosa/generated/librosa.core.stft.html#librosa.core.stft

主要用这两个

matplotlib.pyplot.pcolormesh()

matplotlib.pyplot.colorbar()

二维FFT常用在图像处理上,首先要能理解二维FFT的意义,否则很难明白它到底是怎么工作的。

第一列是原图和对应的频率信息,第二列是去除低频部分后,FFT逆变换得到的图像。第三列是去除高频部分后FFT逆变换得到的图像。

从第二列可以看出高频贡献了图像的细节。从白到黑的边界保留了下来。而原图中大片的白与大片的黑在这个图中没什么区别。

第三列中保留了原图中的亮部与灰部,而由黑到白的临界线却很模糊。细小的白线黑线也没能显示。所以低频贡献了图像的明暗。

2.工作原理理解

二维FFT就是先对行做次一维FFT,这样每个元素都是关于行频率信息了,然后再对列做一维FFT,这样每个元素都包含了行和列的频率信息。每个元素都是个复数,取绝对值可得到振幅,从实部与虚部的比值可等到相位,在二维矩阵的位置信息包含了频率大小和方向。方向在一维FFT中是不用考虑的。

FFT2的结果也是正频率从0到高然后负频率从高到0.fftshift()之后会将低频放到中间位置。

第一幅图的频谱是中间一条白线,也就是说许多个正弦波沿横向传播。纵向上没有变化。

第三幅图的频谱是十字形加一条从左下角到右上角的直线。说明原图在横向,纵向都有变化,变化的方向从左下角到右上角。

从中心到频谱图上某一点构成的向量方向就是这个波传播的方向。

正负对称才能消除虚部,这点与一维FFT原理一致。

语音的时域分析和频域分析是语音分析的两种重要方法,但是都存在着局限性。时域分析对语音信号的频率特性没有直观的了解,频域特性中又没有语音信号随时间的变化关系。而语谱图综合了时域和频域的优点,明显的显示出了语音频谱随时间的变化情况、语谱图的横轴为时间,纵轴为频率,任意给定频率成分在给定时刻的强弱用颜色深浅来表示。颜色深的,频谱值大,颜色浅的,频谱值小。语谱图上不同的黑白程度形成不同的纹路,称之为声纹,不同讲话者的声纹是不一样的,可用作声纹识别。

下面是在python中绘制语谱图:

# 导入相应的包

import numpy, waveimport matplotlib.pyplot as pltimport numpy as npimport os

filename = 'bluesky3.wav'

# 调用wave模块中的open函数,打开语音文件。f = wave.open(filename,'rb')

# 得到语音参数

params = f.getparams()

nchannels, sampwidth, framerate,nframes = params[:4]

# 得到的数据是字符串,需要将其转成int型

strData = f.readframes(nframes)

wavaData = np.fromstring(strData,dtype=np.int16)

# 归一化

wavaData = wavaData * 1.0/max(abs(wavaData))

# .T 表示转置

wavaData = np.reshape(wavaData,[nframes,nchannels]).T

f.close()

# 绘制频谱

plt.specgram(wavaData[0],Fs = framerate,scale_by_freq=True,sides='default')

plt.ylabel('Frequency')

plt.xlabel('Time(s)')

plt.show()