Python 以正弦信号绘制信号

Python 以正弦信号绘制信号

学习目标:
熟悉导入库以扩展 Python 功能的概念
能够生成、可视化（和回放）从 Internet 下载的正弦信号和音频 (WAVE) 文件
对离散时间索引 k 或连续时间 t 上的绘图信号有初步了解
了解 Python 库中函数的基本用法
Matplotlib for graphical output
(like plot()), xlabel(), ylabel(), legend(), subplot(), etc.)
Numpy for basic mathematical functions
(like arange(), linspace(), sin(), log10(), etc.)
interactive Python class IPython.display
(in particular the playback function Audio().)

使用外部库
对于语音处理模块，预期的学习成果之一是了解您通常会使用外部库（例如 librosa、PyAudio 或类似库）的信号处理算法的基本功能。 在了解底层原理之前，我们会尽量避免使用这些库。
但是，以下库（因此我们将经常使用）提供了一些非常基本的功能：

#Let's do some necessary and nice-to-have imports
%matplotlib inline
import matplotlib.pyplot as plt  # plotting
import numpy as np               # math

其他我们不太需要的库将仅在某些笔记本中需要时导入（例如，下面使用的库声音文件允许从/向硬盘加载和保存音频文件）。 如果您在使用此类库时遇到错误，您可以使用命令 !pip install 安装它们。
安装后，您可以使用注释字符 # 再次注释该行，以避免每次执行单元格时重复安装过程。另请注意，如果您安装了库，则可能需要重新启动 iPython 内核。

任务 1：生成、显示和收听正弦信号
正弦信号:x(t)=sin(2πft+ϕ)
可以使用 numpy 的正弦函数 np.sin() 生成正弦信号。
为了实现 \eqref{eq:sin}，我们首先需要一个时间向量 t，它可以使用 numpy 的 arange() 函数创建。

任务1：
以 fs=8000 Hz 的采样频率生成长度为 14 秒且频率为 f=250 Hz 的正弦信号。 显示时域信号并回放。 请注意，正弦信号的播放音量可能很大，因此请注意在播放过程中将音量调整到较低且舒适的水平。
首先生成时间向量 t 是有意义的。 考虑使用命令 np.arange(…) 或 np.linspace(…) 之一（您可以通过访问通用 numpy 文档的链接来熟悉这两个函数）。
播放不同的频率。 通过设置 φ≠0 来玩转相移。

# Let's do some 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

# Time specifications:
fs = 8000                        # samples per second
dt = 1/fs                        # seconds per sample
StopTime = 0.25                  # length in seconds
t = np.arange(0,StopTime,dt)     # time vector in seconds

# sine signal of frequency 250 Hz
f = 250                          # frequency in hertz
x = np.sin(2*np.pi*f*t)          # sine wave signal

plt.plot(t,x)
plt.ylabel('$x(t)$')
plt.xlabel('time $t$ in seconds');

• $(pwd)：表示当前所在目录。
  ${ }：获取变量的结果。
  一般情况下，$var与${var}是没有区别的，但是用${ }会比较精确的界定变量名称的范围。
  
  $0：shell 命令本身。
  
  $1 到 $9： 表示 Shell 的第几个参数。
  

IPython.display.Audio() 创建一个可用于播放信号的播放器对象。 由于我们想将变量 x 作为输入参数，我们还必须提供采样频率 fs（变量速率 =）。

import IPython.display as ipd
ipd.Audio(x/10,rate=fs)           # dividing the level by 10 ensures that playback is not too loud

可以添加两个（或更多）正弦信号（我们称之为假设）。
x(t)=sin(2πf1t+φ1)+sin(2πf2t+φ2)
频率 f1=60 Hz 和 f2=200 Hz。
我们将信号的长度增加一点到 0.5sec 。

StopTime = 0.5                   # length in seconds
t = np.arange(0,StopTime,dt)     # time vector in seconds

# sine signal of frequency 60 Hz
f1 = 60                          # frequency in hertz
phi1 = np.pi/2                   # phase (between 0 and 2 pi)
x1 = np.sin(2*np.pi*f1*t+phi1)   # sine wave signal

# sine signal of frequency 100 Hz
f2 = 200                         # frequency in hertz
phi2 = 0                         # phase (between 0 and 2 pi)
x2 = np.sin(2*np.pi*f2*t+phi2)   # sine wave signal

# Plot the signal over time:
fig = plt.figure(figsize=(12,5)) # create a figure of size 12x5 inch
plt.subplot(2,1,1)               # 1st subplot in an array of 2x1 subplots
plt.plot(t,x1,lw=2,label='$x_1(t)$=sin($2 \cdot \pi \cdot ' + str(f1) + '$ Hz $\cdot t)+ $' + str(phi1))
plt.ylabel('$x_1(t)$')           # we can use LaTeX notation, e.g. LaTeX math mode (using $...$)
plt.legend()                     # display the label information as a legend
plt.subplot(2,1,2)               # 1st subplot in an array of 2x1 subplots
plt.plot(t,x2,'r',lw=2, label='$x_2(t)$=sin($2 \cdot \pi \cdot ' + str(f2) + '$ Hz $\cdot t) + $' + str(phi2))
plt.ylabel('$x_2(t)$')
plt.legend()
plt.xlabel('time $t$ in seconds');

Matplotlib.pyplot.legend()
图例是描述图形元素的区域。 在 matplotlib 库中，有一个名为 legend() 的函数，用于在轴上放置图例。
legend() 中的属性 Loc 用于指定图例的位置。loc 的默认值是 loc=”best”（左上角）。 字符串“左上”、“右上”、“左下”、“右下”将图例放置在轴/图形的相应角上。
legend()函数的属性bbox_to_anchor=(x, y)用于指定图例的坐标，属性ncol表示图例的列数，默认值为1。

• why I comment to legend() it still print same graphs?
• what does {t,x1,lw=2,label=…} mean?

任务 2：正弦信号的叠加
可以添加两个（或更多）正弦信号（我们称之为假设）。
x(t)=sin(2πf1t+φ1)+sin(2πf2t+φ2)
生成两个不同频率的正弦函数 f1 和 f2:
使用 subplot 命令在一个图中显示两个信号以及两个信号的和。
播放不同的频率。 为正弦参数添加一个相位，并使用不同的相位 ϕ1 和 ϕ2 。 考虑收听和信号。

# add the two signals already created before
x3=x1+x2

# plot
fig = plt.figure(figsize=(12,8))
plt.subplot(3,1,1)
plt.plot(t,x1,lw=2, label='$x_1(t)$=sin($2 \cdot \pi \cdot ' + str(f1) + '$ Hz $\cdot t)+ $' + str(phi1))
plt.ylabel('$x_1(t)$')
plt.xlabel('time $t$ in seconds');
plt.legend()
plt.subplot(3,1,2)
plt.plot(t,x2,color='r',lw=2, label='$x_1(t)$=sin($2 \cdot \pi \cdot ' + str(f2) + '$ Hz $\cdot t) + $' + str(phi2))
plt.ylabel('$x_2(t)$')
plt.xlabel('time $t$ in seconds');
plt.legend()
plt.subplot(3,1,3)
plt.plot(t,x3,lw=2,color='g',label='sin1+sin2')
plt.ylabel('$x_1(t)+x_2(t)$')
plt.xlabel('time $t$ in seconds');

任务 3：处理音频文件
为了处理语音信号，我们要加载一个声音文件。 以下命令将文件 Speech.wav 从 Internet 下载到硬盘。
函数 sf.read(…) 将 WAVE 文件的内容加载到内存中。 变量 s 包含语音信号向量 s[k] ，变量 fs 包含采样率，也就是采样频率 fs ，即每秒信号的离散样本数。

# download wave file
!curl https://staffwww.dcs.shef.ac.uk/people/S.Goetze/sound/speech.wav -o speech.wav 

!pip install soundfile

import soundfile as sf
s_file_name = 'speech.wav'
s,fs = sf.read(s_file_name)

print(s.shape) # this prints the dimension (length) of the vector s 
               # (not necessary in the following but to get an idea about the size (shape))

plt.figure(figsize=(6,4)) # create empty figure of size 6x4 inch
plt.plot(s)
plt.ylabel('$s[k]$');
plt.xlabel('sample $k$');

# optional we can create a matching time vector 
plt.figure(figsize=(6,4))
t=np.linspace(0,len(s)/fs,len(s))
plt.plot(t,s)
plt.xlabel('$t$ in seconds');
plt.ylabel('$s(t)$');

上面显示了两种不同的可视化，它们仅在时间轴方面有所不同。 上图显示 s[k] ，即离散时间索引 k 上离散时域中的信号变量。 下图显示了相同的信号，但在连续信号 s(t) 中。 请注意：当然变量 s 在计算机中是离散的，即我们实际上显示 s[k] 。 但是由于采样点足够接近，所以在选择适当的时间轴和轴描述的情况下，它看起来好像是一个时间连续信号 s(t)。

任务 4：分贝 (dB) - 以线性幅度和对数幅度绘制信号
术语分贝（通常缩写为 dB）通常用于讨论或可视化与声学系统和信号（以及语音和音频之外的各种不同领域）相关的电平、振幅、功率等。
信号（例如 s[k] ）可以“正常”显示（即具有线性幅度）或以 dB 为单位。 由于人类以 dB 为单位感知声音，因此熟悉如何以 dB 为单位显示信号是有意义的。
要绘制幅度信号 s[k]（不是信号 s2[k] 的幂），我们可以使用 (\ref{eq:defDBamplitude})。

sdB[k]=20⋅log10(|s[k]+ε|)

变量 ϵ 在上面被介绍为“一个非常小的数字”，以避免未定义的 log(0)。 由于对于实数值，对数仅定义为正数值，因此绝对值 |⋅| 在计算对数之前应用。

以分贝 (dB) 比例绘制上面加载的语音信号。

plt.subplot(2,1,1)
plt.plot(s)
plt.ylabel('$s[k]$');

epsilon = 1e-10 # very small number to avoid log10(0)
s_dB = 20*np.log10(np.abs(s)+epsilon)     

plt.subplot(2,1,2)
plt.plot(s_dB)
plt.ylabel('$s_{\mathrm{dB}}[k]$')
plt.xlabel('sample $k$')
plt.ylim((-40,0));                    # define which part of the y-axis should be 'zoomed in'

e 是科学计数法的一种表示
eN: 10 的 N 次方
1e2=1 * 10^2=100
1.2e-5=1.2 * 10^(-5)=0.000012