python----根据共振峰频率绘制二阶谐振曲线

# 根据共振峰频率绘制二阶谐振曲线

from matplotlib import pyplot as plt
import numpy as np
from scipy import signal
f = np.array([500,1500,2500])
sampleRate = 8000
pitch = 100  #基音
f1,f2,f3= f[0],f[1],f[2]  # 三个共振频率
#冲激函数
yt = np.zeros((8000))
yt[0]=1

if f1>0:
    cft = f1 / sampleRate
    bw = 50
    q = f1 / bw
    rho = np.exp(-np.pi * cft / q)
    theta = 2 * np.pi * cft * np.sqrt(1-1/(4*q*q))
    a2 = -2 * rho*np.cos(theta)
    a3 = rho * rho
    y = signal.filtfilt([1+a2+a3],[1,a2,a3],yt)

plt.figure()
N = len(y)
fn = np.arange(0,N-1)*sampleRate/N
fftg = np.fft.fft(y)
disg = 20 * np.log10(np.abs(fftg))
plt.plot(fn[1:N//2+1],disg[1:N//2+1])
plt.xlabel('fre/hz')
plt.ylabel('amp/db')
plt.title('first formant')

if f2>0:
    cft = f2 /sampleRate
    bw = 50
    q = f2 / bw
    rho = np.exp(-np.pi * cft / q)
    theta = 2 * np.pi * cft * np.sqrt(1-1/(4*q*q))
    a2 = -2*rho*np.cos(theta)
    a3 = rho * rho
    y = signal.filtfilt([1+a2+a3],[1,a2,a3],y)

plt.figure()
N = len(y)
fn = np.arange(0,N-1)*sampleRate/N
fftg=np.fft.fft(y)
disg=20*np.log10(np.abs(fftg))
plt.plot(fn[1:N//2+1],disg[1:N//2+1])
plt.xlabel('fre/hz')
plt.ylabel('amp/db')
plt.title('second formant')


if f3>0:
    cft=f3/sampleRate
    bw=50
    q=f3/bw
    rho=np.exp(-np.pi*cft/q)
    theta = 2*np.pi*cft*np.sqrt(1-1/(4*q*q))
    a2=-2*rho*np.cos(theta)
    a3=rho*rho
    y = signal.filtfilt([1+a2+a3],[1,a2,a3],y)

plt.figure()
N=len(y)
fn = np.arange(0,N-1)*sampleRate/N
fftg = np.fft.fft(y)
disg=20*np.log10(np.abs(fftg))
plt.plot(fn[1:N//2+1],disg[1:N//2+1])
plt.xlabel('fre/hz')
plt.ylabel('amp/db')
plt.title('third formant')
plt.show()

结果是：

可以看到，共振的地方在500，1500，2500hz的频率上面。 

其中滤波器函数用到了scipy库里面的函数。signal.filtfilt()函数。具体的可以进这个网址了解：

利用Python scipy.signal.filtfilt() 实现信号滤波_John-Cao的博客-CSDN博客_signal.filtfilthttps://blog.csdn.net/weixin_37996604/article/details/82864680