基于Python（Numpy）的周期信号傅里叶变换

news/2026/1/21 22:12:47/文章来源:https://www.cnblogs.com/sg-lin/p/19439642

引言

傅里叶变换的作用是将时序信号转变为频率域，而在傅里叶变换中，快速傅里叶变换（Fast Fourier Transform，FFT）是一种高效率的算法。

已知连续信号傅里叶变换公式：

$$X(j\omega)= \int_{-\infty}^{+\infty} {x}(t) e^{-jwt}dt$$

实际中一个信号需要通过等时间间隔上的值或样本来表示，在一定的条件下，这些离散的样本值可恢复出连续时间信号，该性质为连续时间信号的采样定理。

由采样定理可知，一个连续时间信号可通过采样的变成离散时间信号，然后在离散时间系统中处理后，在变换回连续时间信号。因此，在编程处理上，使用灵活的离散时间信号作为傅里叶变换算法基础。离散信号（DFT）的表达式，如下：

$$X(e^{j\omega})=\sum_{n=-\infty}^{+\infty}{x}[n] e^{-jwn} X(k)=\sum_{n=-\infty}^{+\infty}{x}[n]e^{\frac{-jk2\pi n}{N} } -------w=\frac{k2\pi}{N}(n=0,±1,±2...)$$

当频率索引k的取正值时，DFT的表达式可为：

$$X(k)=\sum_{n=0}^{N-1}{x}[n]e^{-jk\frac{2\pi }{N}n } -------(k=0,1,2...,N-1)$$

已知一个采样函数的的频率为f_s，采样的数据点为N，而时域中相邻两个数据点的时间间隔为T_s=1/f_s，则N个数据点的时间长度为NT_s，同时该时间周期长度为离散信号傅里叶变换下的扩展周期。

$$N个离散数据的表示式：\begin{cases}x[nT_{s}],n=0,1,2...N-1 \\0,other\end{cases}$$

可知，N个离散数据函数式代入DFT表达式为：

$$ X(k)=\sum_{n=0}^{N-1}{x}[nT_{s}]e^{\frac{-jk2\pi nT_{s}}{NT_{s}}}$$

则实际频率f_k的表达式为：f_k=(k/NT_s)，f_k同时也表示第k个频率点对应的频率，在Python代码中如下:

import numpy as np
from matplotlib import pyplot as plt
from scipy.fft import fft,fftfreq
#采样点数量
SAMPLE_N = 200
#时域信号的长度
DURATION = 1
def generate_cos(a,k,w):x=np.linspace(0,DURATION,SAMPLE_N,endpoint=False)y=a*np.cos(k*np.pi*5*x)#y=a*np.cos(2*np.pi*k*x)-a*np.cos(2*np.pi*k*x)*w+wreturn x,y

Python中的快速傅里叶变换

numpy实现快速傅里叶变换

某个激励函数为：

$$ f(t)=2+4\cos(2\pi 5t)+8\cos(2\pi 10t) 式中的2\pi，由式可知：\left\{\begin{matrix}T= \frac{2\pi}{w}\\f = \frac{1}{T}\end{matrix}\right. ，则w=2\pi f$$

python代码如下：

FFT_Python01

 import numpy as np
from matplotlib import pyplot as plt
from scipy.fft import fft,fftfreq
#采样点
SAMPLE_N = 200
#时域信号的长度
DURATION = 1def generate_cos(a,k,w):x=np.linspace(0,DURATION,SAMPLE_N,endpoint=False)# y=a*np.cos(k*np.pi*5*x)#a*np.cos(2*np.pi*k*x)*w+w 目的为生成直流分量y=a*np.cos(2*np.pi*k*x)-a*np.cos(2*np.pi*k*x)*w+wreturn x,y
x1,y1 = generate_cos(0,0,2) # 直流分量 2
x2,y2 = generate_cos(4,5,0) # 4cos(2π5t)
x3,y3 = generate_cos(8,10,0)# 4cos(2π10t)plt.subplot(2,3,1)
plt.xlabel("DC component")
plt.plot(x1,y1)plt.subplot(2,3,2)
plt.xlabel("4cos(2π5t)")
plt.plot(x2,y2)plt.subplot(2,3,3)
plt.xlabel("4cos(2π10t)")
plt.plot(x3,y3)#f(t)=2+4cos(2π5t)+4cos(2π10t)
y = y1+y2+y3
plt.subplot(2,3,4)
plt.xlabel("2+4cos(2π5t)+4cos(2π10t)")
plt.plot(x3,y)#重点：np.fft.fft()为numpy中的快速傅里叶变换
y_fft = np.fft.fft(y1+y2+y3)
#获得信号的幅频
amplitude_spectrum = np.abs(y_fft)
#获得信号的相频
phase_spectrum = np.angle(y_fft)# N=Fs采样点的数量
Fs = x1.shape[-1]
#频率轴，上式中DFT变换中第k个频率点的频率数组，x1[1]-x1[0]为采样点之间的时间间隔（频率分辨率）
frequencies = np.fft.fftfreq(Fs,(x1[1]-x1[0]))plt.subplot(2,3,5)
plt.plot(frequencies,amplitude_spectrum)
plt.xlabel("Spectrogram")#重点：numpy中的傅里叶逆变换
recovered_signal = np.fft.ifft(y_fft)
plt.subplot(2,3,6)
plt.plot(x1,recovered_signal)
plt.xlabel("ifft")
plt.show()

Numpy快速傅里叶变换的结果如下：

快速傅里叶变化（FFT）中的问题

从上图（Spectrogram）的频谱图结果知，在频率轴上出现了负频率分量（-5Hz）与正频率分量（5Hz）相同的幅值，同时频域中幅值的大小与时域中并不一样，例如5Hz下，时域中幅度值为：4，频域中幅度值为：400 。

频域中出现负频率值和时域中幅度值不匹配情况，可从离散傅里叶（DFT）表达式获得结果，DFT表达式如下：

$$X(e^{j\omega})=\sum_{n=-\infty}^{+\infty}{x}[n] e^{-jwn} X(k)=\sum_{n=-\infty}^{+\infty}{x}[n]e^{\frac{-jk2\pi n}{N} } -------w=\frac{k2\pi}{N}(k=0,±1,±2...)$$

在DFT X(k)表达式中，k的取值(-∞，+∞ )区间，即离散的时域信号经离散傅里叶变换后，其频谱图（频率轴上的取值）可向负轴和正轴扩展。

傅里叶变换性质（共轭和共轭对称性）：

同时由傅里叶变换的共轭和共轭对称性可知，DFT中若

$$ \left\{\begin{matrix}x[n]\Leftrightarrow X(e^{jw})\\ x^{*}[n] \Leftrightarrow X^{*}(e^{jw}) \end{matrix}\right. *表示共轭，当x[n]为实值时，则变换为共轭对称（共轭对称，即a^{*}_{k}=a_{-k}）：$$

$$ X(e^{-jw}) = X^{*}(e^{jw}) x[n]为实值$$

推导过程：

$$X(e^{j\omega})=\sum_{n=-\infty}^{+\infty}{x}[n] e^{-jwn} $$

$$X^{*}(e^{j\omega})=\sum_{n=-\infty}^{+\infty}{x}^{*}[n] e^{jwn} $$

上式中：以-ω代替ω可得

$$X^{*}(e^{-j\omega})=\sum_{n=-\infty}^{+\infty}{x}^{*}[n] e^{-jwn}，若x[n]为实值时，x[n]=x^{*}[n]，X(e^{-j\omega})=X^{*}(e^{j\omega})$$

若将x[n]为实值信号时，

$$ X(j\omega)用笛卡尔坐标系表示为： X(j\omega)=Re\{X(j\omega)\}+jIm\{X(j\omega)\} $$

同样，共轭：

$$X^{*}(j\omega)=Re\{X(-j\omega)\}-jIm\{X(-j\omega)\} $$

由上可知：

$$ \left\{\begin{matrix} Re\{X(j\omega)\} = Re\{X(-j\omega)\} \\ Im\{X(j\omega)\} = -Im\{X(-j\omega)\} \end{matrix}\right. $$

通过上式推导过程可知，在实值信号的傅里叶变换中实部是频率的偶函数，而虚部是频率的奇函数，因此频谱图中出现负频率分量与正频率分量对称情况，最后研究实值信号傅里叶变换时，只需要给出正频率值（w>0），而负频率值可利用共轭对称性导出。

故在频域中的实值信号变换时，频谱图中只需绘制出正频率轴的信息，设置正频率轴的python代码如下：

Fs = x1.shape[-1]
#生成频率轴
'''n 是总采样数量点，d 相邻采样点的时间间隔如果n时偶数时 f = [0, 1, ...,   n/2-1,     -n/2, ..., -1] / (d*n)   如果n时奇数时 f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (d*n)  Fs: 总采样数量点, x1[1]-x1[0]: 相邻采样点的时间间隔
'''
frequencies = np.fft.fftfreq(Fs,x1[1]-x1[0])
#向下取整，由numpy中频率轴的生成的数组可知数组前（length of frequency axis）/2为正频率数组
N = Fs//2
plt.plot(frequencies[:N],amplitude_spectrum[:N])

频谱图（正频率轴）

频域中幅度值与时域中不匹配情况这是由离散傅里叶变换的本质确定的，X(jω) 函数称为频谱密度函数，因此此时频域中的幅度值不是真实的。

在离散傅里叶逆变换中：

$$ x[n]=\frac{1}{2\pi } \int_{2\pi}^{}X(e^{j\omega})d\omega，其中\frac{1}{2\pi }X(e^{j\omega})d\omega=X(e^{j\omega})df 为各分量的振幅，是无穷小量，所以 X(j\omega)不能表示真实的幅度，改用密度函数表示。$$

已知一个离散周期信号，在区间[N1,N2]为x[n]，而在该区间之外x[n]=0，则该周期信号的傅里叶级数系数a_k（幅度值）和傅里叶变换如下式：

$$ a_{k}=\frac{1}{N}\sum_{n=-N1}^{+N2}x[n]e^{-jk\frac{2\pi}{N}n }=\frac{1}{N}\sum_{n=-\infty}^{+\infty}x[n]e^{-jk\frac{2\pi}{N}n } 傅里叶级数系数$$

$$ X(e^{jk\frac{2\pi}{N} })=\sum_{n=-\infty }^{\infty}x[n]e^{-jk(\frac{2\pi}{N} )n} 傅里叶变换$$

由上可知a_k与X(e^jω)的各值关系如下，与傅里叶级数的系数之间已经存在1/N的差距。

$$ a_{k}=\frac{1}{N}X(e^{jk\omega_{0}}) $$

帕斯瓦尔定理：时域信号的总能量与频域中的信号相等

$$ \sum_{n=0 }^{N-1} |x[n]|^{2} =\frac{1}{N} \sum_{n=0 }^{N-1} |X[k]|^{2}，若x[n]为实值信号，则\sum_{n=0 }^{N-1} x^{2}[n] =\frac{1}{N} \sum_{n=0 }^{N-1} |X[k]|^{2} $$

上式可知，|X[k]|不能作为幅值，其物理意义不明确，要归一化处理。

根据DFT的共轭对称性可知：

$$ X(e^{-j\omega})=X^{*}(e^{j\omega})\Longleftrightarrow X[k]=X^{*}[N-k] $$

X[k]=X^*[N-k] 的证明：

$$ \left\{\begin{matrix}X(k)=\sum_{n=-\infty}^{+\infty}{x}[n]e^{\frac{-jk2\pi n}{N} } \\X(N-k)=\sum_{n=-\infty}^{+\infty}{x}[n]e^{\frac{jk2\pi n}{N} }e^{-j2\pi n}\end{matrix}\right.，其中e^{-j2\pi n}=1(欧拉公式，\cos( 2\pi n) - j\sin( 2\pi n)=1)，故X(N-k)=\sum_{n=-\infty}^{+\infty}{x}[n]e^{\frac{jk2\pi n}{N} } $$

DFT与连续傅里叶系数的关系：

假设有一个单频实值信号：

$$ x[n]=A_{k_{0}}\cos(k_{0} \frac{2\pi}{N}n) $$

由欧拉公式，可得：

$$ \cos x = \frac{ e^{jx}+e^{-jx}}{2} 欧拉公式 $$

代入DFT中公式中：

$$ X[k]=\sum_{n=0}^{N-1} (\frac{A_{k_{0}}}{2}e^{jk_{0}\frac{2\pi}{N}n}+ \frac{A_{k_{0}}}{2}e^{-jk_{0}\frac{2\pi}{N}n})e^{-jk\frac{2\pi}{N}n } $$

利用离散复指数的正交性可得：

$$ \sum_{n=0}^{N-1} e^{j(k_{0}-k)\frac{2\pi}{N}n}=\left\{\begin{matrix}N,k=k_{0}
\\0,k\ne k_{0}
\end{matrix}\right. $$

所以当k=k₀时，正频率中幅度值为:

$$ X[k_{0}]=\frac{A_{k_{0}}}{2} N\Longrightarrow A_{k_{0}}=\frac{2}{N}X[k_{0}]$$

同时对于实值信号的共轭对称性（X[k]=X^*[N-k]）可知，负频率分量中的幅度值为:

$$ X[N-k_{0}]=\frac{A_{k_{0}}}{2} N\Longrightarrow A_{k_{0}}=\frac{2}{N}X[N-k_{0}]$$

因此离散实值信号在傅里叶变换后，获得真实幅度值时，需要乘2/N。

注意特殊情况：

在零频率分量（直流分量，k=0⇒ω=0）时：

$$ x[0]=A_{0}\Longrightarrow X[0]=\sum_{n=0}^{N-1}A_{0}=NA_{0} $$

所以零频率分量(ω=0)的幅度值为：

$$ A_{0}=\frac{X[0]}{N} $$

因此在零频率分量（ω=0）处的真实幅度值需要乘1/N。

奈奎斯特采样定理:

时域信号在傅里叶变换后，频率轴上的最大频率点值为f_s/2（f_s采样频率），频率点之间的间隔（频率分辨率）为f_s/N，由上面离散实值信号的共轭对称性知，只需研究频域中[0,N/2]频率点。

当n为偶数时，频域中将产生N/2+1个频率点（其中包括奈奎斯特频率点k=N/2和零频率点k=0），此时奈奎斯特频率点没有共轭对称点，由帕斯瓦尔定理中可知，时域中的总能量与频域中总能量相等，因此为平衡能量，奈奎斯特频率点集中了共轭点的能量，该频率点对应的幅度值需要乘1/N。

当n为奇数时，频域中将产生(N+1)/2个频率点，此时奈奎斯特频率点存在共轭对称点，因此该频率点的幅度值需要乘1/N。

故在频域中的真实幅度值需要修正，python代码如下：

#获得频谱密度函数值
amplitudeSpectrum = np.abs(y_fft)
#获得相位值
phaseSpectrum = np.angle(y_fft)
#连接处理后的真实的幅度值
yAmpModified = np.concatenate(((amplitudeSpectrum[0]).reshape(1)/Fs,(amplitudeSpectrum[1:N])*2/Fs,(amplitudeSpectrum[N]).reshape(1)*1/Fs))
plt.plot(frequencies[:N],yAmpModified[:N])
plt.xlabel("Spectrogram")

频谱图（修正结果）

总结

综合上面分析可知，快速傅里叶变换中，频域信号幅度值为非真实的，需要对傅里叶变换后的幅度值进行修正。在奈奎斯特采样定理中，频率轴上最大频率点值为f_s/2，频率点之间的间隔（频率分辨率）为f_s/N，对于离散周期实值信号，根据共轭对称性知（X[k]=X^*[N-k] ），频域中只需研究正频率分量0~N/2（频率轴上的取值范围 [0:(N/2) * (f_s/N)]）。根据帕斯瓦尔定理可知，时域信号的总能量与频域中的信号相等的，频域中需要归一化（即除以N），以平衡总能量，由于频域中负频率分量与正频率分量对称相等，故两个对称的频率点均分能量，表明频域中正频率分量中只包含一半的幅度值，另一半的幅度值在负频率分量中。最后综上分析，快速傅里叶变换中，频谱图中的幅度值经过修正才能获得真实的值，即频率点对应的幅度值需要乘以2/N进行修正，除特殊直流分量和奈奎斯频率点（N为偶数）乘以1/N。

完整代码

Completed Code

import numpy as np
from matplotlib import pyplot as plt
from scipy.fft import fft,fftfreq
#采样点
SAMPLE_N = 200
#时域信号的长度
DURATION = 1def generate_cos(a,k,w):x=np.linspace(0,DURATION,SAMPLE_N,endpoint=False)# y=a*np.cos(k*np.pi*5*x)#a*np.cos(2*np.pi*k*x)*w+w 目的为生成直流分量y=a*np.cos(2*np.pi*k*x)-a*np.cos(2*np.pi*k*x)*w+wreturn x,y
x1,y1 = generate_cos(0,0,2) # 直流分量 2
x2,y2 = generate_cos(4,5,0) # 4cos(2π5t)
x3,y3 = generate_cos(8,10,0)# 4cos(2π10t)plt.subplot(2,3,1)
plt.xlabel("DC component")
plt.plot(x1,y1)plt.subplot(2,3,2)
plt.xlabel("4cos(2π5t)")
plt.plot(x2,y2)plt.subplot(2,3,3)
plt.xlabel("4cos(2π10t)")
plt.plot(x3,y3)#f(t)=2+4cos(2π5t)+4cos(2π10t)
y = y1+y2+y3
plt.subplot(2,3,4)
plt.xlabel("2+4cos(2π5t)+4cos(2π10t)")
plt.plot(x3,y)#重点：np.fft.fft()为numpy中的快速傅里叶变换
y_fft = np.fft.fft(y1+y2+y3)
#获得信号的幅频
amplitudeSpectrum = np.abs(y_fft)
#获得信号的相频
phaseSpectrum = np.angle(y_fft)# N=Fs采样点的数量
Fs = x1.shape[-1]
#获得正频率分量的频率轴
N = Fs//2
#频率轴，上式中DFT变换中第k个频率点的频率数组，x1[1]-x1[0]为采样点之间的时间间隔（频率分辨率）
frequencies = np.fft.fftfreq(Fs,(x1[1]-x1[0]))
#连接处理后的真实的幅度值
yAmpModified = np.concatenate(((amplitudeSpectrum[0]).reshape(1)/Fs,(amplitudeSpectrum[1:N])*2/Fs,(amplitudeSpectrum[N]).reshape(1)*1/Fs))
plt.subplot(2,3,5)
plt.plot(frequencies[:N],yAmpModified[:N])
plt.xlabel("Spectrogram")#重点：numpy中的傅里叶逆变换
recoveredSignal = np.fft.ifft(y_fft)
plt.subplot(2,3,6)
plt.plot(x1,recoveredSignal)
plt.xlabel("ifft")
plt.show()

快速傅里变换后频谱图修正的结果如下：