# 频域信号处理

## 观察信号的频谱

# -*- coding: utf-8 -*-
import numpy as np
import pylab as pl

sampling_rate = 8000
fft_size = 512
t = np.arange(0, 1.0, 1.0/sampling_rate)
x = np.sin(2*np.pi*156.25*t)  + 2*np.sin(2*np.pi*234.375*t)
xs = x[:fft_size]
xf = np.fft.rfft(xs)/fft_size
freqs = np.linspace(0, sampling_rate/2, fft_size/2+1)
xfp = 20*np.log10(np.clip(np.abs(xf), 1e-20, 1e100))
pl.figure(figsize=(8,4))
pl.subplot(211)
pl.plot(t[:fft_size], xs)
pl.xlabel(u"时间(秒)")
pl.title(u"156.25Hz和234.375Hz的波形和频谱")
pl.subplot(212)
pl.plot(freqs, xfp)
pl.xlabel(u"频率(Hz)")
pl.subplots_adjust(hspace=0.4)
pl.show()


t = np.arange(0, 1.0, 1.0/sampling_rate)


x = np.sin(2*np.pi*156.25*t)  + 2*np.sin(2*np.pi*234.375*t)


N点FFT能精确计算的频率

xs = x[:fft_size]
xf = np.fft.rfft(xs)/fft_size


rfft函数的返回值是N/2+1个复数，分别表示从0(Hz)到sampling_rate/2(Hz)的N/2+1点频率的成分。于是可以通过下面的np.linspace计算出返回值中每个下标对应的真正的频率：

freqs = np.linspace(0, sampling_rate/2, fft_size/2+1)


xfp = 20*np.log10(np.clip(np.abs(xf), 1e-20, 1e100))


>>> xfp[10]
-6.0205999132796251
>>> xfp[15]
-9.6432746655328714e-16


x = np.sin(2*np.pi*200*t)  + 2*np.sin(2*np.pi*300*t)


>>> t = np.arange(0, 1.0, 1.0/8000)
>>> x = np.sin(2*np.pi*50*t)[:512]
>>> pl.plot(np.hstack([x,x,x]))
>>> pl.show()


50Hz正弦波的512点FFT所计算的频谱的实际波形

### 窗函数

>>> import pylab as pl
>>> import scipy.signal as signal
>>> pl.figure(figsize=(8,3))
>>> pl.plot(signal.hann(512))


hann窗函数

>>> np.sin(np.arange(0, 2*np.pi, 2*np.pi/10))
array([  0.00000000e+00,   5.87785252e-01,   9.51056516e-01,
9.51056516e-01,   5.87785252e-01,   1.22464680e-16,
-5.87785252e-01,  -9.51056516e-01,  -9.51056516e-01,
-5.87785252e-01])
>>> np.sin(np.linspace(0, 2*np.pi, 10))
array([  0.00000000e+00,   6.42787610e-01,   9.84807753e-01,
8.66025404e-01,   3.42020143e-01,  -3.42020143e-01,
-8.66025404e-01,  -9.84807753e-01,  -6.42787610e-01,
-2.44929360e-16])


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


50Hz正弦波与窗函数乘积之后的重复波形如下：

>>> t = np.arange(0, 1.0, 1.0/8000)
>>> x = np.sin(2*np.pi*50*t)[:512] * signal.hann(512, sym=0)
>>> pl.plot(np.hstack([x,x,x]))
>>> pl.show()


>>> np.sum(signal.hann(512, sym=0))/512
0.5


### 频谱平均

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

def average_fft(x, fft_size):
n = len(x) // fft_size * fft_size
tmp = x[:n].reshape(-1, fft_size)
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)


average_fft(x, fft_size)对数组x进行fft_size点FFT运算，以dB为单位返回其平均后的幅值。由于x的长度可能不是fft_size的整数倍，因此首先将其缩短为fft_size的整数倍，然后用reshape函数将其转换为一个二维数组tmp。tmp的第1轴的长度为fft_size：

n = len(x) // fft_size * fft_size
tmp = x[:n].reshape(-1, fft_size)


tmp *= signal.hann(fft_size, sym=0)


xf = np.abs(np.fft.rfft(tmp)/fft_size)


avgf = np.average(xf, axis=0)


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


>>> b,a=signal.iirdesign(1000/4000.0, 1100/4000.0, 1, -40, 0, "cheby1")
>>> x = np.random.rand(100000) - 0.5
>>> y = signal.filtfilt(b, a, x)


## 快速卷积

# -*- coding: utf-8 -*-
import numpy as np

def fft_convolve(a,b):
n = len(a)+len(b)-1
N = 2**(int(np.log2(n))+1)
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)))


N = 2**(int(np.log2(n))+1)


>>> import timeit
>>> setup="""import numpy as np
a=np.random.rand(10000)
b=np.random.rand(10000)
from spectrum_fft_convolve import fft_convolve"""
>>> timeit.timeit("np.convolve(a,b)",setup, number=10)
1.852900578146091
>>> timeit.timeit("fft_convolve(a,b)",setup, number=10)
0.19475575806416145


### 分段运算

overlap-add的计算方法如下图所示：

1. 建立一个缓存，其大小为N+M-1，初始值为0
2. 每次从x中读取N个数据，和h进行卷积，得到N+M-1个数据，和缓存中的数据进行求和，并放进缓存中，然后输出缓存前N个数据
3. 将缓存中的数据向左移动N个元素，也就是让缓存中的第N个元素成为第0个元素，后面的N个元素全部设置为0
4. 跳转到2重复运行

# -*- coding: utf-8 -*-
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 xrange(len(x)/N):
#从输入信号中读取N个数据
xslice = x[i*N:(i+1)*N]
#计算卷积
yslice = np.convolve(xslice, h)
#将卷积的结果加入到缓冲中
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度，简单地说就是能将正弦波形转换为余弦波形，下面的程序验证这一特性：

# -*- coding: utf-8 -*-
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转换函数在scipy.fftpack函数库中
• 为了生成完美的正弦波，需要计算整数个周期，因此调用linspace时指定endpoint=False，这样就不会包括区间的结束点：8*np.pi

Hilbert变换将正弦波变为余弦波

Hilbert正变换的相角偏移符号

• 直流分量为0
• 正频率成分偏移+90度
• 负频率成分偏移-90度

>>> x = np.random.rand(16)
>>> y = fftpack.hilbert(x)
>>> X = np.fft.fft(x)
>>> Y = np.fft.fft(y)
>>> np.imag(Y/X)
array([ 0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,
0., -1., -1., -1., -1., -1., -1., -1.])


• 下标为0的频率分量表示直流分量
• 下标为N/2的的频率分量为取样频率/2的频率分量
• 1到N/2-1为正频率分量
• N/2+1到N为负频率分量

Hilbert变换可以用作包络检波。具体算法如下所示：

# -*- coding: utf-8 -*-
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()


>>> run filter_lfilter_example01.py # 运行滤波器的例子
>>> hy = fftpack.hilbert(y)
>>> pl.plot( np.sqrt(y**2 + hy**2),"r", linewidth=2)
>>> pl.plot(y)
>>> pl.title(u"频率扫描波的包络")
>>> pl.show()