离散傅里叶变换用于频率分析
在本文中，我们将探讨离散傅里叶变换的基础知识，以及它如何用作信号频率分析的工具。
引言
傅里叶分析通常关注函数的分析和合成。将信号分解为易于分析的组件，并从这些组件进行重构。
在本文中，我们将讨论离散时间信号的傅里叶分析。
离散时间信号仅在特定的时间点定义，并表示为具有连续值范围的实数序列。
信号的傅里叶级数表示
离散时间的复指数信号本质上是周期性的。
\( e^{j k\omega (n+N)} = e ^{j k \omega n} \)
其中 , \(\omega=\frac{2\pi}{N}\)
所有周期为 N 的离散时间复指数信号集为 \(\phi_k[n] = e ^{j k \frac{2\pi}{N} n} , \forall k \in \mathcal{Z}\)
当 k 改变整数倍 N 时，我们得到相同的序列。
因此，只有 N 个 k 值可以定义唯一的离散时间复指数序列。
\(e^{j k \omega n}\) 仅在 k 的 N 个连续值范围内是唯一的，求和仅在此范围内考虑。
(x[n]) 具有周期 N。
离散时间傅里叶级数
我们考虑一个周期为 N 的任意周期信号 (x[n])，并且 (x[n]) 可以表示为：
因此，我们需要确定 (x[n]) 是否可以表示为上述形式，如果可以，则找到使此表示有效的系数 (x[k])。
所有信号处理参考资料中都可以找到相同的理论推导。结果如下：
和
(X[k]) 称为傅里叶级数系数。这些系数也称为信号 (x[n]) 的频谱系数。
信号 (x[n]) 的频谱系数表示称为信号的傅里叶级数表示。
这些系数指定了信号 (x[n]) 分解为 N 个谐波相关的复指数之和。
任何周期性的离散时间信号 (x[n]) 都可以使用傅里叶级数表示来表示，傅里叶级数表示使我们能够将任何周期信号表示为复指数的加权和。
让我们看一些例子来理解傅里叶级数表示能提供什么信息。
```
def FourierSeries(input,N=None):
""" computes the fourier series coefficients of input signal
Parameters
-----------
input : numpy array
input discrete time signal
Returns
--------
out : complex,list
fourier series coefficients
"""
N=len(input);
w=2*cmath.pi/N;
input=input[0:N];
n=numpy.arange(0,N);
r=cexp(-1j*w*n);
output = [complex(0)] * N
for k in range(N):
r=input*cexp(-1j*w*n*k)
output[k]=np.sum(r);
return output;
```
我们考虑一个频率为 (f_{s}=10Hz) 的正弦信号，持续时间为 1 秒，并以 (F_{s}=150Hz) 对波形进行采样，然后观察信号的傅里叶级数系数。
在频率谱中，仅考虑余弦波形的单个周期，即 (Fs/f_{s}=15) 个样本。结果可以在子图 2 中看到。
如果我们获取完整的信号并计算傅里叶级数系数，结果可以在子图 3 中看到。
让我们考虑傅里叶级数的幅度图。
```
def FourierSinusoids(F,w,Fs,synthesis=None):
""" the function generates a discrete time sinusoid,computes
the fourier series coefficients and plots the time and frequency
data
Paraneters
----------
F : numpy array
sinuoidal frequency components
w : numpy array
the weights associated with freuency components
Fs : Integer
sampling frequency
synthesis : int
if 1 ,reconstructs signal and plots the original and reconstructed
signal
"""
if synthesis==None:
synthesis=0;
Ts=1.0/Fs;
xs=numpy.arange(0,1,Ts)
signal=numpy.zeros(np.shape(xs));
for i in range(len(F)):
omega=2*np.pi*F[i];
signal = signal+ w[i]*numpy.cos(omega*xs);
#plot the time domain signal
subplot(2,1,1)
plt.plot(range(0,len(signal)),signal)
xlabel('Time')
ylabel('Amplitude')
title('time doman')
#plt.ylim(-2, 2)
#compute fourier series coefficients
r1=FourierSeries(signal)
a1=cabs(r1)
if synthesis==0:
#plot the freuency domain signal
L=len(a1);
fr=np.arange(0,L);
subplot(2,1,2)
plt.stem(fr,a1,'r') # plotting the spectrum
xlabel('Freq (Hz)')
ylabel('|Y(freq)|')
title('complete signal')
ticks=np.arange(0,L+1,25);
plt.xticks(ticks,ticks);
show()
if synthesis==1:
rsignal=IFourierSeries(r1);
print np.allclose(rsignal, signal)
subplot(2,1,2)
plt.stem(xs,signal)
xlabel('Time')
ylabel('Amplitude')
title('reconstructed signal')
show()
if __name__ == "__main__":
mode =0
if mode==0:
F=[10];
F=np.array(F);
w=numpy.ones(F.shape);
#plot the time domain signal and fourier series component
FourierSinusoids(F,w,150);
```
我们可以看到，在对应于 10Hz 的频率处观察到峰值。傅里叶级数系数的总数等于输入样本的总数。
现在我们考虑两个频率为 10Hz 和 15Hz 的正弦信号的组合。该信号的周期是 10 和 15 的 LCM，即 5Hz。
```
if mode==1:
F=[10,20];
F=np.array(F);
w=numpy.ones(F.shape);
FourierSinusoids(F,w,150);
```
现在我们继续添加频率分量。下面的信号周期为 1 秒。
```
if mode==2:
F=range(1,10);
F=np.array(F);
w=numpy.ones(F.shape);
FourierSinusoids(F,w,150);
```
信号合成
在上一节中，我们看到了将周期信号分解为其频率分量的过程。
现在我们将讨论从傅里叶级数系数合成信号。
下方是原始信号和重构信号的图以及相应的 Python 代码。
```
def IFourierSeries(input):
""" function reconstructs the signal from fourier series coefficients
Parameters
---------
input : cmath list
fourier series coefficients
Returns
-----------
out : numpy arrary
reconstructed signal
"""
N=len(input);
w=2*cmath.pi/N;
k=numpy.arange(0,N);
output = [complex(0)] * N
for n in range(N):
r=input*cexp(-1j*w*n*k);
output[n]=np.mean(r);
print output.__class__
return output;
if mode==3:
F=range(1,10);
F=np.array(F);
w=numpy.ones(F.shape);
FourierSinusoids(F,w,150,1);
```
离散傅里叶变换
现在出现了一个问题，对于非周期信号，我们该怎么办？经过对离散时间傅里叶变换和频率域采样的大量理论分析，我们发现只需假设非周期信号的周期扩展，然后按照上述方法计算傅里叶级数。
周期信号的傅里叶级数系数也是周期性的，周期相同为 N。
如果我们考虑傅里叶级数系数的单个周期 N 个值，我们就得到一个有限长度的序列。
这被称为离散傅里叶变换 (DFT)。
因此，通过计算 DFT，我们得到单个周期的傅里叶级数系数。
选择信号的周期取决于我们。我们考虑一个长度为 150、占空比为 5 的非周期冲激信号。
我们考虑 N=150, 450 并观察结果。
当增加信号的周期时，我们可以看到频率域的分辨率增加。
当 $N -> \infty$ 时，频率域的采样点会越来越密集。如果频率域的采样点无限密集，则可以将其视为连续信号。
这就得到了信号的离散时间傅里叶变换表示。
```
def FourierRect(N):
""" the function generates rectangular aperiodic pulse and computer the DFT coefficients
Parameters
----------
N : period of aperiodic signal
"""
x = np.zeros((1,N))
x[:,0:30]=1;
x=x.flatten();
#compute the DFT coefficients
r1=FourierSeries(x)
#magnitude of DFT coefficients
a1=cabs(r1)
#plot the time domain signal
subplot(2,1,1)
plt.plot(range(0,len(x)),x)
xlabel('Time')
ylabel('Amplitude')
title('time doman')
plt.ylim(-2,2);
#plot the DFT coefficients
L=len(a1);
fr=np.arange(0,L);
subplot(2,1,2)
plt.stem(fr,a1,'r') # plotting the spectrum
xlabel('Freq (Hz)')
ylabel('|Y(freq)|')
title('complete signal')
ticks=np.arange(0,L+1,25);
plt.xticks(ticks,ticks);
show()
if mode==4:
FourierRect(150);
if mode==5:
FourierRect(150*3);
```
相邻采样点之间的距离是 (\omega )，当 (w \rightarrow ;0) 时，采样点无限密集。
傅里叶变换本质上是连续的，不能用于数值计算。
离散傅里叶变换是信号的离散时间傅里叶变换的采样版本，并且是以适合在信号处理单元上进行数值计算的形式。
快速傅里叶变换 (FFT) 是一种计算离散傅里叶变换 (DFT) 及其逆运算的算法。它是计算信号 DFT 的有效方法。
我们将使用 Python 的 FFT 例程，并将其性能与朴素实现进行比较。
使用内置 FFT 例程： elapsed time was 6.8903e-05 seconds
使用朴素代码： elapsed time was 0.0653119 seconds
我们可以看到大约 1000 倍的性能提升。
朴素算法的复杂度为 (o(N^2))，FFT 算法的复杂度为 (O(N log N))。
```
if mode==6:
Fs=150;
F=range(1,10);
F=np.array(F);
w=numpy.ones(F.shape);
Ts=1.0/Fs;
xs=numpy.arange(0,1,Ts)
signal=numpy.zeros(np.shape(xs));
for i in range(len(F)):
omega=2*np.pi*F[i];
signal = signal+ w[i]*numpy.cos(omega*xs);
start_time = time.time()
FourierSeries(signal)
end_time = time.time()
print("Elapsed time naive algo %g seconds" % (end_time - start_time))
start_time = time.time()
fft(signal)
end_time = time.time()
print("Elapsed time of fft algo %g seconds" % (end_time - start_time))
```
循环卷积与 DFT
序列的循环移位
令 \(x[n]\) 表示有限长度的时域序列。
当我们进行 DFT 时，时域实际上是构建了信号的周期扩展。
因此，信号的时域移位实际上意味着信号的循环移位。
信号处理中最基本应用之一是线性卷积。
它可用于表示离散时间 LTI 系统的输出、相关和互相关、
滤波以及许多其他信号处理操作。
离散时间系统的线性卷积表示为：
信号 (x[n]) 和 (h[n]) 是有限长度的离散时间信号。
两个长度分别为 M 和 L 的序列的线性卷积会产生一个长度为 M+L-1 的序列。
一个相关的离散时间傅里叶变换性质是：
我们考虑两个序列，计算它们的 DFT，相乘，然后进行逆 DFT。
`时域的循环卷积导致频域中 DFT 系数的乘法。`
我们取序列 `[1,1,1,1,1]` 和 `[1,2,3,4]`。
The result of linear convolution is `[ 1 3 6 10 10 9 7 4 ]` The result of inverse of product of DFT coefficients is `[ 10. 10. 10. 10. 10.]`
我们对信号进行 N 点 DFT，N=max(M,L)。
...|1 1 1 1 1|1 1 1 1 1|1 1 1 1 1|.... ...|1 2 3 4 0|1 2 3 4 0|1 2 3 4 0|.... Result ...|10|10|10| ...
进行 DFT 会导致时域的周期性，现在我们像往常一样进行卷积。\(h[n-n_{o}]\) 将被循环移位。
...|1 1 1 1 1|1 1 1 1 1|1 1 1 1 1|.... ...|0 1 2 3 4|0 1 2 3 4|0 1 2 3 4|.... Result ...|10|10|10| ...
...|1 1 1 1 1|1 1 1 1 1|1 1 1 1 1|.... ...|4 0 1 2 3 |4 0 1 2 3|4 0 1 2 3|.... Result ...|10|10|10| ...
经过 N 次移位后，序列会重复。
因此，我们执行 N 次移位等同于序列在时域的周期进行卷积。
Circular convolution [ 10. 10. 10. 10. 10.]
循环移位是 DFT 引入的周期性的结果，并导致循环卷积运算。
我们可以将此归因于频率域采样引入的时间混叠。
采样定理指出，为了避免频率域的混叠，采样率必须大于被采样信号最高频率的两倍。
我们之前提到 DFT 是 DTFT 的采样版本。问题是如何提高采样率以避免时域混叠。
我们发现对序列进行零填充会导致傅里叶级数的采样点更密集。这相当于在频率域增加了 DTFT 的采样率。
这让我们直观地认为，如果我们对序列进行零填充，它将导致信号的 DTFT 采样率增加。
为了使循环卷积的结果等于线性卷积，我们只需将序列零填充到 M+L-1 的长度。
...|1 1 1 1 1 0 0 0 0 |1 1 1 1 1 0 0 0 0 | ...|1 2 3 4 0 0 0 0 0 |1 2 3 4 0 0 0 0 0 | result : 10 ...|1 1 1 1 1 1 0 0 0 |1 1 1 1 1 0 0 0 0 | ...|0 1 2 3 4 0 0 0 0 |0 1 2 3 4 0 0 0 0 | result : 10 ...|1 1 1 1 1 1 0 0 0 |1 1 1 1 1 0 0 0 0 | ...|0 0 1 2 3 4 0 0 0 |0 0 1 2 3 4 0 0 0 | result : 9
因此，我们可以看到，由于零填充，循环移位的分量为零，得到的结果等同于线性卷积结果。
上述示例的代码如下。
```
if mode == 7:
x=np.array([1,1,1,1,1]);
h=np.array([1,2,3,4,0]);
r=np.convolve(x,h)
print 'Linear convolution'
print r
#take 5 point DFT of the sequence
f1=fft(x,5);
f2=fft(h,5);
s=ifft(f1*f2);
print 'Circular convolution'
print abs(s)
#take 9 point DFT of sequence
f1=fft(x,9);
f2=fft(h,9);
s=ifft(f1*f2);
print 'zero padded Circular convolution'
print r
```
代码
所有示例的代码都可以在 GitHub 存储库中找到：
https://github.com/pi19404/pyVision/tree/master/pySignalProc/tutorial/fourierSeries.py
您可以在文件中将 mode 变量从 1-7 更改，以生成文章中显示的所有图。
文件 fourierSeries.py 也可以从以下链接下载：