【上海校区】频域信号处理 - 黑马程序员技术交流社区

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

[Python] 纯文本查看 复制代码

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类似，丢掉最后一个点。

[Python] 纯文本查看 复制代码

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计算，最后对每个频率分量的幅值求其平均值可以准确地测量信号的频谱，测试随机数序列频谱如下所示

[Python] 纯文本查看 复制代码

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)

结果为

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

快速卷积

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

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

[Python] 纯文本查看 复制代码

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在两个数组的分段长度相当时最为有效，因此在实时性要求很强的系统中，采用直接卷积会更好一些。

[Python] 纯文本查看 复制代码

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

[Python] 纯文本查看 复制代码

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度。它也可用来进行包络检波。

[Python] 纯文本查看 复制代码

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()