《Python数字信号处理应用》学习笔记——第二章谐波_信号波形图中如何截取波形信号片段 python-优快云博客

本文链接：https://blog.youkuaiyun.com/lstef/article/details/139068840

专栏总目录

本章了解多种波形，观察频谱理解他们的谐波结构，也就是构成他们的三角函数稽核的结构。

介绍数字信号处理中最重要的现象之一：混叠，同时解释Spectrum类的工作原理

2.1 三角波

一个正弦波，仅包含一个频率元素，所以其频谱只有一个尖峰。

复杂波形，例如，下图小提琴录音片段，其离散傅里叶变换的结果中有很多尖峰。

本章继续研究波形和频谱之间的关系。

从三角波开始，三角波就像是正弦波的直线版本。图2-1显示的三角波的频率是200Hz。

图 2-1 200Hz的三角波片段

2.1.1 创建TriangleSignal类

可以使用thinkdsp.TriangleSignal生成三角波：

class TriangleSignal(Sinusoid):
    def evaluate(self, ts):
        # 获得周期数
        cycles = self.freq * ts + self.offset / PI2
        # 获取周期数的小数部分，忽略整数部分
        frac, _ = np.modf(cycles)
        # 减小振幅
        ys = np.abs(frac - 0.5)
        # unbias——将波形中心放置0的位置
        # normalize最大化波形（音量）
        ys = normalize(unbias(ys), self.amp)
        return ys

TriangleSignal从Sinusoid继承了__init__，所以它们使用的参数是一样的，都有freq、amp和offset。

唯一的差别在于evalute。如我们在前面看到过的，ts是我们相对信号求值的采样时间的序列。

产生三角波的方法很多。其细节并不重要，但是我们需要知道evaluate的工作原理：

1. cycles是自起始时间的周期数。np.modf是将周期的数量分解成小数部分和整数部分，前者被存于frac，后者被忽略了。

2. frac序列介于0~1，其频率是给定的。减去0.5之后其值介于-0.5~0.5。对其取绝对值就得到了在0~0.5穿梭的波形。

3.unbias将波形下移，使得其中心位于0，然后normalize将波形归一化到给定的振幅amp。

2.3.2 创建三角波信号

以下是产生图2-1的代码（可以在jupyter notebook中直接运行获得下图）：

from thinkdsp import TriangleSignal
from thinkdsp import decorate

# 创建三角波信号
signal = TriangleSignal(200)
# 截取三个周期
duration = signal.period*3
# 创建wave数据
segment = signal.make_wave(duration, framerate=10000)

#图形化数据
segment.plot()
decorate(xlabel='Time (s)')

如果在idle或.py文件中执行，需要在上述脚本下面追加代码：

#将上面数据保存成图像trianglesignal.png
plt.savefig('trianglesignal.png')

2.3.3 获取三角波声音

如果想要听到该声音：

将上述三角波转换为声音数据

# 获得采样率10000，持续时间0.5秒的波形数据 
wave = signal.make_wave(duration=0.5, framerate=10000)
# 音量最大化
wave.normalize()
# 淡入淡出
wave.apodize()

jupyter notebook上直接执行即可

# 播放
wave.make_audio()

如果在idle或.py文件中执行，需要在转换为声音数据之后追加如下代码：

# 保存音频为WAV文件
from scipy.io.wavfile import write
scaled = np.int16(wave.ys/np.max(np.abs(wave.ys)) * 32767)
write('output.wav', wave.framerate, scaled)
 
# 现在，您可以使用音频播放器打开并播放output.wav文件。

2.3.4 获取三角波频谱

使用上面的signal创建wave，从而使用创建的wave来生成spectrum：

（jupyter notebook中可以直接运行下方代码）

spectrum = wave.make_spectrum()
spectrum.plot()
decorate(xlabel='Frequency (Hz)')

在idle或.py中运行，则需要增加下面代码，运行后通过察看输出文件的方式看到该图像。

#将上面数据保存成图像trianglesignal_spectrum.png
plt.savefig('trianglesignal_spectrum.png')

显示的是其结果的两个视图。下面图中的右图经过了缩放，使得谐波显示得更清楚。最高峰位于基频——200Hz，频率是基频200的倍数的谐波频率峰也被显示在图上。

图2-2 频率200Hz的三角波信号的频谱，以下两种不同的纵坐标尺度显示。右图中切去了基频部分，从而更清楚地显示谐波。

奇怪的是：在偶数倍频率，比如400Hz、800Hz等的位置没有峰值。三角波的谐波全都是基频的奇数倍，比如600Hz、1000Hz、1400Hz等。

上面频谱的另一个特征：谐波的振幅与频率的关系——各谐波的振幅随着频率的平方等比例地衰减。

如：最前面的两个谐波之间的频率（200Hz和600Hz）比率为3，幅度比为9。接下来的两个谐波（600Hz和1000Hz）频率比为1.7，对应的幅度比大约是1.7^2 = 2.9。这种关系被称为谐波结构。

2.2 方波

2.2.1 创建SquareSignal类

thinkdsp还提供了SquareSignal，用来处理方波。其类的定义如下：

class SquareSignal(Sinusoid):
    def evaluate(self, ts):
        cycles = self.freq * ts + self.offset / PI2
        frac, _ = np.modf(cycles)
        ys = self.amp * np.sign(unbias(frac))
        return ys

和SquareSignal一样，它从Sinusoid继承了__init__，所以他们的参数是一样的。同时evaluate地方式也类似。同样，cycles是从起始时间开始的周期数目，而frac是小数部分，其值在每个周期介于0~1。

Unbias将frac移至-0.5~0.5，然后np.sign将负值投射到-1，而正值投射到1。乘以amp之后就得到取值于-amp到+amp之间的方波。

图 2-3 100Hz频率的方波的一部分

图 2-4 100Hz频率方波的频谱

图2-3显示的是频率100Hz的方波的三个周期，图2-4显示的是其频谱。

和三角波一样，方波也只含有奇数谐波，这也就是其频谱中在300Hz、500Hz和700Hz这样的值处有尖峰的原因。但相比来说，它的谐波衰减得要慢一些。也就是，其幅度随频率（而不是频率的平方）等比例地降低。

2.3.2 创建方波信号

生成上述方波的代码如下（jupyter notebook直接运行）：

from thinkdsp import SquareSignal

signal = SquareSignal(200)
duration = signal.period*3
segment = signal.make_wave(duration, framerate=10000)
segment.plot()
decorate(xlabel='Time (s)')

#在idle或.py文件中使用如下方式输出图像
#将上面数据保存成图像squareesignal.png
plt.savefig('squaresignal.png')

2.3.3获取方波声音

将信号转换为声音，聆听

wave = signal.make_wave(duration=0.5, framerate=10000)
wave.apodize()
#jupyter notebook中运行到此步骤即可
wave.make_audio()


# 保存音频为WAV文件
from scipy.io.wavfile import write
scaled = np.int16(wave.ys/np.max(np.abs(wave.ys)) * 32767)
write('output.wav', wave.framerate, scaled)
 
# 现在，您可以使用音频播放器打开并播放output.wav文件。

2.3.4 获取方波频谱

显示频谱图

spectrum = wave.make_spectrum()
spectrum.plot()
decorate(xlabel='Frequency (Hz)')

#在idle或.py文件中使用如下方式输出图像
#将上面数据保存成图像squareesignal.png
plt.savefig('squaresignal.png')

2.3锯齿波信号

2.3.1 创建锯齿波信号

绘制3个周期的片段

from thinkdsp import SawtoothSignal

signal = SawtoothSignal(200)
duration = signal.period*3
segment = signal.make_wave(duration, framerate=10000)
segment.plot()
decorate(xlabel='Time (s)')

#在idle或.py文件中使用如下方式输出图像
#将上面数据保存成图像sawtoothsignal.png
plt.savefig('sawtoothsignal.png')

2.3.2 获取锯齿波声音

wave = signal.make_wave(duration=0.5, framerate=10000)
wave.apodize()
wave.make_audio()


# 保存音频为WAV文件
from scipy.io.wavfile import write
scaled = np.int16(wave.ys/np.max(np.abs(wave.ys)) * 32767)
write('output.wav', wave.framerate, scaled)
 
# 现在，您可以使用音频播放器打开并播放output.wav文件。

2.3.3 获取锯齿波频谱

spectrum = wave.make_spectrum()
spectrum.plot()
decorate(xlabel='Frequency (Hz)')

#在idle或.py文件中使用如下方式输出图像
#将上面数据保存成图像sawtoothsignal_spectrum.png
plt.savefig('sawtoothsignal_spectrum.png')

2.4 混叠

Aliasing在信号处理领域通常指的是“混叠”现象。当一个信号被采样时，如果采样频率不足以捕捉信号中的高频成分，那么这些高频成分在重建过程中可能会错误地显示为低频成分，这种现象就称为混叠。为了防止混叠，采样频率必须至少是信号最高频率成分的两倍，这个原则被称为奈奎斯特采样定理。

图 2-5 一个频率为1100Hz的三角波的频谱，采样率为10000帧每秒。右图显示的是放大之后的谐波。

上图这个波的谐波应该位于3300Hz、5500Hz、7700Hz和9900Hz。上图中和预期的一样，1100Hz和3300Hz处有尖峰，但是第三个尖峰在4500Hz，而不是5500Hz。第4个尖峰在2300Hz而不是7700Hz。如果仔细看的话，应对位于9900Hz的尖峰现在实际上在100Hz处。这是什么问题？

如果对信号在离散的时间点上进行取值，就会失去在采样点之间的信息。对于低频元素，这不是问题，因为会在每个周期中采样多次。

但是如果对频率为5000Hz的信号每秒采样10000次，那在每个周期之中就仅有两个采样。这仅仅是恰好够用，但如果频率更高的画，就不够了。

2.4.1 创建余弦信号

为了演示，我们来产生频率4500Hz和5500Hz的余弦信号，然后对它们每秒采样10000次：

from thinkdsp import CosSignal

# 采样率
framerate = 10000
# 产生4500Hz信号
signal = CosSignal(4500)
# 选取样长度为5个信号周期的时长
duration = signal.period*5
# 生成wave
segment = signal.make_wave(duration, framerate=framerate)
segment.plot()

#在idle或.py文件中使用如下方式输出图像
#将上面数据保存成图像cossignal_4500.png
plt.savefig('cossignal_4500.png')

# 产生5500Hz信号
signal = CosSignal(5500)
# 生成wave
segment = signal.make_wave(duration, framerate=framerate)
segment.plot()

#在idle或.py文件中使用如下方式输出图像
#将上面数据保存成图像cossignal_5500.png
plt.savefig('cossignal_5500.png')

4500Hz余弦信号，波形图如下：

5500Hz余弦信号，波形图如下：

上述脚本，输出结果如图 2-6所示。在绘制信号时，使用的灰色线，采样用的是垂直线表示，以便于比较这两个波形。这样问题清楚了：虽然两个信号是不同的，但是波形是一样的！

图 2-6 频率为4500Hz和5500Hz的余弦信号，每秒采样10000次。信号是不同的，但采样结果是一样的

当我们对一个5500Hz的信号进行10000帧每秒的采样时，其结果与4500Hz的无法区分。同样 7700Hz的信号也无法与2300Hz的信号区分，9900Hz的信号无法与100Hz的区分。

这种效应被称为混叠，也就是采样高频信号后，其结果和采样低频信号的一样。

在这个例子中，我们能测定的最高频率是5000Hz，也就是采样率的一半。高于5000Hz的信号被折叠到5000Hz以下，这也就是这个频率有时被称为 “折叠频率” 的原因。另外一些时候也被称为奈奎斯特频率，参见：http://en.wikipedia.org/wiki/Nyquist_frequency。

如果混叠频率低于0，折叠的模式一样。比如，1100Hz的三角波信号的第5个谐波位于12100Hz。以5000Hz折叠，就会出现在-2100Hz，但是它在0Hz处被再次折叠，从而回到2100Hz。实际上，能在图2-4的2100Hz处观察到一个小峰值，而下一个在4300Hz。

2.4.2 制作一个三角波信号并绘制其频谱。观察谐波如何折叠。

from thinkdsp import TriangleSignal


signal = TriangleSignal(1100)
segment = signal.make_wave(duration=0.5, framerate=10000)
spectrum = segment.make_spectrum()
spectrum.plot()
decorate(xlabel='Frequency (Hz)')


#在idle或.py文件中使用如下方式输出图像
#将上面数据保存成图像sawtoothsignal.png
plt.savefig('sawtoothsignal.png')

2.5 计算频谱

make_spectrum的Wave方式，具体实现方式：

from np.fft import rfft, rfftfreq

#classWave:
    def mack_spectrum(self):
        n = len(self,ys)
        d = 1 / self.framerate

        hs = rfft(self.ys)
        fs = rfftfreq(n, d)
        
        return Spectrum(hs, fs, self.framerate)

参数 self 是一个Wave对象。 n 是波形的采样数，而 d 是采样频率的倒数，也就是样本间的时间差。

np.fft 是NumPy模块，提供快速傅里叶变换（FFT）的相关函数，这是一种计算离散傅里叶变化的高效算法。

make_spectrum 使用了 rfft ，其意义是“傅里叶变换市属部分（realFFT)”，因为Wave只有实数值，而没有复数值。之后会了解完整的FFT，它能处理复数信号（见7.9节），将rfft的结果赋值给hs，这是一个NumPy复数数值，表示的是波形中各个频率元素的振幅和相位差。

将 rfftfreq 的结果赋值给 fs ，这是一个包含对应 hs 各元素的频率的数组。

为了便于理解hs的值，可以参考以下两种理解复数的方式。

1. 复数是实数部分和虚数部分的和，通常写做 x + iy 的形式，其中 i 是虚数单元。

2. 复数还可以理解成一个量值与复指数的乘积，其中A 是量值，是弧度角，也被称为 “幅角”。可以在极角坐标系下理解A 和。

hs 中的每一个值对应一个频率元素。其量值与对应的元素的振幅成比例，其角度为相位差。

Spectrum 类提供了两个制度属性： amps 和 angles 。它们返回的是NumPy数组，代表的分别是 hs 的振幅和角度。在绘制 Spectrum 对象时，通常绘制的是 fs 对应的 amps 。有时绘制 fs 对应的 angles 也是有用的。

实际上，几乎完全不需要察看 hs 的实数和虚数部分。

可以将DFT想象为振幅和角度的向量，只不过恰好以复数的形式来表示了。

可以直接通过使用hs来修改Spetrum，例如：

spectrum.hs *= 2
spectrum.hs[spectrum.fs > cutoff] = 0

第一行将hs的元素乘以2，这将所有均匀的振幅都变成原来的两倍。第二行将hs元素中频率超过截至频率的所有元素都置零。

不过Spectrum还提供了如下的操作方式：

spectrum.scale(2)
spectrum.low_pass(cutoff)

可以在http://greenteapress.com/thinkdsp.html阅读以上所述方法相关文档。

2.5.1 振幅和相位

创建一个锯齿波信号

from thinkdsp import SawtoothSignal


signal = SawtoothSignal(500)
wave = signal.make_wave(duration=1, framerate=10000)
segment = wave.segment(duration=0.005)
segment.plot()
decorate(xlabel='Time (s)')


#在idle或.py文件中使用如下方式输出图像
#将上面数据保存成图像sawtoothsignal.png
plt.savefig('sawtoothsignal.png')

2.5.2 获取该波声音

wave.make_audio()


# 保存音频为WAV文件
from scipy.io.wavfile import write
scaled = np.int16(wave.ys/np.max(np.abs(wave.ys)) * 32767)
write('output.wav', wave.framerate, scaled)
 
# 现在，您可以使用音频播放器打开并播放output.wav文件。

2.5.3 计算hs——提取波形数组并计算实数FFT（这是一种针对实数输入优化的FFT）。

（jupyter notebook可以直接执行）

import numpy as np

hs = np.fft.rfft(wave.ys)
hs

返回结果：

2.5.4 计算fs ——计算与FFT元素相对应的频率。

n = len(wave.ys)                 # number of samples
d = 1 / wave.framerate           # time between samples
fs = np.fft.rfftfreq(n, d)
fs

返回结果：

2.5.5绘制振幅与频率的关系图

import matplotlib.pyplot as plt

magnitude = np.absolute(hs)
plt.plot(fs, magnitude)
decorate(xlabel='Frequency (Hz)')


#在idle或.py文件中使用如下方式输出图像
#将上面数据保存成图像testsignal_spectrum.png
plt.savefig('testsignal_spectrum.png')

2.5.6 绘制相位与频率的关系图

import numpy as np

angle = np.angle(hs)
plt.plot(fs, angle)
decorate(xlabel='Phase (radian)')

2.5.7 随机打乱相位

import random
random.shuffle(angle)
plt.plot(fs, angle)
decorate(xlabel='Phase (radian)')

#在idle或.py文件中使用如下方式输出图像
#将上面数据保存成图像random_phase.png
plt.savefig('random_phase.png')

将打乱的相位重新放回频谱中。hs中的每个元素都是一个复数，其幅度为A，相位为φ，我们可以用这个来计算Ae^(iφ)。

i = complex(0, 1)
spectrum = wave.make_spectrum()
spectrum.hs = magnitude * np.exp(i * angle)

2.5.8 使用irfft将频谱转换回波形。

wave2 = spectrum.make_wave()
wave2.normalize()
segment = wave2.segment(duration=0.005)
segment.plot()
decorate(xlabel='Time (s)')


#在idle或.py文件中使用如下方式输出图像
#将上面数据保存成图像wave.png
plt.savefig('wave.png')

2.5.9 获取该波形声音

输出output1.wav

wave2.make_audio()


# 保存音频为WAV文件
from scipy.io.wavfile import write
scaled = np.int16(wave2.ys/np.max(np.abs(wave2.ys)) * 32767)
write('output1.wav', wave2.framerate, scaled)
 
# 现在，您可以使用音频播放器打开并播放output1.wav文件。

获取原始波形用于比较，输出output.wav

wave.make_audio()


# 保存音频为WAV文件
from scipy.io.wavfile import write
scaled = np.int16(wave.ys/np.max(np.abs(wave.ys)) * 32767)
write('output.wav', wave.framerate, scaled)
 
# 现在，您可以使用音频播放器打开并播放output.wav文件。

尽管这两个信号的波形不同，但它们具有相同的频率成分和相同的振幅。它们只在相位上有所不同。

2.6 混叠相互作用

下面的交互程序探索了混叠对锯齿信号谐波的影响。

该代码只能在jupyter notebook中运行，运行后，随着对freq及framerate的调整，产生声音及频谱图，用于对比。

def view_harmonics(freq, framerate):
    """Plot the spectrum of a sawtooth signal.
    
    freq: frequency in Hz
    framerate: in frames/second
    """
    signal = SawtoothSignal(freq)
    wave = signal.make_wave(duration=0.5, framerate=framerate)
    spectrum = wave.make_spectrum()
    spectrum.plot(color='C0')
    decorate(xlabel='Frequency (Hz)', ylabel='Amplitude')
    display(wave.make_audio())

from ipywidgets import interact, interactive, fixed
import ipywidgets as widgets

slider1 = widgets.FloatSlider(min=100, max=10000, value=100, step=100)
slider2 = widgets.FloatSlider(min=5000, max=40000, value=10000, step=1000)
interact(view_harmonics, freq=slider1, framerate=slider2);

上一篇：《Python数字信号处理应用》学习笔记——第一章声音和信号-优快云博客

下一篇：《Python数字信号处理应用》学习笔记——第三章非周期信号