来源:EETOP 行者无疆(论坛usrname:ICNO.1) 的博客 地址:http://www./blog/?xhsir520 Python 是目前的热门语言,一直觉得掌握一门编程语言对作为搞技术的来说还是很有必要的,结合工作中能用到的一些数据处理和分析的内容,觉得从数据分析入手,争取能够掌握Python在数据处理领域的一些应用。下面是基于Python的numpy进行的数字信号的频谱分析介绍 一、傅里叶变换 傅里叶变换是信号领域沟通时域和频域的桥梁,在频域里可以更方便的进行一些分析。傅里叶主要针对的是平稳信号的频率特性分析,简单说就是具有一定周期性的信号,因为傅里叶变换采取的是有限取样的方式,所以对于取样长度和取样对象有着一定的要求。 二、基于Python的频谱分析 # _*_ coding:utf-8 _*_ import numpy as np #导入一个数据处理的模块 import pylab as pl #导入一个绘图模块,matplotlib下的模块 sampling_rate = 8000 ##取样频率 fft_size =512 #FFT处理的取样长度 t = np.arange(0,1.1,1.0/sampling_rate) #np.arange(起点,终点,间隔)产生1s长的取样时间 x = np.sin(2*np.pi*156.25*t)+2*np.sin(2*np.pi*234.375*t) #两个正弦波叠加,156.25HZ和234.375HZ,因此如上面简单 #的介绍FFT对于取样时间有要求, #N点FFT进行精确频谱分析的要求是N个取样点包含整数个 #取样对象的波形。 #因此N点FFT能够完美计算频谱对取样对象的要求 #是n*Fs/N(n*采样频率/FFT长度), #因此对8KHZ和512点而言, #完美采样对象的周期最小要求是8000/512=15.625HZ, #所以156.25的n为10,234.375的n为15。 xs = x[:fft_size]# 从波形数据中取样fft_size个点进行运算 xf = np.fft.rfft(xs)/fft_size # 利用np.fft.rfft()进行FFT计算,rfft()是为了更方便 #对实数信号进行变换,由公式可知/fft_size为了正确显示波形能量 # rfft函数的返回值是N/2+1个复数,分别表示从0(Hz) #到sampling_rate/2(Hz)的分。 #于是可以通过下面的np.linspace计算出返回值中每个下标对应的真正的频率: freqs = np.linspace(0,sampling_rate/2,fft_size/2+1) # np.linspace(start, stop, num=50, endpoint=True, retstep=False, dtype=None) #在指定的间隔内返回均匀间隔的数字 xfp = 20*np.log10(np.clip(np.abs(xf),1e-20,1e1000)) #最后我们计算每个频率分量的幅值,并通过 20*np.log10() #将其转换为以db单位的值。为了防止0幅值的成分造成log10无法计算, #我们调用np.clip对xf的幅值进行上下限处理 pl.figure(figsize=(8,4)) pl.subplot(211) pl.plot(t[:fft_size], xs) pl.xlabel(u"时间(秒)") pl.title(u"The Wave and Spectrum 156.25Hz234.375Hz") pl.subplot(212) pl.plot(freqs, xfp) pl.xlabel(u"Hz") pl.subplots_adjust(hspace=0.4) pl.show() #绘图显示结果 现在来看看频谱泄露,将采样对象的频率改变 x = np.sin(2*np.pi*100*t)+2*np.sin(2*np.pi*234.375*t) 我们明显看出,第一个对象的频谱分析出现“泄露”,能量分散到其他频率上, 没法准确计算到计算对象的频谱特性。 窗函数 上面我们可以看出可以通过加“窗”函数的方法来处理,尽量保证FFT长度内 的取样对象是对称的。 import pylab as pl import scipy.signal as signal pl.figure(figsize=(8,3)) pl.plot(signal.hann(512))#汉明窗函数 pl.show() 对上述出现频谱泄露的函数进行加窗处理,后面会介绍一下各种加窗函数的原理和效果。 |
|
来自: 毕杰lb7q1kq7pr > 《待分类》