python fft ifft

条件

任何一个满足狄利克雷条件的函数都可以通过傅里叶基数展开。
numpy和scipy中都有fft变换，且效果都是一样的。

代码

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       #显示负号

#原始数据生成
x=np.linspace(0,2,3000)#设置自变量在0、2上采样3000个点，那么采样率为1500点/x轴单位
y=10+7*np.sin(2*np.pi*30*x) + 5*np.sin(2*np.pi*60*x+10/180*np.pi)+3*np.sin(2*np.pi*100*x)#常量幅值是10，在30Hz，60Hz，100Hz上，对应的幅值未7、5、3对应的相位为0度、10度、0度
#原始数据绘制
plt.subplot(221)
plt.plot(x,y,'-')
plt.title('原始数据')


fft_y=fft(y)#傅里叶变换

yreal =fft_y.real#获取实部
yimag = fft_y.imag#获取虚部，虚部和实部是要计算相角。
ydegree = np.angle(fft_y)*180/np.pi#复数的辐角就是相角，输出的是弧度，这里转为角度

yabs = abs(fft_y)#计算频率对应的幅值
yabs = yabs/len(x)*2#幅值归一化，除了常数项以外，都除以采样点数的一半
yabs[0]/=2#常数项除以采样点数

#计算幅值和相角对应的频率，
#幅值和相角会计算采样点数个数据，这里就是3000个点，第一个是常量剩下2999个点，幅值的2999点关于中间对称，相角的2999点关于中间奇对称，即1号数据=-2999号数据
#所以有用的幅值和相角就是1500个数据，而采样率是1500Hz，那么一个幅值对应一个整数频率
#如果x=np.linspace(0,1,3000)那么，一个幅值对应两个整数频率，所以要想幅值和频率对应上要获得频率坐标轴如下：
#x轴上的一个1表示一个样周期，所以如果x取值1到3取5000个点，那么就经过了两个周期，相当于一个采样周期采集2500个点，即采样率是2500Hz，所以这里要除以范围，才能和计算得到的频率幅值对应
xfft = range(len(x))/(max(x)-min(x))
#绘制相频图
plt.subplot(222)
plt.plot(xfft,ydegree, '-')
plt.title('相频图')
#绘制幅频图
plt.subplot(223)
plt.plot(xfft,yabs, '*')
plt.title('幅值归一化')

#开始分析
#发现相频图数据很乱，因为大部分相角对应的幅值很小，所以不管其相角是0还是80在整体上意义不大，我们只关心大幅值对应的相角，只有这些相角才会影响整体数据走向，所以提取大的幅值对应的相角
for i in range(len(yabs)):
    if yabs[i]>2:print(xfft[i],ydegree[i]+90)
"""
打印数据：
0.0 0.0
30.0 -86.39368585633036
60.0 -72.7941709180644
100.0 -78.00623508482248
1400.0 78.00623508482248
1440.0 72.7941709180644
1470.0 86.39368585633036
"""
#可以看到幅值是找对了，但是相角好像差的很大，但是这些相角加上90度会发现较为接近设置值，因为FFT展开默认展开成cos，刚好和sin相差90度


#开始根据展开反算原始函数，看看是否能对应上

#通过自带函数反变换
fft_y_n=ifft(fft_y)
plt.subplot(221)
plt.plot(x,fft_y_n,'r-',lw="0.5")

#通过幅值和角度复原函数
y_n= []
for i in xfft:
    temp = [ yabs[j]*np.cos(i*2*np.pi*x[j]+ydegree[j]/180*np.pi) for j in range(int(len(yabs)/2))]#这里需要注意两点，第一点是i从xfft中取，第二点是range(int(len(yabs)/2))用一半数据叠加
    y_n.append(sum(temp))
    
fft_y_n=ifft(fft_y)
plt.plot(x,y_n,'g-',lw="0.5")
plt.show()
#三条曲线重合

实例

一段没有横坐标是的数据，（其实数据一定会有横坐标的，最不行可以用index作为横坐标。）如每个小时的温度数据，假如说有50个小时，那么我认为这50个小时对应从x坐标是0-1那么时间和x坐标就有对应关系了，1小时对应0.2，根据上面进行FFT展开，可以获得x和温度的不同频率三角函数，把x用时间替换就是时间关于频率的三角函数。

import numpy as np
from scipy.fftpack import fft,ifft
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl
import sys

mpl.rcParams['font.sans-serif'] = ['SimHei']   #显示中文
mpl.rcParams['axes.unicode_minus']=False       #显示负号

h=np.linspace(0,59,60)#改为np.linspace(0,47,48)
tem =10+7*np.sin(2*np.pi/24*h)#频率是1/24
#定义所有的数据构成一个周期,与上面的tem数据对应，x上一个是1/60小时
x=np.linspace(0,1,60)#改为np.linspace(0,1,48)后清晰可见


fft_y=fft(tem)#傅里叶变换

yreal =fft_y.real#获取实部
yimag = fft_y.imag#获取虚部，虚部和实部是要计算相角。
ydegree = np.angle(fft_y)*180/np.pi#复数的辐角就是相角，输出的是弧度，这里转为角度

yabs = abs(fft_y)#计算频率对应的幅值
yabs = yabs/len(x)*2#幅值归一化，除了常数项以外，都除以采样点数的一半
yabs[0]/=2#常数项除以采样点数

#绘制相频图
plt.subplot(221)
plt.plot(x,ydegree, '-')
plt.title('相频图')
#绘制幅频图
plt.subplot(222)
plt.plot(x,yabs, '*')
plt.title('幅值归一化')

#幅值最大值应该在1/24处，显然没有1/24这个值，那么就在2/60和3/60之间，所以在第2、4个数上大道最大，修改采样点数到频率到 组成的每个函数频率可以整除采样点数分之1即可实现，每个频率上清晰
plt.show()