频域信号处理

代码来源于

　　文章从6个方面来写，首先是观察频谱的特征，第二部分是加上窗函数之后的特征，第三部分是频谱平均，第四部分是比较fft与直接卷积时间效率区别，第五部分是由于fft对输入信号的长度有要求，因此介绍了overlap-add分段运算，最后一部分是hilbert变换的实现。

• 观察信号的频谱

　　数据通过fft转换成频域信号，对频域信号进行分析，再通过ifft转换成时域信号。

import numpy as np
import pylab as pl
import matplotlib as mpl
mpl.rcparams['font.sans-serif'] = ['kaiti']
mpl.rcparams['font.serif'] = ['kaiti']
mpl.rcparams['axes.unicode_minus']=false
sampling_rate = 8000   #取样频率
fft_size = 512            #fft长度
t = np.arange(0, 1.0,  1.0/sampling_rate)   
#假设取样频率为fs, 取波形中的n个数据进行fft变换。那么这n点数据包含整数个周期的波形时，fft所计算的结果是精确的。于是能精确计算的波形的周期是: n*fs/n。
#对于8khz取样，512点fft来说，8000/512.0 = 15.625hz，前面的156.25hz和234.375hz正好是其10倍和15倍。
#选取整数倍的数据，查看当fft后的数据在频谱中形成整数周期时的情况。
x = np.sin(2*np.pi*156.25*t)  + 2*np.sin(2*np.pi*234.375*t)
#选取非整数倍的数据，查看当fft后的数据在频谱中没有形成非整数周期时的情况。
x = np.sin(2*np.pi*200*t)  + 2*np.sin(2*np.pi*300*t)
xs = x[:fft_size]  #取数据
xf = np.fft.rfft(xs)/fft_size    #rfft：对实数信号进行fft变换。/fft_size是为了正确显示波形能量
freqs = np.linspace(0, sampling_rate/2, fft_size/2+1) #fft_size/2+1个点，后面的是与前面的共轭
#计算每个频率分量的幅值，并通过 20*np.log10() 将其转换为以db(分贝)为单位的值
xfp = 20*np.log10(np.clip(np.abs(xf), 1e-20, 1e100))  #clip将数据限制在最小值和最大值之间
pl.figure(figsize=(8,4))   #新建一个8*4英寸的图纸
pl.subplot(211)   #绘制2行1列的图纸，这个图形占据第一行
pl.plot(t[:fft_size], xs)
pl.xlabel(u"时间(秒)")
pl.title(u"156.25hz和234.375hz的波形和频谱")
pl.subplot(212)     #绘制2行1列的图纸，这个图形占据第二行
pl.plot(freqs, xfp)
pl.xlabel(u"频率(hz)")
pl.subplots_adjust(hspace=0.4)
pl.show()

　　运行结果如下图所示：

　　整数周期情况：

 

　　非整数周期情况：

　

如上两图可看出非整数周期情况下，第二种情况可能会发生信号泄露状况。由于fft的假设前提是测量之外的信号是所测量信号的不断重复，如代码中那样，取8k个样例，从信号中取出的512个数据就是fft的测量范围，也就是说每个波形周期是15.625hz，它计算的是这512个数据一直重复的波形的频谱。显然如果512个数据包含整数个周期的话，那么得到的结果就是原始信号的频谱，而如果不是整数周期的话，得到的频谱就是如下波形的频谱，这里假设对50hz的正弦波进行512点fft，这种波形会发生跳变。

• 窗函数

　　为了减少fft所截取的数据段前后的跳变，可以对数据先乘以一个窗函数，使得其前后数据能平滑过渡。例如常用的hann窗函数的定义如下：

hann窗的曲线如下所示：

可看到，最前和最后的值都是0，如果直接去乘的前后会出现两个0，因此可考虑将这个波形用n+1个点表示，而取前n个点，这样第n+1个点就是下一个波形的第一个点，也就是0,通过设置sym参数解决。这与调用linspace时指定endpoint=false类似，丢掉最后一个点。

signal.hann(8)
#array([0.        , 0.1882551 , 0.61126047, 0.95048443, 0.95048443,
#       0.61126047, 0.1882551 , 0.        ])
signal.hann(8, sym=0)
#array([0.        , 0.14644661, 0.5       , 0.85355339, 1.        ,
#       0.85355339, 0.5       , 0.14644661])

 将hann窗与50hz相乘，它的曲线会更加平滑。

之前的非整数周期加了hann窗之后的结果如下图所示：

• 频谱平均

　　对于频谱特性不随时间变化的信号，例如引擎、压缩机等机器噪声，可以对其进行长时间的采样，然后分段进行fft计算，最后对每个频率分量的幅值求其平均值可以准确地测量信号的频谱，测试随机数序列频谱如下所示

def average_fft(x, fft_size):
    n = len(x) // fft_size * fft_size
    tmp = x[:n].reshape(-1, fft_size)  #将n行转成n列，转为一个二维数组
    tmp *= signal.hann(fft_size, sym=0)
    xf = np.abs(np.fft.rfft(tmp)/fft_size)
    avgf = np.average(xf, axis=0)
    return 20*np.log10(avgf)

x = np.random.rand(100000) - 0.5
xf = average_fft(x, 512)
pl.plot(xf)
pl.show()

 结果为

这个频谱的能量趋近于一条直线，每个窗的能量相差不大，被称为白色噪声。

•  快速卷积

　　信息处理可看作是将原始信号与一个信号进行卷积，也就需要考虑运算效率。这部分主要是比较fft和直接卷积的运算效率

def fft_convolve(a,b):
    n = len(a)+len(b)-1
    n = 2**(int(np.log2(n))+1)          #fft的长度：大于n的最小的2的整数次幂
    a = np.fft.fft(a, n)
    b = np.fft.fft(b, n)
    return np.fft.ifft(a*b)[:n]

if __name__ == "__main__":
    a = np.random.rand(128)
    b = np.random.rand(128)
    c = np.convolve(a,b)
    
    print (np.sum(np.abs(c - fft_convolve(a,b))))

 结果为1.865645261656436e-12，也就是说fft和普通卷积的结果相差很小，但速度却快很多。

• 分段运算

　　对于输入信号 x 和系统向量（eg:fir滤波器）h而言，x的长度不固定，h的长度固定。为了加快卷积效率， 我们需要x和h的长度相当，也就是说对x进行分段处理，这种分段算法被称为overlap-add运算。但是由于fft在两个数组的分段长度相当时最为有效，因此在实时性要求很强的系统中，采用直接卷积会更好一些。

 

import numpy as np
x = np.random.rand(1000)
h = np.random.rand(101)
y = np.convolve(x, h)

n = 50 # 分段大小
m = len(h) # 滤波器长度

output = []

#缓存初始化为0
buffer = np.zeros(m+n-1,dtype=np.float64)

for i in range(int(len(x)/n)):
    #从输入信号中读取n个数据
    xslice = x[i*n:(i+1)*n]
    #计算卷积
    yslice = np.convolve(xslice, h)
    pl.cla()
    pl.plot(yslice)
    #将卷积的结果加入到缓冲中
    buffer += yslice
    #输出缓存中的前n个数据，注意使用copy，否则输出的是buffer的一个视图
    output.append( buffer[:n].copy() )
    #缓存中的数据左移动n个元素
    buffer[0:m-1] = buffer[n:]
    #后面的补0
    buffer[m-1:] = 0

#将输出的数据组合为数组
y2 = np.hstack(output)
#计算和直接卷积的结果之间的误差
print (np.sum(np.abs( y2 - y[:len(x)] ) ))

 

• hilbert 变换

　　hilbert变换能在振幅保持不变的情况下将输入信号的相角偏移90度，简单地说就是能将正弦波形转换为余弦波形。

from scipy import fftpack
import numpy as np
import matplotlib.pyplot as pl

# 产生1024点4个周期的正弦波
t = np.linspace(0, 8*np.pi, 1024, endpoint=false)
x = np.sin(t)

# 进行hilbert变换
y = fftpack.hilbert(x)
pl.plot(x, label=u"原始波形")
pl.plot(y, label=u"hilbert转换后的波形")
pl.legend()
pl.show()

结果如下所示：

　　hilbert变换后可将直流分量变为0，正频率成分偏移+90度，负频率成分偏移-90度。它也可用来进行包络检波。

import numpy as np
import pylab as pl
from scipy import fftpack

t = np.arange(0, 0.3, 1/20000.0)
x = np.sin(2*np.pi*1000*t) * (np.sin(2*np.pi*10*t) + np.sin(2*np.pi*7*t) + 3.0)
hx = fftpack.hilbert(x)

pl.plot(x, label=u"载波信号")
pl.plot(np.sqrt(x**2 + hx**2), "r", linewidth=2, label=u"检出的包络信号")
pl.title(u"使用hilbert变换进行包络检波")
pl.legend()
pl.show()