正弦波信号是频率成分最为单一的一种信号，因这种信号的波形是数学上的正弦曲线而得名‌。任何复杂信号，例如光谱信号，音频信号，都可以通过‌傅里叶变换分解为许多频率不同、幅度不等，相位不同的正弦信号的叠加。傅里叶变换主要用在对波形信号进行频域分析。

        正弦曲线可表示为y=Asin(ωx+φ)+k，定义为函数y=Asin(ωx+φ)+k在直角坐标系上的图象，其中sin为正弦符号，x是直角坐标系x轴上的数值，y是在同一直角坐标系上函数对应的y值，k、ω和φ是常数（k、ω、φ∈R且ω≠0）

A——振幅，当物体作轨迹符合正弦曲线的直线往复运动时，其值为行程的1/2。

(ωx+φ)——相位，反映变量y所处的状态。

φ——初相，x=0时的相位；反映在坐标系上则为图像的左右移动。

k——偏距，反映在坐标系上则为图像的上移或下移。

ω——角速度， 控制正弦周期(单位弧度内震动的次数)。

本文使用python语言演示对波形的时频分析，并通过傅立时变换对其进行频域分析。

1、生成正弦波信号进行时域分析

 首先对一段正弦波信号进行采样，正弦波信号需要确定其振幅、角速度和初始相位，以一定的时间间隔对该波形进行采样，生成采样点数据，每组采样点数据为（time，amplitude)，以下python代码显示了以600HZ的采样频率，连续采集1200个点，然后显示波形的振幅随时间变化的时域图：

import numpy as np
import matplotlib.pyplot as plt

# 采样点个数
N = 1200  
# 采样频率 600HZ ，采样周期1.0/600
T = 1.0 / 600.0  # 采样周期  
# 总采样时长 1200* 1.0/600 == 2
# 生成采样时间点序列，在2个时间内，共生成1200个采样时间点
times = np.linspace(0.0, N*T, N, endpoint=False)  

# 第一个子图：时域波形
# 正弦波1：频率为1HZ，振幅为2，初始相位为 0.5*pi 弧度
# 按times的时间采样点序列对正弦波1进行采样，得到的振幅采样点序列
amplitudes1 = 2*np.sin(2*np.pi*times+0.5*np.pi)  
plt.subplot(211)
plt.plot(times,amplitudes1 )
plt.title('Time-Amplitude domain signal')  
plt.xlabel('Time [sec]')  
plt.ylabel('Amplitude')  

# 第二个子图：时域波形
# 正弦波2：频率5HZ，振幅为1，相位为0
# 按times的时间采样点序列对正弦波2进行采样，得到的振幅采样点序列
amplitudes2 = np.sin(2*np.pi*times*5 )  
plt.subplot(212)
plt.plot(times,amplitudes2)
plt.xlabel('Time [sec]')  
plt.ylabel('Amplitude')  

plt.tight_layout()
plt.show()

运行以上程序，结果如下： 

在上面的图中可以看到，通过采样的方式显示的两个波形，其频率、相位、振幅皆不相同。

2、快速傅里叶变换(fft)进行频域分析

什么是傅里叶变换？

        法国科学家傅里叶提出，任何一条周期曲线，无论多么跳跃或不规则，都能表示成一组光滑正弦曲线叠加之和。

        假设有一时间域函数：y = f(x)，根据傅里叶的理论它可以被分解为一系列正弦函数的叠加，他们的振幅A，频率ω或初相位φ不同：

        所以傅里叶变换可以把一个比较复杂的函数转换为多个简单函数的叠加，看问题的角度也从时间域转到了频率域，有些的问题处理起来就会比较简单。

        傅里叶变换的目的是可将时域（即时间域）上的信号转变为频域（即频率域）上的信号。

        将一个复杂的波形信号进行分解，即可看出该信号是由哪些正弦波组成，每个正弦波的频率及该频率上的能量是多少。

        随着域的不同，对同一个事物的了解角度也就随之改变，因此在时域中某些不好处理的地方，在频域就可以较为简单的处理。这就可以大量减少处理信号存储量。

        本文所示代码为python代码。

2.1、由多个正弦波叠加一个复杂波形

        在用python编程时，很多模块都提供了傅里叶变换的相关函数，比如numpy、scipy、mapplotlib等，在使用前要先导入相关模块。

        如下python代码生成一个由两个正弦波叠加组成的复杂波形，并显示其时域分析图：

import numpy as np
import matplotlib.pyplot as plt

B = 600     # 1时间单位的采样点数，即采样频率为600HZ
fs = 1.0/B  # 采样周期
T = 2       # 采样的时间长度，对2个时间单位长的波形进行采样
N = B*T     # 在2时间单位内的采样点个数

# 总采样时长 1200* 1.0/600 == 2
# 生成采样时间点序列，在2个时间内，共生成1200个采样时间点
times = np.linspace(0, T, N)  

# 第一个子图：时域分析图
# 波形为两个正弦波叠加的波形：
#   一个正弦波是：频率为3HZ，振幅为1，初始相位为 0；
#   另一个正弦波是：频率为7HZ，振幅为0.5，初始相位为 0.5*pi ；
# 按times的时间序列对叠加波进行采样，得到的叠加波形的振幅序列
amplitudes1  = np.sin(3 * 2.0*np.pi*times) + 0.5*np.sin(7 * 2.0*np.pi*times+0.5*np.pi)   
plt.plot(times,amplitudes1 )
plt.title('Time-Amplitude domain signal')  
plt.xlabel('Time [sec]')  
plt.ylabel('Amplitude')  

        运行以上代码，可以看到叠加波形的其振幅随时间的变化曲线，可以看到此波形不再是简单的正弦波。 

 

        接着，我们对此波形进行频域分析。也就是看看此波形可以拆解成哪些正弦波的叠加。

2.2、对叠加波形使用傅里叶变换进行频域分析

        在做频域分析的步骤及相关函数有：         

1、通过采样数与采样周期求得傅里叶变换分解所得曲线的频率数组

freqs = np.fft.fftfreq(采样数量, 采样周期)

        在上面的例子中，采样数量是1200，采样周期是1/600，即一个时间单位采600个数据，共采集两个时间单位的数据。

        由此函数得到的是该复合波形拆解出了多少个正弦波，频率数组中的每一项对应一个正弦波的频率，频率数组中的元素作为频域图的X轴坐标。

2、对原波形采样点的振幅序列amplitudes1，经过快速傅里叶变换得到一个复数数组

np.fft.fft(原函数值序列) -> 目标函数值序列(复数)

        其中，在前面的代码中，原函数采样序列即为波形的采样点振幅序列amplitudes1 。

        得到的复数数组中的每个元素是一个复数，该复数对应在第1步得到的相应频率正弦波，复数的模代表的是该频率正弦波的振幅（某频率上的能量），复数的辐角代表该频率正弦波的相位。

3、对复数数组进行求模运算，得到对应频率正弦波的的能量

np.abs(复数数组）

        得到复数数组中每个复数的模组成的数组

        此数组中每个元素的值代表对应频率上能量（振幅），作为频域图的Y轴坐标。

        将第1步得到的频率数组作为X轴，第3步得到的振幅数组作为Y轴坐标，画频域图，代码如下：

import numpy as np
import matplotlib.pyplot as plt

# 采样点个数
N = 1200  
# 采样频率 600HZ ，采样周期1.0/600
T = 1.0 / 600.0  # 采样周期  
# 总采样时长 1200* 1.0/600 == 2
# 生成采样时间点序列，在2个时间内，共生成1200个采样时间点
times = np.linspace(0, N*T, N)  

# 第一个子图：时域分析图
# 波形为两个正弦波叠加的波形：
#   一个正弦波是：频率为3HZ，振幅为1，初始相位为 0；
#   另一个正弦波是：频率为7HZ，振幅为0.5，初始相位为 0.5*pi ；
# 按times的时间序列对叠加波进行采样，得到的叠加波形的振幅序列
amplitudes1  = np.sin(3 * 2.0*np.pi*times) + 0.5*np.sin(7 * 2.0*np.pi*times+0.5*np.pi)   
plt.subplot(311)
plt.plot(times,amplitudes1 )
plt.title('Time-Amplitude domain signal')  
plt.xlabel('Time [sec]')  
plt.ylabel('Amplitude')  

# 第二个子图：频域分析图
sig_freq = np.fft.fftfreq(N,T)           # 得到分解出的频率数组
sig_fft = np.fft.fft(amplitudes1 )       # 得到分解出的各正弦波
sig_fft_abs = np.abs(sig_fft)            # 求模，得到分解出的各正弦波的能量
plt.subplot(312)
plt.plot(sig_freq,sig_fft_abs)
plt.title('Frequency-Amplitude domain signal')  
plt.xlabel('Frequency [Hz]')  
plt.ylabel('Amplitude')  

# 第三个子图：频域分析，只显示频率-10到10之间的图形
plt.subplot(313)
plt.plot(sig_freq,sig_fft_abs)
ax = plt.gca()                  # 获得当前激活的子图对象
ax.set_xlim(-8, 8)            # 为了看的更清晰，限定只显示频率为-10到10之间的图形
plt.title('Frequency[-10,10]-Amplitude domain signal')  
plt.xlabel('Frequency [Hz]')  
plt.ylabel('Amplitude')  

plt.tight_layout()
plt.show()

        运行以上代码，其图形如下：

        其中：
        1）第一张子图是叠加波形的时域图，即振幅随时间的变化曲线
        2）第二张子图是叠加波形的频域图，即叠加波形通过傅立叶变换后得到的各正弦波的频率及频率上的振幅（能量）
        3）第三张子图是将频域图的-8hz到8hz之间的部分的放大图

3、波形幅值的归一化和分贝（dB）

        在上图的频域图中可以看到，一个时间周期内的采样点数是600，其频率为-300到300，在0Hz左右对称，振幅的最大值为600。可以做以下改进：
        1）频率只需要考虑0到300之间的就可以，另一边的频率是对称的。
        2）振幅是一个相对概念，将其归一化到0-1之间。
        改进后的代码如下：

import numpy as np
import matplotlib.pyplot as plt

B = 600     # 1时间单位的采样点数，即采样频率为600HZ
fs = 1.0/B  # 采样周期
T = 2       # 采样的时间长度，对2个时间单位长的波形进行采样
N = B*T     # 在2时间单位内的采样点个数

# 总采样时长 1200* 1.0/600 == 2
# 生成采样时间点序列，在2个时间内，共生成1200个采样时间点
times = np.linspace(0, T, N)  

# 第一个子图：时域分析图
# 波形为两个正弦波叠加的波形：
#   一个正弦波是：频率为3HZ，振幅为1，初始相位为 0；
#   另一个正弦波是：频率为7HZ，振幅为0.5，初始相位为 0.5*pi ；
# 按times的时间序列对叠加波进行采样，得到的叠加波形的振幅序列
amplitudes1  = np.sin(3 * 2.0*np.pi*times) + 0.5*np.sin(7 * 2.0*np.pi*times+0.5*np.pi)   
plt.subplot(311)
plt.plot(times,amplitudes1 )
plt.title('Time-Amplitude domain signal')  
plt.xlabel('Time [sec]')  
plt.ylabel('Amplitude')  

# 第二个子图：频域分析图
sig_freq = np.fft.fftfreq(N,fs)           # 得到分解出的频率数组
sig_freq = sig_freq[:N//2]                # 取一半的频率
sig_fft = np.fft.fft(amplitudes1 )       # 得到分解出的各正弦波
sig_fft_abs = np.abs(sig_fft)            # 求模，得到分解出的各正弦波的能量
sig_fft_abs = sig_fft_abs[:N//2]        # 取一半频率的振幅
sig_fft_abs = sig_fft_abs / B        # 归一化到[0,1]之间
# plt.plot(xf, np.abs(Y[:N//2])/B)  

plt.subplot(312)
plt.plot(sig_freq,sig_fft_abs)
plt.title('Frequency-Amplitude domain signal')  
plt.xlabel('Frequency [Hz]')  
plt.ylabel('Amplitude')  

# 第三个子图：频域分析，只显示频率0到10之间的图形
plt.subplot(313)
plt.plot(sig_freq,sig_fft_abs)
ax = plt.gca()                  # 获得当前激活的子图对象
ax.set_xlim(0, 8)            # 为了看的更清晰，限定只显示频率为-10到10之间的图形
plt.title('Frequency[0,8]-Amplitude domain signal')  
plt.xlabel('Frequency [Hz]')  
plt.ylabel('Amplitude')  

plt.tight_layout()
plt.show()

        运行以上程序，结果为：

        从上图中可以清晰的看到，将叠加波形通过傅里叶变换，拆解成两个正弦波，一个正弦波的频率为3Hz，振幅为1，另一个正弦波的频率为7Hz，振幅为0.5，拆解是成功的。

        上面的程序中，是通过对得到的频率数组和频率所对应的能量数组取一半，来只考虑正数部分，也可以通过下面的代码来只显示正数部分的频率：

# 也可以用以下办法筛选掉频率为负的部分
sig_freq = np.fft.fftfreq(N,fs)           # 得到分解出的频率数组
sig_fft = np.fft.fft(amplitudes1 )       # 得到分解出的各正弦波
sig_fft_abs = np.abs(sig_fft)            # 求模，得到分解出的各正弦波的能量
mask = np.where(sig_freq >= 0)
sig_freq = sig_freq[mask]
sig_fft_abs = sig_fft_abs[mask]/B
plt.subplot(322)
plt.plot(sig_freq,sig_fft_abs)
plt.title('Frequency-Amplitude domain signal')  
plt.xlabel('Frequency [Hz]')  
plt.ylabel('Amplitude')  

        因为波形数据其实是一个相对值，在程序中，对波形的采样频率是600Hz，所以在做傅里叶变换时得到各频率能量的最大值是600，此处的600是一个相对值，有时会将其归一化到-1.0到1.0之间，将数组中的各项除以600就可以了，因为频率正负对称，如果只考虑正的部分，则其最大值为300，归一化到-1.0到1.0之间，只需除以300就可以了。

        在一个复杂的波形中，比如音频，其采样率非常大（44100），通过傅里叶变换分出的波形分量很多，每个波形分量的振幅相差很大，为了更直观的显示波形的分量，对波形分量级差比较大的复杂波形，也采用分贝（dB）来表示其振幅。

        分贝是国家选定的非国际单位制单位, 是我国法定计量单位中的级差单位，表示为dB，其定义为：“两个同类功率量或可与功率类比的量之比值的常用对数乘以10等于1时的级差” 。 

        比如在处理音频处理中，当计算的对象是功率时，音频的分贝定义如下：

        当计算的对象是振幅时，音频的分贝定义为：

        所以在处理音频时，经常对幅值运用以上公式显示其分贝数。如下程序所示：

import numpy as np
import matplotlib.pyplot as plt

B = 600     # 1时间单位的采样点数，即采样频率为600HZ
fs = 1.0/B  # 采样周期
T = 2       # 采样的时间长度，对2个时间单位长的波形进行采样
N = B*T     # 在2时间单位内的采样点个数

# 总采样时长 1200* 1.0/600 == 2
# 生成采样时间点序列，在2个时间内，共生成1200个采样时间点
times = np.linspace(0, T, N)  

# 第一个子图：时域分析图
# 波形为两个正弦波叠加的波形：
#   一个正弦波是：频率为3HZ，振幅为1，初始相位为 0；
#   另一个正弦波是：频率为7HZ，振幅为0.5，初始相位为 0.5*pi ；
# 按times的时间序列对叠加波进行采样，得到的叠加波形的振幅序列
amplitudes1  = np.sin(3 * 2.0*np.pi*times) + 0.5*np.sin(7 * 2.0*np.pi*times+0.5*np.pi)   
plt.subplot(311)
plt.plot(times,amplitudes1 )
plt.title('Time-Amplitude domain signal')  
plt.xlabel('Time [sec]')  
plt.ylabel('Amplitude')  

# 第二个子图：频域分析图，用分贝表示
sig_freq = np.fft.fftfreq(N,fs)           # 得到分解出的频率数组
sig_fft = np.fft.fft(amplitudes1 )       # 得到分解出的各正弦波
sig_fft_abs = np.abs(sig_fft)            # 求模，得到分解出的各正弦波的能量

sig_freq = sig_freq[:N//2]
sig_fft_abs = sig_fft_abs[:N//2]/B
# print(sig_fft_abs[sig_fft_abs.argmax()])    # argmax()求最大分量的位置

plt.subplot(312)
plt.plot(sig_freq, sig_fft_abs, label="real")
ax = plt.gca()                  # 获得当前激活的子图对象
plt.title('Frequency-Amplitude ')  
plt.xlabel('Frequency [Hz]')  
plt.ylabel('Amplitude')  

# 第三个子图：频域分析图

plt.subplot(313)
plt.plot(sig_freq, 20*np.log10(abs(sig_fft_abs)), label="real")
ax = plt.gca()                  # 获得当前激活的子图对象
#ax.set_ylim(-10, 0)
plt.title('Frequency-Amplitude(dB)')  
plt.xlabel('Frequency [Hz]')  
plt.ylabel('Amplitude(dB)')  

plt.tight_layout()
plt.show()

      

4、滤波和快速傅里叶逆变换ifft 

        逆向快速傅里叶变换(IFFT)的计算原理是将频域（注意频域是复数）数据进行取共轭复数（虚部取反），然后再进行FFT变换，这样便将频域信号转换到时域。因为FFT变换的结果是复数，所以从频域进行FFT变换过来的结果也是复数，而此时只需取复数的实部，便是原时域信号。

np.fft.ifft(频率复数序列) -> 叠加波形采样点幅值序列

        经过傅里叶逆变换后，得到叠加波形的采样点幅值序列。

        这样，就可以通过筛选频率复数数组，将该数组中的某些频率上满足一定条件的对应的值设为0，就可以过滤掉某些频率的正弦波分量，进而达到滤波的目的，然后再通过傅里叶逆变换，就可以得到剩余波形的采样点数据。

        在第2节程序的基础上再增加三个子图，以显示不同过滤条件下的波形图。程序如下所示：

import numpy as np
import matplotlib.pyplot as plt

B = 600     # 1时间单位的采样点数，即采样频率为600HZ
fs = 1.0/B  # 采样周期
T = 2       # 采样的时间长度，对2个时间单位长的波形进行采样
N = B*T     # 在2时间单位内的采样点个数

# 总采样时长 1200* 1.0/600 == 2
# 生成采样时间点序列，在2个时间内，共生成1200个采样时间点
times = np.linspace(0, T, N)  

# 第一个子图：时域分析图
# 波形为两个正弦波叠加的波形：
#   一个正弦波是：频率为3HZ，振幅为1，初始相位为 0；
#   另一个正弦波是：频率为7HZ，振幅为0.5，初始相位为 0.5*pi ；
# 按times的时间序列对叠加波进行采样，得到的叠加波形的振幅序列
amplitudes1  = np.sin(3 * 2.0*np.pi*times) + 0.5*np.sin(7 * 2.0*np.pi*times+0.5*np.pi)   
plt.subplot(321)
plt.plot(times,amplitudes1 )
plt.title('Time-Amplitude domain signal')  
plt.xlabel('Time [sec]')  
plt.ylabel('Amplitude')  

# 第二个子图：频域分析图
sig_freq = np.fft.fftfreq(N,fs)           # 得到分解出的频率数组
sig_freq = sig_freq[:N//2]                # 取一半的频率
sig_fft = np.fft.fft(amplitudes1 )       # 得到分解出的各正弦波
sig_fft_abs = np.abs(sig_fft)            # 求模，得到分解出的各正弦波的能量
sig_fft_abs = sig_fft_abs[:N//2]        # 取一半频率的振幅
# plt.plot(xf, 2.0/N * np.abs(Y[:N//2]))  

plt.subplot(322)
plt.plot(sig_freq,sig_fft_abs)
plt.title('Frequency-Amplitude domain signal')  
plt.xlabel('Frequency [Hz]')  
plt.ylabel('Amplitude')  

# 第三个子图：频域分析，只显示频率0到10之间的图形
plt.subplot(323)
plt.plot(sig_freq,sig_fft_abs)
ax = plt.gca()                  # 获得当前激活的子图对象
ax.set_xlim(0, 8)            # 为了看的更清晰，限定只显示频率为-10到10之间的图形
ax.grid()
plt.title('Frequency[-8,8]-Amplitude domain signal')  
plt.xlabel('Frequency [Hz]')  
plt.ylabel('Amplitude')  

# 第四个子图：根据频率作滤波，分解出频率为3的正弦波，利用傅里叶逆逆变换得到其波形采样数组
sig_freq = np.fft.fftfreq(N,fs)           # 得到分解出的频率数组
sig_fft = np.fft.fft(amplitudes1 )       # 得到分解出的各正弦波

sig_fft_filtered = sig_fft * (abs(sig_freq) == 3)   # 筛选频率数组中频率为3Hz的分量
sig_amplitude_filtered = np.fft.ifft(sig_fft_filtered)    # 利用傅里叶逆变换得到频率为3Hz的正弦波采样点数组
plt.subplot(324)
plt.plot(times,sig_amplitude_filtered.real )

plt.title('filtered-signal(3Hz)')  

# 第五个子图：根据频率作滤波，分解出频率为7的正弦波，利用傅里叶逆逆变换得到其波形采样数组
sig_freq = np.fft.fftfreq(N,fs)           # 得到分解出的频率数组
sig_fft = np.fft.fft(amplitudes1 )        # 得到分解出的各正弦波

sig_fft_filtered = sig_fft * (abs(sig_freq) == 7)   # 筛选频率数组中频率为7Hz的分量
sig_amplitude_filtered = np.fft.ifft(sig_fft_filtered)    # 利用傅里叶逆变换得到频率为7Hz的正弦波采样点数组
plt.subplot(325)
plt.plot(times,sig_amplitude_filtered.real )
plt.title('filtered-signal(7Hz)')  

# 第六个子图：根据能量做滤波，分解出频率为7的正弦波，利用傅里叶逆逆变换得到其波形采样数组
sig_freq = np.fft.fftfreq(N,fs)           # 得到分解出的频率数组
sig_fft = np.fft.fft(amplitudes1 )        # 得到分解出的各正弦波
sig_fft_abs = np.abs(sig_fft)            # 求模，得到分解出的各正弦波的能量


sig_fft_filtered = sig_fft * (sig_fft_abs<500)   # 筛选频率数组中频率为7Hz的分量
sig_amplitude_filtered = np.fft.ifft(sig_fft_filtered)    # 利用傅里叶逆变换得到频率为7Hz的正弦波采样点数组

plt.subplot(326)
plt.plot(times,sig_amplitude_filtered.real )
plt.title('filtered-signal(7Hz)') 

plt.tight_layout()
plt.show()

运行上面的程序，其效果如下：

在上面的程序中，子图4是通过代码：

sig_fft_filtered = sig_fft * (abs(sig_freq) == 3)   # 筛选频率数组中频率为3Hz的分量
sig_amplitude_filtered = np.fft.ifft(sig_fft_filtered)    # 利用傅里叶逆变换得到频率为3Hz的正弦波采样点数组

其中：第一行代码只取sig_fft中与频率为3Hz对应的元素的值，与其它频率对应的值都设为0，然后在第二行代码中运行傅里叶逆变换将3Hz的频率波形换算成波形幅值的采样数据。

        子图5，只过滤剩下了频率为7Hz的波形。

        子图6，是通过每种功能的能量大小来进行过率的，代码如下：

sig_fft_abs = np.abs(sig_fft)            # 求模，得到分解出的各正弦波的能量

sig_fft_filtered = sig_fft * (sig_fft_abs<500)   # 筛选频率数组中频率为7Hz的分量
sig_amplitude_filtered = np.fft.ifft(sig_fft_filtered)    # 利用傅里叶逆变换得到频率为7Hz的正弦波采样点数组

其中，第一行代码求出sig_fft中每个频率上的能量（幅值），第二行代码求取sig_fft数组中，对应能量数组sig_fft_abs中少于500的元素所在位置上的值，其它值设为0，因为频率为7Hz的正弦波分量的能量小于500，所以sig_fft中就只剩下7Hz正弦波分量，再通过2傅里叶逆变量，求取此分量的采样点幅值数组。