1、使用matlab进行傅里叶变换
1、傅里叶级数
本人本科不是电气自动化这类专业的，没学过信号与系统，这方面就不太懂，考研的时候只复习了傅里叶级数相关的知识，傅里叶级数在考研中的应用跟泰勒级数的作用基本一样，但是出现的概论远没有泰勒级数出现的概率大，所以其实也就是有一个基本概念而已，这里把这个基本概念再看一下吧：
以下内容均为大一高等数学内容，基本就是抄课本
泰勒级数：
f
(
x
)
=
f
(
a
)
+
f
′
(
a
)
1
!
(
x
−
a
)
+
f
′
′
(
a
)
2
!
(
x
−
a
)
2
+
f
′
′
′
(
a
)
3
!
(
x
−
a
)
3
+
⋯
f(x)=f(a)+\frac{f^{\prime}(a)}{1 !}(x-a)+\frac{f^{\prime \prime}(a)}{2 !}(x-a)^{2}+\frac{f^{\prime \prime \prime}(a)}{3 !}(x-a)^{3}+\cdots
f(x)=f(a)+1!f′(a)(x−a)+2!f′′(a)(x−a)2+3!f′′′(a)(x−a)3+⋯
傅里叶级数：
f
(
x
)
=
a
0
2
+
∑
n
=
1
∞
(
a
n
cos
2
π
n
x
T
+
b
n
sin
2
π
n
x
T
)
f(x)=\frac{a_{0}}{2}+\sum_{n=1}^{\infty}\left(a_{n} \cos \frac{2 \pi n x}{T}+b_{n} \sin \frac{2 \pi n x}{T}\right)
f(x)=2a0+n=1∑∞(ancosT2πnx+bnsinT2πnx)
其中
a
n
=
2
T
∫
x
0
x
0
+
T
f
(
x
)
cos
2
π
n
x
T
d
x
b
n
=
2
T
∫
x
0
x
0
+
T
f
(
x
)
sin
2
π
n
x
T
d
x
\begin{aligned} &a_{n}=\frac{2}{T} \int_{x_{0}}^{x_{0}+T} f(x) \cos \frac{2 \pi n x}{T} d x \\ &b_{n}=\frac{2}{T} \int_{x_{0}}^{x_{0}+T} f(x) \sin \frac{2 \pi n x}{T} d x \end{aligned}
an=T2∫x0x0+Tf(x)cosT2πnxdxbn=T2∫x0x0+Tf(x)sinT2πnxdx
简单来说，泰勒级数是将函数用幂函数的方式进行展开，然后傅里叶级数是将函数用三角函数的方式来进行展开。
2、傅里叶变换
根据对傅里叶变换的了解，我们可以理解为将任何周期哈数分解成一堆正弦和余弦函数，如果出现非周期函数怎么办呢，这就是傅里叶变换的问题了，结合欧拉公式，对傅里叶变换的公式进行进一步简化
e
i
φ
=
cos
φ
+
i
sin
φ
e^{i \varphi}=\cos \varphi+i \sin \varphi
eiφ=cosφ+isinφ
最终的结果如下
f
(
x
)
=
∑
n
=
0
∞
a
n
cos
n
x
+
b
n
sin
n
x
=
∑
−
∞
∞
c
n
e
i
n
x
(
n
=
1
,
2
,
⋯
)
f(x)=\sum_{n=0}^{\infty} a_{n} \cos n x+b_{n} \sin n x=\sum_{-\infty}^{\infty} c_{n} e^{i n x} \quad(n=1,2, \cdots)
f(x)=n=0∑∞ancosnx+bnsinnx=−∞∑∞cneinx(n=1,2,⋯)
其中
c
n
=
1
T
∫
x
0
x
0
+
T
f
(
x
)
e
−
i
2
n
π
x
T
d
x
c_{n}=\frac{1}{T} \int_{x_{0}}^{x_{0}+T} f(x) e^{-i \frac{2 n \pi x}{T}} d x
cn=T1∫x0x0+Tf(x)e−iT2nπxdx
将 w 0 = 2 π T w_0=\frac{2\pi}{T} w0=T2π放到公式中进行替换
{
f
(
x
)
=
∑
n
=
−
∞
∞
c
n
e
i
ω
0
n
x
c
n
=
1
T
∫
x
0
x
0
+
T
f
(
x
)
e
−
i
ω
0
n
x
d
x
\left\{\begin{array}{l} f(x)=\sum_{n=-\infty}^{\infty} c_{n} e^{i \omega_{0} n x} \\ c_{n}=\frac{1}{T} \int_{x_{0}}^{x_{0}+T} f(x) e^{-i \omega_{0} n x} d x \end{array}\right.
{f(x)=∑n=−∞∞cneiω0nxcn=T1∫x0x0+Tf(x)e−iω0nxdx
简单理解下，当T趋于无穷大的时候，
w
0
w_0
w0趋于0，这个时候令
w
=
n
w
0
w=nw_0
w=nw0，就可以得到：
F
(
ω
)
=
1
T
∫
−
∞
∞
f
(
x
)
e
−
i
ω
x
d
x
F(\omega)=\frac{1}{T} \int_{-\infty}^{\infty} f(x) e^{-i \omega x} d x
F(ω)=T1∫−∞∞f(x)e−iωxdx
同时可以得到
f
(
x
)
=
1
2
π
∫
−
∞
∞
F
(
ω
)
e
i
ω
x
d
ω
f(x)=\frac{1}{2 \pi} \int_{-\infty}^{\infty} F(\omega) e^{i \omega x} d \omega
f(x)=2π1∫−∞∞F(ω)eiωxdω
这就分别得到了频率的函数以及由频率函数的进行积分得到的原来的函数，就对应了求频域还有得到频域后进行的逆变换。
该处参考：
https://zhuanlan.zhihu.com/p/104079068
3、尝试使用matlab进行傅里叶变换
那就直接来一段吧，生成一个叠加的正弦波，然后直接调用fft看看啥情况
这里绘出的图像如下所示，也不太对啊，乱七八糟的，这是因为求出来的都是含虚数的，绘图的时候会自动把虚部省略掉，这样就成了这样的图像。
这里我们还需要对生成的信号进行取模和对半，因为复数部分的频率是不需要的，还有就是超过采样频率的另一半也是不需要的，如下图所示
这样变换之后再次绘制结果，就是下面的样子了
完整代码如下：
Fs = 100; % 采样频率，即1s采多少个点
t = (0:1/Fs:10-1/Fs)'; % 1000个采样点
% 生成原始信号
w1 = 10;w2 = 5;
x = 5*sin(2*pi*w1*t)+10*sin(2*pi*w2*t); % 生成图像
% 绘制原始信号
figure(1)
subplot(2,1,1);
plot(x)
title('原始信号','FontSize',25)
xlabel('采样点','FontSize',25)
ylabel('幅值','FontSize',25)
% 进行傅里叶变换
L = length(x); % 总共的数据点数目
y = fft(x);
P2 = abs(y/L); % 取一半就行
P1 = P2(1:L/2+1); % 取交流部分
P1(2:end-1) = 2*P1(2:end-1); % 交流部分模值x2
f = Fs*(0:(L/2))/L;
% 绘制傅里叶变换结果
subplot(2,1,2);
plot(f,P1)
title('幅频图','FontSize',25)
xlabel('频率','FontSize',25)
ylabel('频域幅值','FontSize',25)
如果需要访问其他信号，只要再这里注释掉自己的生成信号修改为其他信号即可
2、使用python进行傅里叶变换
python中实现fft变换还是经典的调库操作，不过就是有两种不同的库，分别是numpy还有scipy这两个库，numpy中有一个fft的库，scipy中也有一个fftpack的库，各自都有fft函数，两者的用法基本是一致的。
numpy这个库比较熟悉了，平时用的也比较多，然后就是scipy这个其实就是基于numpy基础上来的这个库，因此他们的用法才会大同小异也是理所当然的，下面就用python的方法来实现傅里叶变换。
为了方便后面的绘图，在这里加上下面的代码，不然绘图显示中文会异常
mpl.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文
mpl.rcParams['axes.unicode_minus'] = False # 显示负号
这里还是生成一段信号，信号和之前在matlab中的一样
下面是傅里叶变换取模之后在取一半进行显示的结果
运行可以看到结果如下
如果使用采集到的信号的话，将原来的读取函数注释掉，换成我们采集的信号即可
下面尝试对真实情况下采集的信号进行测试
实际情况对比一般两个就够了
源码如下：
import numpy as np
from scipy.fftpack import fft
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文
mpl.rcParams['axes.unicode_minus'] = False # 显示负号
# t = np.linspace(0,1000,1000)
# w1 = 10
# w2 = 5
# y = 5*np.sin(2*np.pi*w1*t)+10*np.sin(2*np.pi*w2*t)
# N = len(t)
y = []
f = open('3-3-1.txt')
for line in f:
y.append(float(line))
N = len(y)
fft_y = fft(y)
abs_y = np.abs(fft_y) # 复数的模
angle_y = np.angle(fft_y) # 复数的角度
normalization_y = 2*abs_y/N # 归一化
normalization_half_y = normalization_y[range(int(N/2))] # 前面的一半归一化
x=np.arange(N) # 时间周期
half_x = x[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.subplot(211)
plt.plot(x, y)
plt.title('原始波形')
plt.subplot(212)
plt.plot(half_x, normalization_half_y, 'blue')
plt.title('单边振幅谱(归一化)')
plt.show()
3、逆变换
逆变换只要再原来的代码上加上ifft的函数即可，如下所示
再进行绘图就可以看到效果了
4、调用matlab工具箱
工具箱如下所示，点击即可进入
打开可以看到下面页面
这里已经准备好了两段数据，数据可以自己生成，也可以用采集的数据
将数据拖进去
设置采样率
可以勾选查看波形
点击查看傅里叶分析结果
查看短时傅里叶分析结果
平移器可以滑动查看波形
还可以进行信号的截取等操作，该处参考使用信号分析器
扫一扫
7674
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
