sprintf小数点后不补零_零、信号处理基础知识快速介绍（后）

Tux ZZ：零、信号处理基础知识快速介绍（前）​zhuanlan.zhihu.com

这篇文章主要讲述以下内容：

前：

1. 声音的时域表示
2. 正弦余弦信号
3. 声音的频域表示
4. 相关性、DFT、补零

后：

1. 窗函数、STFT
2. 滤波器

---

1. 窗函数、STFT

前篇中我们介绍了DFT的使用和实现。DFT输入一段时域信号，输出一段频域信号。但是直接使用DFT有两个问题：一、如果输入的信号在非整周期处被截断，DFT输出结果会出现能量泄漏。二、从频谱只能看出频率信息，难以看出时间信息。

　　问题一

截断信号，补零后（橘色）和补零前（蓝色）对比

　　上图是前文最后一个例子生成的图像。我们输入的信号为440Hz+880Hz两个频率，但是很明显在峰周围还有一些起伏，这种现象叫做能量泄漏，因为输入信号在非整周期处截断了。为了缓解这个问题，我们可以为输入信号乘一个中间大两边小的窗函数，使其截断处变得平滑些。

import numpy as np
import pylab as pl

if __name__ == "__main__":
  sr = 2000 # 使用2000Hz采样率
  n = 20 # 选个很短的长度，但是至少要2个周期
  freq_list = [440, 880] # 生成这两个频率的信号的叠加
  out = np.zeros(n, dtype=np.float32) # float32存储
  t = np.arange(n) / sr # 生成时间序列
  for freq in freq_list:
    out += 0.25 * np.sin(2 * np.pi * freq * t)
  
  # 生成一个窗，这里选hanning
  window = np.hanning(n)
  # 补零FFT，得到标准化后的振幅
  magn = np.abs(np.fft.rfft(out, n=n * 10)) * 2.0 / n
  # 加窗补零FFT，得到标准化后的振幅
  # 加窗后标准化系数变成了 2 / np.sum(window)
  magn_windowed = np.abs(np.fft.rfft(out * window, n=n * 10)) * 2.0 / np.sum(window)

  # 绘制并显示
  f = np.fft.rfftfreq(n * 10, d=1 / sr)

  pl.figure()
  pl.title("Time Domain")
  pl.plot(t, out, label="before")
  pl.plot(t, out * window, label="after")
  pl.plot(t, window, label="window")
  pl.legend()

  pl.figure()
  pl.title("DFT")
  pl.plot(f, magn, label="before")
  pl.plot(f, magn_windowed, label="after")
  pl.legend()

  pl.show()

加Hanning窗前后的时域和频域结果

　　从上图可明显看出，加窗后频谱峰周围的波动起伏变得小了不少，能量泄漏被减轻。但是也可以看出峰的宽度变得更宽，说明窗在提高信噪比的同时一定程度上降低了频率分辨能力。

　　窗的重要性质主要有三点：一、主瓣（Main-lobe）宽度：主瓣越宽，频率分辨能力越差。也叫扇形损失。二、旁瓣（Side-lobe）幅度：旁瓣幅度越小，信噪比越高，也叫噪声带宽、B值。三、第一旁瓣幅度：虽然有些窗旁瓣幅度衰减不如另外一些窗，但是其第一旁瓣幅度特别小，这是对分辨频谱中的峰值很有帮助的。

　　一个好的窗应该同时有窄主瓣和低旁瓣幅度，这样的窗无论是信噪比还是频率分辨能力都是极好的。但是鱼和熊掌并不能兼得，大部分主瓣特性好的窗旁瓣特性都不好，旁瓣特性好的主瓣特性不好，或者两个都不好（极少数情况特殊需要）。所以根据实际应用，对窗的选择做出取舍是很重要的。

大多数常见的窗函数的介绍都可以在维基百科找到。现在，我介绍几个我平时喜欢使用的窗：

　　* Hanning窗（也叫Hann窗）

主瓣宽度2.0，第一旁瓣-31.47dB，八度衰减18dB

　　这个窗因为性质均衡十分常用，遇到需要加窗的问题为什么不先试试它呢？

　　* Blackman窗

主瓣宽度3.03，第一旁瓣-58.11dB，八度衰减18dB

　　这个窗整体信噪比很好，又不像Blackman-Harris大幅增加主瓣宽度，是个比较常用的窗。

　　* Nuttall Minimum 3-Term窗

主瓣宽度3.0，第一旁瓣-70.83dB，八度衰减6dB

　　这是窗是Nuttall家族的一员，因为其不错的信噪比、主瓣宽度、极其漂亮的第一旁瓣，最近做谐波分析时我特别喜欢他。这个窗在Nuttall窗论文[1]中被提及。

　　* Boxcar窗

主瓣宽度1.0，旁瓣……

　　这个系数全为1的窗加了和没加有什么区别？要什么盒子车？其实这个窗就相当于没加窗，只不过在数学形式上，这个窗宽度之外全为0，达到把信号硬生生截断的效果。这个窗一般不会特意去用，因为它的旁瓣特性实在是太差了。这里单独拎出来讲是因为它很特殊。

　　任何一个窗的特性都可以用这段代码去检查。这段代码也附赠了一个通用的Nuttall窗生成器，只需要把Nuttall论文[1]中的系数抄下来，就可以生成相应的窗了，用法和scipy.signal（我喜欢简称sp）的窗函数们类似。

import numpy as np
import scipy.signal as sp
import pylab as pl

nuttall_min3_coeff = np.array([0.4243801, 0.4973406, 0.0782793], dtype=np.float32)
def nuttall(M, coeff, sym):
    if M < 1:
        return np.array([], dtype=np.float32)
    if M == 1:
        return np.ones(1, dtype=np.float32)
    odd = bool(M % 2)
    if not sym and not odd:
        M = M + 1
    n = np.arange(0, M)
    sign_list = [-1.0, 1.0]
    w = coeff[0]
    for (i, x) in enumerate(coeff[1:]):
        sign = sign_list[i % 2]
        w += sign * x * np.cos((i + 1) * 2.0 * np.pi * n / (M - 1))
    if not sym and not odd:
        w = w[:-1]
    return w

def test_window(x):
  (n,) = x.shape
  log_magn = np.log(np.abs(np.fft.rfft(x, n=n * 1000)) * 2 / n + 1e-8)
  pl.subplot(121)
  pl.title("Window")
  pl.plot(x)
  pl.subplot(122)
  pl.title("Window Lobes")
  pl.plot(np.arange(n * 1000 // 2 + 1) / 1000, log_magn)
  pl.show()

test_window(nuttall(100, nuttall_min3_coeff, sym=True))

　　问题二

　　前面我们缓解了DFT能量泄漏的问题。DFT还有另一个问题：从频谱只能看出频率信息，难以看出时间信息。为了能轻松地同时看出时间和频率信息，STFT出现了。其思想十分简单：将原始信号切成一个个小窗口，窗口按一定步幅沿着时间轴滑动，将每个窗口的时域信号进行DFT，得到时频谱（spectrogram，中文语境经常直接简称频谱）

　　在Audacity中，我们可以在轨道上切换查看频谱。在之前我们生成的440Hz+880Hz+1320Hz三成分信号的频谱中，我们可以轻松地看到三组信号。你也可以拖进一首音乐，然后查看它的频谱。如果你觉得频谱很模糊，你可以试着在设置中调大FFT尺寸，这样可以提高频率分辨率，但是会损失时间分辨率。

切换到频谱

　　我们也可以在Python中实现STFT：

import numpy as np
import pylab as pl
import scipy.io.wavfile as wavfile

def ensureEven(x: int):
  # 确保x为偶数
  return x + 1 if x % 2 == 1 else x

def load_wav(path):
  # wav加载函数
  samprate, w = wavfile.read(path)
  if(w.dtype == np.int8):
    w = w.astype(np.float32) / 127.0
  elif(w.dtype == np.short):
    w = w.astype(np.float32) / 32767.0
  elif(w.dtype == np.int32):
    w = w.astype(np.float32) / 2147483647.0
  elif(w.dtype == np.float32):
    w = w.copy()
  elif(w.dtype == np.float64):
    w = w.astype(np.float32)
  else:
    raise ValueError("Unsupported sample format: %s" % (str(w.dtype)))
  return w, samprate

def roundUpToPowerOf2(v):
  # 返回至少为v的2的整次幂
  return int(2 ** np.ceil(np.log2(v)))

def getNFrame(inputSize, hopSize):
  # 获得帧数/跳数
  return int(round(inputSize / hopSize))

def getFrame(input, center, size):
  # 获得一帧/一跳
  out = np.zeros(size, input.dtype)
  outBegin, outEnd, inputBegin, inputEnd = getFrameRange(len(input), center, size)
  out[outBegin:outEnd] = input[inputBegin:inputEnd]
  return out

def getFrameRange(inputLen, center, size):
  # 获得帧范围
  leftSize = size // 2
  rightSize = size - leftSize
  inputBegin = min(inputLen, max(center - leftSize, 0))
  inputEnd = max(0, min(center + rightSize, inputLen))
  outBegin = max(leftSize - center, 0)
  outEnd = outBegin + (inputEnd - inputBegin)
  return outBegin, outEnd, inputBegin, inputEnd

def stft(x, window, hop_size, fft_size):
    # 获得一些常量
    (window_size,) = window.shape
    (n_x,) = x.shape
    n_hop = getNFrame(n_x, hop_size)
    n_bin = fft_size // 2 + 1

    # 进行一些必要检查
    assert window_size <= fft_size
    assert window_size % 2 == 0
    assert hop_size > 0
    # FFT只有在2的整次幂时才能使用
    assert fft_size > 0 and roundUpToPowerOf2(fft_size) == fft_size

    # 计算标准化过的窗
    window_norm_fac = 2.0 / np.sum(window)
    window = window * window_norm_fac
    del window_norm_fac

    # 取出一个个帧，做FFT
    out = np.zeros((n_hop, n_bin), dtype=np.complex64)
    for i_hop in range(n_hop):
      i_center = int(round(i_hop * hop_size))
      frame = getFrame(x, i_center, window_size)
      frame *= window
      out[i_hop] = np.fft.rfft(frame, n=fft_size)
    return out

if __name__ == "__main__":
  w, sr = load_wav("suiheisen.wav")
  w = np.atleast_2d(w.T)[0] # 无论有几个声道，只要第一个
  window = np.hanning(ensureEven(int(sr * 0.05))) # 50ms Hanning窗
  hop_size = sr * 0.005 # 5ms 帧间距

  # 计算STFT，取对数模得到幅度
  spec = stft(w, window, hop_size, roundUpToPowerOf2(window.size))
  magn = np.log(np.abs(spec) + 1e-5) # 取log前加一个较小的数防止除零

  # 绘图
  pl.imshow(magn.T, interpolation="nearest", aspect="auto", origin="lower", extent=[0.0, w.size / sr, 0.0, sr / 2])
  pl.show()

运行结果

　　因为与物理学中测不准原理类似的原因，无论是使用STFT、还是更高级的小波变换，频率分辨率和时间分辨率是无法兼得的。应该根据需求和信号本身选择合适的窗口大小。需要时间分辨率高时（比如识别快速变化的信号）使用小窗口，需要频率分辨率高时（比如需要的信号中混合了大量其他噪声，或者信号足够平稳）使用大窗口。

2.滤波器

使用滤波器我们可以消除不想要的频率成分。这里我们对FIR（有限脉冲响应）滤波器和IIR（无限脉冲响应）滤波器这两种常用的滤波器仅仅做一下简单介绍。

　　脉冲响应的连续表示是狄拉克δ函数。但是在离散+带限（有限采样率）情况下，其时域一个Sinc函数，变换到频域是一个Boxcar函数。将这样一个信号输入进FIR或者IIR滤波器，可以测试滤波器的特性。

　　FIR滤波器在同样阶数时频率衰减性能通常不如IIR滤波器，其优势是设计容易、不会有不稳定问题、且可以保证线性相位。IIR滤波器设计相对困难，容易不稳定、通常不是线性相位，但是阶数相同时比IIR频率衰减性能好。随着电脑计算速度的飞速发展，大多数情况可以无脑上高阶FIR滤波器。对于计算性能受限、需要超低延迟（阶数越高的滤波器需要引入更高的延迟）、或者不需要那么好的相位特性的地方，IIR滤波器也是个很好的选择。IIR也可以使用全通滤波器或者filtfilt等工具校准相位特性，很多时候是够用的。

　　通过设计得到一套FIR滤波器系数后，可以直接利用和原信号卷积的方式进行滤波。IIR滤波器拿到的系数则需要用差分公式（或者lfilter一类的包装好的函数）去使用了。

　　设计数字滤波器通常使用Z变换。FIR滤波器可以直接用Sinc函数加窗实现。

　　具体可以参考Scipy的firwin（设计FIR滤波器的函数）和butter（设计巴特沃夫IIR滤波器的函数）的文档。

[1] Nuttall, Albert. "Some windows with very good sidelobe behavior." IEEE Transactions on Acoustics, Speech, and Signal Processing 29.1 (1981): 84-91.

下一篇文章将开始讲述语音信号的基本结构，并介绍这些年来常用的分析和合成方法。

再下一篇开始基频分析。