有关傅里叶变换的知识整理_自相关函数的傅里叶变换

有关傅里叶变换的知识整理

傅里叶变换的含义

傅里叶变换是信号领域沟通时域和频域的桥梁，在频域里可以更方便的进行一些分析。傅里叶主要针对的是平稳信号的频率特性分析，简单说就是具有一定周期性的信号，因为傅里叶变换采取的是有限取样的方式，所以对于取样长度和取样对象有着一定的要求。

快速傅里叶变换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
 
mpl.rcParams['font.sans-serif'] = ['SimHei']   #显示中文
mpl.rcParams['axes.unicode_minus']=False       #显示负号
 
#采样点选择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()
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

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)的分。
#于是可以通过下面的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()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

信号的四种频率特性：频谱、频谱密度、能量谱密度、功率谱密度 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) + 3*np.cos(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()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

功率谱