python在信号处理的应用PDF python做信号频率分析

文章目录

• 题目
• 完整代码

---

题目

录制一段自己的语音（单通道），编写Python代码，分析其采样频率、数据长度。画出语音信号的时域波形，提取出其中清音和浊音的帧（帧长为N），画出帧的时域波形，并加窗，分析其频谱。 描述清音和浊音在时域、频域的差别。 改变帧长、窗函数，观察不同帧长和窗函数对频谱带来的影响。
注意：本题使用的是单通道音频，若使用双通道下面代码会报错，需要提取一个通道分析，x[:,0]或x[:,1]。如过不想自己录，可以用一些网站现成的，比如下面这个https://spandh.dcs.shef.ac.uk//gridcorpus/ 获取音频基本信息：

import numpy as np
from scipy.io import wavfile
import matplotlib.pylab as plt

#%%
Fs,x =wavfile.read('swio1s.wav') #读取音频文件
print('Fs =',Fs,'Hz')
print('Shape:', x.shape)
print('Datatype:', x.dtype)
print('Duration:', round(x.shape[0] / float(Fs), 3),'seconds')
#%%
print(len(x))

画出时域波形

t = np.arange(0,len(x),1)
#%%

plt.plot(t/Fs, x/np.power(2,15))
#plt.plot(t/Fs,x/np.power(2,15))
plt.xlabel('t/s')
plt.grid()
plt.show()

• 先找一段清音

截取其中一段波形，从第12400采样点开始，长度为256

N = 256
Tstart = 12400
s=x/np.power(2,15)
stemp = s[np.arange(Tstart,Tstart+N,1)]
plt.plot(stemp)
plt.show()

傅里叶变换

Stemp = np.fft.fft(stemp) #傅里叶变换
plt.plot(20*np.log10(np.abs(Stemp)))
plt.show()

加窗汉宁窗之后

#%% 加汉宁窗
n = np.arange(0,N,1)
w = 0.5*(1-np.cos(2*np.pi*n/(N-1)))
#plt.plot(w)
#plt.show()
#%% 加窗之后的频谱
Stemp = np.fft.fft(stemp*w) #傅里叶变换
n = np.arange(0,256,1)
f = n*Fs/1000/N
plt.plot(f, 20*np.log10(np.abs(Stemp[n])))
plt.grid()
plt.xlabel('Frequency/kHz')
plt.ylabel('Amplitude/dB')
plt.show()

• 再找一段浊音，浊音的波形具有周期性

起始采样点

N = 1000
Tstart = 30500
s=x/np.power(2,15)
stemp = s[np.arange(Tstart,Tstart+N,1)]
plt.plot(stemp)
plt.show()

傅里叶变换

Stemp = np.fft.fft(stemp) #傅里叶变换
plt.plot(20*np.log10(np.abs(Stemp)))
plt.show()

加窗之后

#%% 加汉宁窗
n = np.arange(0,N,1)
w = 0.5*(1-np.cos(2*np.pi*n/(N-1)))
#plt.plot(w)
#plt.show()
#%% 加窗之后的频谱
Stemp = np.fft.fft(stemp*w) #傅里叶变换
n = np.arange(0,256,1)
f = n*Fs/1000/N
plt.plot(f, 20*np.log10(np.abs(Stemp[n])))
plt.grid()
plt.xlabel('Frequency/kHz')
plt.ylabel('Amplitude/dB')
plt.show()

简要分析：

浊音:时域上为短时周期信号，波形呈现周期性重复

频域上信号能量大，衰减快，共振峰明显

清音:时域上波形不具有周期性

频域上能量小，衰减相对慢。

窗函数主瓣窄，旁瓣大，频率识别精度高，幅值识别精度低

完整代码

import numpy as np
from scipy.io import wavfile
import matplotlib.pylab as plt

#%%
Fs,x =wavfile.read('swio1s.wav') #读取音频文件
print('Fs =',Fs,'Hz')
print('Shape:', x.shape)
print('Datatype:', x.dtype)
print('Duration:', round(x.shape[0] / float(Fs), 3),'seconds')
#%%
print(len(x))
#%%
t = np.arange(0,len(x),1)
#%%

plt.plot(t/Fs, x/np.power(2,15))
#plt.plot(t/Fs,x/np.power(2,15))
plt.xlabel('t/s')
plt.grid()
plt.show()

#%%
#
N = 256
Tstart = 12400
s=x/np.power(2,15)
stemp = s[np.arange(Tstart,Tstart+N,1)]
plt.plot(stemp)
plt.show()

Stemp = np.fft.fft(stemp) #傅里叶变换
plt.plot(20*np.log10(np.abs(Stemp)))
plt.show()

#%% 加汉宁窗
n = np.arange(0,N,1)
w = 0.5*(1-np.cos(2*np.pi*n/(N-1)))
#plt.plot(w)
#plt.show()
#%% 加窗之后的频谱
Stemp = np.fft.fft(stemp*w) #傅里叶变换
n = np.arange(0,256,1)
f = n*Fs/1000/N
plt.plot(f, 20*np.log10(np.abs(Stemp[n])))
plt.grid()
plt.xlabel('Frequency/kHz')
plt.ylabel('Amplitude/dB')
plt.show()



N = 1000
Tstart = 30500
s=x/np.power(2,15)
stemp = s[np.arange(Tstart,Tstart+N,1)]
plt.plot(stemp)
plt.show()

Stemp = np.fft.fft(stemp) #傅里叶变换
plt.plot(20*np.log10(np.abs(Stemp)))
plt.show()

#%% 加汉宁窗
n = np.arange(0,N,1)
w = 0.5*(1-np.cos(2*np.pi*n/(N-1)))
#plt.plot(w)
#plt.show()
#%% 加窗之后的频谱
Stemp = np.fft.fft(stemp*w) #傅里叶变换
n = np.arange(0,256,1)
f = n*Fs/1000/N
plt.plot(f, 20*np.log10(np.abs(Stemp[n])))
plt.grid()
plt.xlabel('Frequency/kHz')
plt.ylabel('Amplitude/dB')
plt.show()