有关傅⾥叶变换的知识整理
有关傅⾥叶变换的知识整理
傅⾥叶变换的含义
傅⾥叶变换是信号领域沟通时域和频域的桥梁，在频域⾥可以更⽅便的进⾏⼀些分析。傅⾥叶主要针对的是平稳信号的频率特性分析，简单说就是具有⼀定周期性的信号，因为傅⾥叶变换采取的是有限取样的⽅式，所以对于取样长度和取样对象有着⼀定的要求。
快速傅⾥叶变换FFT
1.假设采样频率为Fs，信号频率fs，采样点数为N。那么FFT之后结果就是⼀个为N点的复数。每⼀个点就对应着⼀个频率点。这个点的模值，就是该频率值下的幅度特性
假设FFT之后某点n⽤复数a+bi表⽰，那么这个复数的模就是An=sqrt（aa+bb）（某点处的幅度值An = A*(N/2)
采样频率： 采样频率，也称为采样速度或者采样率，定义了每秒从连续信号中提取并组成离散信号的采样个数，它⽤赫兹（Hz）来表⽰。采样定理： 采样定理指出，如果信号是带限的，并且采样频率⾼于信号带宽的两倍，那么，原来的连续信号可以从采样样本中完全重建出来。
在进⾏模拟/数字信号的转换过程中，当采样频率fs⼤于信号中最⾼频率fmax的2倍时，即fs>2*fmax⼀般实际应⽤中保证采样频率为信号最⾼频率的2.56～4倍；采样定理⼜称奈奎斯特定理。
2.FFT得到的复数的模（即绝对值）就是对应的“振幅谱”，复数所对应的⾓度，就是所对应的“相位谱”
在傅⾥叶分析中，把各个分量的幅度|Fn|或 Cn 随着频率nω1的变化称为信号的幅度谱。⽽把各个分量的相位 φn 随⾓频率 nω1 变化称为信号的相位谱。幅度谱和相位谱统称为信号的频谱。
三⾓形式的傅⾥叶级数频率为⾮负的，对应的频谱⼀般称为单边谱；指数形式的傅⾥叶级数频率为整个实轴，所以称为双边谱。
3. 将振幅谱进⾏归⼀化和取半处理
考虑到数量级较⼤，⼀般进⾏归⼀化处理，既然第⼀个峰值是A1的N倍，那么将每⼀个振幅值都除以N即可；FFT具有对称性，⼀般只需要⽤N的⼀半，前半部分即可。
完整代码
import numpy as np
from scipy.fftpack import fft,ifft
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl
#采样点选择1400个，因为设置的信号频率分量最⾼为600赫兹，根据采样定理知采样频率要⼤于信号频率2倍，所以这⾥设置采样频率为1400赫兹（即⼀秒内有1400个采样点，⼀样意思的）
x=np.linspace(0,1,1400)
#设置需要采样的信号，频率分量有200，400和600
y=7*np.sin(2*np.pi*200*x)+5*np.sin(2*np.pi*400*x)+3*np.sin(2*np.pi*600*x)
fft_y=fft(y)                          #快速傅⾥叶变换
N=1400
x = np.arange(N)            # 频率个数
half_x = x[range(int(N/2))]  #取⼀半区间
abs_y=np.abs(fft_y)                # 取复数的绝对值，即复数的模(双边频谱)
angle_y=np.angle(fft_y)            #取复数的⾓度
normalization_y=abs_y/N            #归⼀化处理（双边频谱）
normalization_half_y = normalization_y[range(int(N/2))]      #由于对称性，只取⼀半区间（单边频谱）
plt.subplot(231)
plt.plot(x,y)
plt.title('原始波形')
plt.subplot(232)
plt.plot(x,fft_y,'black')
plt.title('双边振幅谱(未求振幅绝对值)',fontsize=9,color='black')
plt.subplot(233)
plt.plot(x,abs_y,'r')
plt.title('双边振幅谱(未归⼀化)',fontsize=9,color='red')
plt.subplot(234)
plt.plot(x,angle_y,'violet')
plt.title('双边相位谱(未归⼀化)',fontsize=9,color='violet')
plt.subplot(235)
plt.plot(x,normalization_y,'g')
plt.title('双边振幅谱(归⼀化)',fontsize=9,color='green')
plt.subplot(236)
plt.plot(half_x,normalization_half_y,'blue')
plt.title('单边振幅谱(归⼀化)',fontsize=9,color='blue')
plt.show()
import numpy as np#导⼊⼀个数据处理模块
import pylab as pl#导⼊⼀个绘图模块，matplotlib下的模块
sampling_rate =8000#采样频率为8000Hz
fft_size =512 #FFT处理的取样长度
t = np.arange(0,1.0,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
# 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)的分。
matplotlib中subplot#于是可以通过下⾯的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,1e100))
#最后我们计算每个频率分量的幅值，并通过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"Time(S)")
pl.title(u"156.25Hz and 234.375Hz WaveForm And Freq")
pl.subplot(212)
pl.plot(freqs, xfp)
pl.xlabel(u"Freq(Hz)")
pl.subplots_adjust(hspace=0.4)
pl.show()
信号的四种频率特性：频谱、频谱密度、能量谱密度、功率谱密度 Energy Spectra, Power Spectra
功率谱
功率谱是功率谱密度函数（PSD）的简称，它定义为单位频带内的信号功率。
功率谱是针对功率信号来说的。功率谱的推导公式相对复杂，不过幸运的是维纳-⾟钦定理证明了：⼀段信号的功率谱等于这段信号⾃相关函数的傅⾥叶变换。
所以求功率谱就有了两种⽅法：1.(傅⽴叶变换的平⽅)/(区间长度)；2.⾃相关函数的傅⾥叶变换。这两种⽅法分别叫做直接法和相关函数法。
功率谱表⽰了信号功率随着频率的变化关系。常⽤于功率信号(区别于能量信号)的表述与分析，其曲线(即功率谱曲线)⼀般横坐标为频率，纵坐标为功率。由于功率没有负值，所以功率谱曲线上的纵坐标也没有负数值，功率谱曲线所覆盖的⾯积在数值上等于信号的总功率(能量)
代码实现
from scipy.fftpack import fft, fftshift, ifft
from scipy.fftpack import fftfreq
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")
fs =1000
#采样点数
num_fft =1024;
"""
⽣成原始信号序列
在原始信号中加上噪声
np.random.randn(t.size)
"""
t = np.arange(0,1,1/fs)
f0 =100
f1 =200
x = np.cos(2*np.pi*f0*t)+s(2*np.pi*f1*t)+ np.random.randn(t.size)
plt.figure(figsize=(15,12))
ax=plt.subplot(511)
ax.set_title('original signal')
plt.tight_layout()
plt.plot(x)
"""
FFT(Fast Fourier Transformation)快速傅⾥叶变换
"""
Y=fft(x, num_fft)
Y= np.abs(Y)
ax=plt.subplot(512)
ax.set_title('fft transform')
plt.plot(20*np.log10(Y[:num_fft//2]))
"""
功率谱 power spectrum
直接平⽅
"""
ps =Y**2/ num_fft
ax=plt.subplot(513)
ax.set_title('direct method')
plt.plot(20*np.log10(ps[:num_fft//2]))
"""
相关功谱率 power spectrum using correlate
间接法
"""
cor_x = np.correlate(x, x,'same')
cor_X =fft(cor_x, num_fft)
ps_cor = np.abs(cor_X)
ps_cor = ps_cor / np.max(ps_cor)
ax=plt.subplot(514)
ax.set_title('indirect method')
plt.plot(20*np.log10(ps_cor[:num_fft//2]))
plt.tight_layout()
plt.show()

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。