python 频谱图的应用和递归过滤器

预期的学习目标（ILOS）：
您应该
演示频谱图的应用以及各自的numpy命令的参数如何影响时间和频率分辨率，并通过thet，您可以分析的内容。
能够使用Butterworth，Cheby I和II以及Cauer / Cauer / Elliptic Filter Designs设计递归过滤器。
线性绘制滤波器传输函数，并在DB中过滤信号并通过分析频谱图来感知过滤前后过滤之间的差异
了解Python库的基本用途

# Let's do the ususal necessary and nice-to-have imports
%matplotlib inline
import matplotlib.pyplot as plt  # plotting
import numpy as np               # math!pip install seaborn 
import seaborn as sns; sns.set() # styling (optional)# download speech and noise example files
s_file_name = 'speech_8kHz_murder.wav'
!curl https://staffwww.dcs.shef.ac.uk/people/S.Goetze/sound/{s_file_name} -o {s_file_name} import soundfile as sf
from IPython import display as ipd# load speech wave into variable
s, fs = sf.read(s_file_name)print('File "' + s_file_name + '" loaded. Its sampling rate is ' + str(fs) + ' Hz.')# listen to the sound file (if you want)
ipd.Audio(s, rate=fs)

频谱图（任务1）

频谱图可用于同时可视化时间和频率信息。 Matplotlib的Specgram（）命令具有可选参数，例如 对于FFT长度。 看看文档，以使自己与参数相关。

# visualise
fig=plt.figure(figsize=(12,6)) # create a figure of size 12 x 6 inchesLDFT=512plt.subplot(1,2,1)
plt.specgram(s, Fs=fs, NFFT=LDFT); # NFFT default: 256
plt.title('spectrogram with DFT length of ' + str(LDFT) + ' (high time res.)')
plt.xlabel('time $t$')
plt.ylabel('frequency $f$')
plt.grid(False) # no grid (in case you used seaborn)plt.colorbar(label='dB');
plt.clim(-150,0)LDFT=8192plt.subplot(1,2,2)
plt.specgram(s, Fs=fs, NFFT=LDFT); # NFFT default: 256
plt.title('spectrogram with DFT length of ' + str(LDFT) + ' (high freq. res.)')
plt.xlabel('time $t$')
plt.ylabel('frequency $f$')
plt.grid(False)  # no grid (in case you used seaborn)plt.colorbar(label='dB');
plt.clim(-150,0)plt.tight_layout()

上面的频谱图可视化相同的语音信号，但是对于左侧，选择了较小的DFT长度（在这种情况下为FFT长度），即每个分析块的样品较少。 与右图中的频谱图相比，这导致时间方向的分辨率更高。 可以在左图中的频谱图中更好地分析时间结构，而光谱信息（例如共振（右图中的光水平线））在右光谱图中更为突出。

无限脉冲响应（IIR）
是适用于许多线性时间不变的系统的属性，通过具有脉冲响应 h（t）而不同，这些系统的区别是，它并没有完全零以超过某个点的零，但仍在继续无限期。这与有限的脉冲响应（FIR）系统形成鲜明对比，在该系统中，脉冲响应有时{\ displayStyleT t> t} t> t在某些有限的{\ displayStyleT t} t时确实完全为零，因此有限持续时间。线性时间不变系统的常见示例是大多数电子和数字过滤器。具有此属性的系统被称为IIR系统或IIR过滤器。

在实践中，脉冲响应，即使是IIR系统，通常都接近零，并且可以被忽略到特定点。然而，产生IIR或FIR反应的物理系统是不同的，因此区分的重要性。例如，由电阻器，电容器和/或电感器（也许是线性放大器）组成的模拟电子过滤器通常是IIR过滤器。另一方面，基于使用无反馈的延迟线的离散时间过滤器（通常是数字过滤器）必然是FIR过滤器。模拟滤波器中的电容器（或电感器）具有“记忆”，并且它们的内部状态永远不会在冲动后完全放松（假设忽略量子效应的电容器和电感器的经典模型）。但是在后一种情况下，在脉冲达到了点击延迟线的末端之后，系统没有进一步的记忆力，并恢复了其初始状态。它的冲动响应超出了这一点，完全为零。

Butterworth滤波器
是一种信号处理过滤器，旨在具有在PassBand中尽可能平坦的频率响应。 它也称为最大平坦的幅度滤波器。 这是由英国工程师兼物理学家史蒂芬·巴特沃思（Stephen Butterworth）在1930年首次描述的，题为“关于过滤器放大器的理论”。

Chebyshev滤波器
是模拟或数字过滤器，其滚动量比Butterworth过滤器更陡，并且具有PassBand Ripple（I型）或挡侧纹波（II型）。 Chebyshev过滤器具有最小化过滤器范围内理想化和实际过滤器特性之间的误差的属性. 这种类型的过滤器以Pafnuty Chebyshev的名字命名，因为其数学特征源自Chebyshev多项式。 I型Chebyshev过滤器通常称为“ Chebyshev过滤器”，而II型过滤器通常称为“逆Chebyshev过滤器”。

由于Chebyshev滤波器固有的PassBand涟漪，在PassBand中具有更平滑的响应过滤器，但对于某些应用程序，首选停止带中的响应更加不规则的响应

椭圆滤波器
（以Wilhelm Cauer或Zolotarev滤波器的名字命名，在Yegor Zolotarev之后）是一个信号处理滤波器，在Passband和stopband中都具有均衡的纹波（Equiripple）行为。 每个频段中的连锁量都是独立调节的，对于给定的纹波值，没有其他相等顺序的过滤器可以在通带和停止带之间的增益更快地过渡（是否均衡）。 或者]或者，人们可以放弃能够独立调整通过带和停止带纹波的能力，而是设计一个对组件变化最大程度不敏感的过滤器。

当停止带中的涟漪接近零时，过滤器变为I型Chebyshev滤波器。 随着通带中的涟漪接近零，过滤器变为II型Chebyshev滤波器，最后，随着两个连锁值接近零，滤波器变为Butterworth滤波器。

公差方案
我们想设计一个零件线性滤波器，例如低通，带通，高通或类似的线性滤波器，我们希望这样的过滤器至少满足以下条件：
它应该在传球区域1或接近一个
在停止频段区域中应为0或接近零
以下函数plot_tolerance_scheme（）是一个辅助函数，可视化这种公差方案:

def plot_tolerance_scheme(Wp=0.25,Ws=0.3,Rp_lin=0.9,Rs_lin=0.1):"""Plots a tolerance scheme for a (low-pass) filter design.Parameters----------Wp : float, optional (but recommended)pass limit frequency (normalised) $W_p = \Omega / \pi$Ws : float, optional (but recommended)stop limit frequency (normalised)Rp_lin : float, optional (but recommended)allowed ripple range in pass band (linear), default 0.9Rs_lin : float, optional (but recommended)allowed ripple range in pass band (linear), default 0.1Example use:-------signalplot_tolerance_scheme(Wp=0.25,Ws=0.3,Rp_lin=0.9,Rs_lin=0.1)"""dh1x=[0,Ws];  dh1y=[1,1];            # (x,y) coordinates of linesdh2x=[0,Wp];  dh2y=[Rp_lin,Rp_lin]; dv2x=[Wp,Wp]; dv2y=[0,Rp_lin];   sh1x=[Ws,1];  sh1y=[Rs_lin,Rs_lin]; sh2x=[Wp,1];  sh2y=[0,0]; svx=[Ws,Ws];  svy=[Rs_lin,1];  # plot the actual linesplt.plot(dh1x,dh1y,'k--',dh2x,dh2y,'k--',dv2x,dv2y,'k--',sh1x,sh1y,'k--',sh2x,sh2y,'k--',svx,svy,'k--');plt.xlabel('Frequency $\Omega/\pi$');

在下面的公差方案中，对于Exmample值WP = 0.25，WS = 0.3，以下条件可视化。 允许的公差是通过参数RP_LIN定义的传递区域，允许的公差是RS_LIN的停止频段区域。 请注意，它们具有结尾… lin，因为我们以后需要这些数量来对数缩放（那时我们会关心）。

Wp=0.25;    # passband edge frequency 
Ws=0.3;     # stopband edge frequency 
Rp_lin=0.9; # allowed ripples in the pass band area
Rs_lin=0.1; # allowed ripples in the stop band areaplot_tolerance_scheme(Wp,Ws,Rp_lin,Rs_lin)

请注意，低通差方案是根据归一化频率Ω 给出的。
提醒我们，我们已经熟悉了频率的三种不同符号
Ω = w/f = 2pi*f

零极图（z域中过滤器的传输函数）
对于语音和音频处理，系统通常具有实值系数BM∈R和An∈R，从而导致了实现系统。 对于实值系统，电线杆和零是实现的，要么在复杂的共轭对中发生。

对于因果和稳定系统H [Z]的杆子，所有极点都必须位于单位圆内，即| z | = 1。

以下代码实现了另一个辅助功能Zplane（），该ZPLANE（）在复杂的z平面中绘制了极点和零。

def zplane(z, p, title='Poles and Zeros'):"Plots zeros and poles in the complex z-plane"ax = plt.gca()ax.plot(np.real(z), np.imag(z), 'bo', fillstyle='none', ms=10)ax.plot(np.real(p), np.imag(p), 'rx', fillstyle='none', ms=10)unit_circle = plt.Circle((0, 0), radius=1, fill=False,color='black', ls='--', alpha=0.9)ax.add_patch(unit_circle)plt.title(title)plt.xlabel('Re{$z$}')plt.ylabel('Im{$z$}')plt.axis('equal')

下面的代码示例性显示了位置的两个（共轭复合物）零，Z0,1 = -0.9和Z0,2 = 0.8（由蓝色圆圈表示）以及位置Z∞的杆子，1 = 0.5+0.5J和Z ∞，2 = 0.5-0.5J（由红色十字表示）z平面。 上面的代码是绘制零和极线的辅助功能。

poles = [0.5-0.5j, 0.5+0.5j] # two (comjugate complex) poles
zeros = [-0.9, 0.8]          # two (real-valued) zeros
zplane(zeros,poles)

Butterworth过滤器设计（任务3）
要实际设计一个过滤器，我们可以使用Scipy的信号库中的function Butter（），该信号库可以计算过滤器系数A和BM。

为了确保我们必须知道我们的tolerace方案，我们必须知道，对于0≤m≤q，对于0≤n≤p和bm，多少个系数A（定义）butter（ ）方法需要实现这一目标。如果我们假设具有与反馈路径系数AN一样多的正向系数BM，即p = q，我们将滤波器的顺序定义为P（或Q）。函数buttord（）采用我们的公差方案（值），并计算为所需的过滤顺序n。此外，函数buttord（）返回过滤器的截止频率（又称固有频率），因此在下面使用了可变名称wn（如果我们用FC表示此频率，这将导致与中心频率FC的混淆）。截止频率定义为3db衰减频率响应的频率

from scipy import signal as sig         # filter designs (if not already imported)Wp=0.25;    # passband edge frequency 
Ws=0.3;     # stopband edge frequency 
Rp_lin=0.9; # allowed ripples in the pass band area
Rs_lin=0.1; # allowed ripples in the stop band area# the ripple hights in dB are needed in later cells of this Jupyter Notebook 
Rp=-20*np.log10(Rp_lin);  # max passband ripple in dB
Rs=-20*np.log10(Rs_lin);  # min stopband attenuation in dB# get lowest filter order N to fullfill requirements above
N, Wn = sig.buttord(Wp, Ws, Rp, Rs)print('The minimum possible filter order to fulfil the tolerance scheme is '+str(N)+'.')
print('The cut-off frequency which will be {:.2f}.'.format(Wn)) # format number - two digits after decimal pt.# design low-pass filter of order N using butterworth method
b, a = sig.butter(N, Wn, 'low')print(b)
print(a)

我们可以使用Scipy的信号处理库中的FREQZ（）函数，该库计算了频域中的传输函数H（f）。 我们还可以通过利用freqz（）来计算频率响应

f,h=sig.freqz(b,a)
omega=np.linspace(0,1,len(f))# the following code does actually the same as the freqz command:
#h=np.abs(np.fft.fft(b,1024))/np.abs(np.fft.fft(a,1024));
#h=h[0:513] # only show first half (positive frequencies)
#omega=np.linspace(0,1,513)# plot frequency response
plot_tolerance_scheme(Wp,Ws,Rp_lin,Rs_lin)
plt.plot([Wn,Wn],[0,1],color='red',ls=':',label='cutoff frequency')
plt.plot(omega,np.abs(h), lw=2, label='Butterworth low-pass')
plt.title('Butterworth low-pass filter of order ' + str(N))
plt.ylabel('Amplitude $|h(e^{j \Omega})|$')
plt.legend();

要查看相应的极点和零，我们可以使用Numpy的Roots（）函数将系数转换为零和极点。

请注意：我们不计算系数BM的根部以获取以下图的零位置。 Butterworth滤波器理论表明，所有零均处于z = -1，并且Roots（）函数将引入一些计算误差，显示零不完全是z = -1

# plot zeros and poles in the z plane
zplane(-1, np.roots(a))
plt.text(-0.95,0.1,'('+str(N)+')'); # adding text with info how many zeros are located a z=-1

任务3（更改斜率 /过渡频率）被求解如下：

wp=0.475;    # passband edge frequency 
ws=0.525;     # stopband edge frequency 
Rp_lin=0.95; # allowed ripples in the pass band area
Rs_lin=0.05; # allowed ripples in the stop band area# the ripple hights in dB are needed in later cells of this Jupyter Notebook 
Rp=-20*np.log10(Rp_lin);  # max passband ripple in dB
Rs=-20*np.log10(Rs_lin);  # min stopband attenuation in dB# get lowest filter order N to fullfill requirements above
N, Wn = sig.buttord(wp, ws, Rp, Rs)
# design low-pass filter of order N using butterworth method
b, a = sig.butter(N, Wn, 'low')h=np.abs(np.fft.fft(b,1024))/np.abs(np.fft.fft(a,1024));
h=h[0:513] # only show first half (positive frequencies)
omega=np.linspace(0,1,513)# plot frequency response
plot_tolerance_scheme(wp,ws,Rp_lin,Rs_lin)
plt.plot(omega, abs(h), lw=2)
plt.title('Butterworth low-pass filter of order ' + str(N))
plt.ylabel('Amplitude $|h(e^{j \Omega})|$');

过滤（任务4）
Scipy的信号库中的功能Filtfilt（B，A，S）实现了一个过滤器。
我们可以在过滤前后比较信号s [k]。 以下代码生成的图显示了在应用过滤器（左图）和相同信号后（右图；仅在X轴（时间 - ）方向上放大的信号。

s_filtered = sig.filtfilt(b, a, s)plt.figure(figsize=(12,4))plt.subplot(1,2,1) # plot signal before and after filtering in time domain
plt.plot(s,label='$s[k]$ (before filtering)')
plt.plot(s_filtered,label='after filtering')
plt.legend()plt.subplot(1,2,2) # zoom into the signal 
start_idx=16000
end_idx=start_idx+1500
plt.plot(np.arange(start_idx,end_idx),s[start_idx:end_idx],label='$s[k]$ (before filtering)')
plt.plot(np.arange(start_idx,end_idx),s_filtered[start_idx:end_idx],label='after filtering')
plt.title('Cut-out showing a harmonic (voiced) piece and a noisy (unvoiced) one')
plt.legend();

左面板显示，我们可以看到信号的微小变化，并在过滤前后进行比较。 但是，信号似乎没有太多发生。 如果我们对16,000和17,500（右图）的样本之间的零件进行更详细的看法，该零件在此切割开始时显示了周期性的（声音）零件），则是一个嘈杂的部分（相当无用）和（在 结束）另一个过渡到一个不太嘈杂的部分，我们观察到右面板中间的嘈杂部分的更多变化。

聆听原始和过滤的信号如下所示，表明低通效应显然是可感知的。

# listen to the sound file (if you want)
ipd.Audio(s, rate=fs)
# listen to the sound file (if you want)
ipd.Audio(s_filtered, rate=fs)

由于语音在低频范围内具有大部分能量，因此大多数能量仍然是保存的（因此上面的时间域相似）。 如果我们计算镜头（下面的代码），我们清楚地看到过滤后缺少高频（下图右图）。 高频率是指时域信号中的快速振荡，这应该解释为什么上述时域图中的嘈杂部分会受到更大的影响。

# visualise as spectrograms
fig=plt.figure(figsize=(8,4)) # create a figure of size 12 x 6 inchesplt.subplot(1,2,1)
plt.specgram(s, Fs=fs)
plt.title('spectrogram of signal before filtering')
plt.xlabel('time $t$')
plt.ylabel('frequency $f$')
plt.grid(False) # no grid (in case you used seaborn)plt.colorbar(label='dB');
plt.clim(-180,-30) # select the 'color'/amplitude range (in dB)plt.subplot(1,2,2)
plt.specgram(s_filtered, Fs=fs)
plt.title('spectrogram of filtered signal')
plt.xlabel('time $t$')
plt.ylabel('frequency $f$')
plt.grid(False)  # no grid (in case you used seaborn)plt.colorbar(label='dB');
plt.clim(-180,-30)plt.tight_layout() # to see all axis descriptions

Chebyshev I和Chebyshev II过滤设计（任务5）
Chebyshev滤波器设计分别允许在传球带（Chebyshev I Filter Design）或Stop Band（Chebyshev II过滤器设计）中进行涟漪，从而更好地利用了公差方案。 与Butterworth滤波器相比，可以减少必要的过滤系数N的数量。

wp=0.25;    # passband edge frequency 
ws=0.3;     # stopband edge frequency 
Rp_lin=0.9; # allowed ripples in the pass band area
Rs_lin=0.1; # allowed ripples in the stop band area# get lowest filter order N to fullfill requirements above
N, Wn = sig.cheb1ord(wp, ws, Rp, Rs)
# design low-pass filter of order N using Chebyshef I method
b, a = sig.cheby1(N, Rp, Wn, 'low')h=np.abs(np.fft.fft(b,1024))/np.abs(np.fft.fft(a,1024));
h=h[0:513] # only show first half (positive frequencies)
omega=np.linspace(0,1,513)plt.figure(figsize=(12, 12))# plot frequency response
plt.subplot(2,2,1)
plot_tolerance_scheme()
plt.plot(omega, abs(h), lw=2)
plt.title('Chebyshev type I low-pass filter of order ' + str(N))
plt.ylabel('Amplitude $|h(e^{j \Omega})|$')
#plt.show()# plot zeros and poles in the z plane
plt.subplot(2,2,2)
zplane(np.roots(b), np.roots(a))# get lowest filter order N to fullfill requirements above
N, Wn = sig.cheb2ord(wp, ws, Rp, Rs)
# design low-pass filter of order N using Chebyshef II method
b, a = sig.cheby2(N, Rs, Wn, 'low')h=np.abs(np.fft.fft(b,1024))/np.abs(np.fft.fft(a,1024));
h=h[0:513] # only show first half (positive frequencies)
omega=np.linspace(0,1,513)# plot frequency response
plt.subplot(2,2,3)
plot_tolerance_scheme()
plt.plot(omega, abs(h), lw=2)
plt.title('Chebyshev type II low-pass filter of order ' + str(N))
plt.ylabel('Amplitude $|h(e^{j \Omega})|$')
#plt.show()# plot zeros and poles in the z plane
plt.subplot(2,2,4)
zplane(np.roots(b), np.roots(a))

我们可以再次应用过滤器，并在过滤前后分析频谱图。

s_filtered = sig.filtfilt(b, a, s)# visualise as spectrograms
fig=plt.figure(figsize=(8,4)) # create a figure of size 12 x 6 inchesplt.subplot(1,2,1)
plt.specgram(s, Fs=fs)
plt.title('spectrogram of signal before filtering')
plt.xlabel('time $t$')
plt.ylabel('frequency $f$')
plt.grid(False) # no grid (in case you used seaborn)plt.colorbar(label='dB');
plt.clim(-180,-30)plt.subplot(1,2,2)
plt.specgram(s_filtered, Fs=fs)
plt.title('spectrogram of filtered signal')
plt.xlabel('time $t$')
plt.ylabel('frequency $f$')
plt.grid(False)  # no grid (in case you used seaborn)plt.colorbar(label='dB');
plt.clim(-180,-30)plt.tight_layout() # to see all axis descriptions

在上方的频谱图中，由单位圆上的三个零引起的Chebyshev II滤波器设计的涟漪设计在右图中清晰可见。

椭圆滤波器设计（任务6）
与Butterworth相比，可以进一步减少必要的滤波器系数n，而Chebyshev I和II过滤器可以通过允许波带和停止带中的涟漪来进行。 这是由CAUER过滤器完成的。

# get lowest filter order N to fullfill requirements above
N, Wn = sig.ellipord(wp, ws, Rp, Rs)
# design low-pass filter of order N using butterworth method
b, a = sig.ellip(N, Rp, Rs, Wn, 'low')h=np.abs(np.fft.fft(b,1024))/np.abs(np.fft.fft(a,1024));
h=h[0:513] # only show first half (positive frequencies)
omega=np.linspace(0,1,513)# plot frequency response
plot_tolerance_scheme()
plt.plot(omega, abs(h), lw=2)
plt.title('Cauer (elliptical) low-pass filter of order ' + str(N))
plt.ylabel('Amplitude $|h(e^{j \Omega})|$')
plt.show()

绘图DB（任务7）
过滤器（尤其是用于处理语音或音乐）通常不仅考虑了对数振幅的设计，而且还以线性振幅绘制。 滤波器设计方法也可以用于设计其他零件线性过滤器，例如 乐队。

让我们首先定义所需的特征，并确定设计过滤器所需的顺序。

Wp = [0.25, 0.5]       # pass-band frequency limits (normalised)
Ws = [0.2,  0.6]       # stop-band frequency limits (normalised)
Rp = 1                 # we allow 1 dB ripple in pass-band
Rs = 40                # we's like to have 40dB attenuation
Rs_lin= 10**(-Rs/20)   # transforming dB back to linear
Rp_lin= 10**(-Rp/20)   # transforming dB back to linear# determine necessary filter order as well as cut-off frequencies
N, Wn = sig.buttord(Wp, Ws, Rp, Rs)
print('The minimum possible filter order to fulfil the tolerance scheme is '+str(N)+'.')
print('The 1st cut-off frequency which will be {:.2f}.'.format(Wn[0])) 
print('The 2nd cut-off frequency which will be {:.2f}.'.format(Wn[1])) # design filter
b, a = sig.butter(N, Wn, 'band')# create frequency response
f,h=sig.freqz(b,a)
omega=np.linspace(0,1,len(f))# plot filter
plt.figure(figsize=(12,4))
plt.subplot(1,3,1)
plt.plot(omega, np.abs(h))
#plt.title('Butterworth bandpass filter fit to constraints')
plt.ylabel('Amplitude $|h(e^{j \Omega})|$')
plt.xlabel('Frequency $\Omega / \pi$')
plt.grid(True,which='both', axis='both')plt.fill([0,     Ws[0], Ws[0], 0],     [Rs_lin, Rs_lin, 1.1, 1.1], '0.7', lw=0) # stop
plt.fill([Wp[0], Wp[0], Wp[1], Wp[1]], [-.1, Rp_lin, Rp_lin, -.1], '0.7', lw=0) # pass
plt.fill([Ws[1], Ws[1], 1, 1],         [1.1, Rs_lin, Rs_lin, 1.1], '0.7', lw=0) # stop
plt.axis([0, 1, -0.1, 1.1]);# plot filter again (this time in dB)
plt.subplot(1,3,2)
plt.plot([Wn[0],Wn[0]],[-100,2],color='r',ls=':',label='cutoff frequency1')
plt.plot([Wn[1],Wn[1]],[-100,2],color='g',ls=':',label='cutoff frequency2')
plt.plot(omega, 20*np.log10(np.abs(h)),label='filter transfer function')
plt.ylabel('Amplitude $|h(e^{j \Omega})|$ in dB')
plt.xlabel('Frequency $\Omega / \pi$')
plt.legend(loc='lower right')
plt.fill([0,     Ws[0],  Ws[0],   0],    [-Rs, -Rs, 2, 2],       '0.7', lw=0) # stop
plt.fill([Wp[0], Wp[0],  Wp[1],  Wp[1]], [-100, -Rp, -Rp, -100], '0.7', lw=0) # pass
plt.fill([Ws[1], Ws[1], 1, 1],           [2, -Rs, -Rs, 2],       '0.7', lw=0) # stop
plt.axis([0, 1, -70, 2]);# plot zeros and poles in the z plane
plt.subplot(1,3,3)
zplane(np.roots(b), np.roots(a))plt.tight_layout()s_filtered = sig.filtfilt(b, a, s)# visualise as spectrograms
fig=plt.figure(figsize=(8,4)) # create a figure of size 12 x 6 inchesplt.subplot(1,2,1)
plt.specgram(s, Fs=fs)
plt.title('spectrogram of signal before filtering')
plt.xlabel('time $t$')
plt.ylabel('frequency $f$')
plt.grid(False) # no grid (in case you used seaborn)plt.colorbar(label='dB');
plt.clim(-180,-30)plt.subplot(1,2,2)
plt.specgram(s_filtered, Fs=fs)
plt.title('spectrogram of filtered signal')
plt.xlabel('time $t$')
plt.ylabel('frequency $f$')
plt.grid(False)  # no grid (in case you used seaborn)plt.colorbar(label='dB');
plt.clim(-180,-30)plt.tight_layout() # to see all axis descriptions