python和傅里叶级数 FFT频谱分析

上学的时候玩心很重，信号与系统课程学的很差。这段时间在学习python的FFT，重新把课程内容捡起来看了看，感触颇多。记录如下：

一、理论准备：

对一个普通信号而言，可能会有直流分量；信号可能呈现奇函数性质，也可能呈现偶函数性质。故可将一个普通信号表示为直流分量和一系列正弦函数、余弦函数叠加的形式。

    

令

    

根据欧拉公式，将三角形式的傅里叶级数转化为式指数形式的傅里叶级数

    

为方便和fft算法对比，需要找到2个教材上已经分析完的周期信号，最好是一个只有正弦分量、一个只有余弦分量。目标分析函数如下：

信号1：占空比为50%的方波，偶函数。显然，该信号为奇谐信号，且是偶信号，故分解后只有奇次余弦

    

信号2：奇函数锯齿波。由于函数为奇函数，故只有正弦分量。分解后的信号如下：

   

二、python代码和分析

python代码如下：(以信号2为例)

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

def showfft(y):
    for i in range(len(y)):
        absy=np.round(np.abs(y[i])/len(y),3)
        angley=np.round(np.angle(y[i]),3)
        realy=np.round(y[i].real/len(y),3)
        xandib="\t"+str(realy)+"+j*"+str(np.round(y[i].imag/len(y),3))
        print(i,absy,"∠",angley,xandib)
    return

plt.figure(figsize=(9,6), dpi=200)
plt.figure(1)
x= np.linspace(0,1,100)
y1=12*np.sin(2*np.pi*2*x) -6*np.sin(2*np.pi*4*x)+4*np.sin(2*np.pi*6*x)-3*np.sin(2*np.pi*8*x)+2.4*np.sin(2*np.pi*10*x)

ax1=plt.subplot(311)
ax1.plot(x,y1)
ffty=fft(y1)
showfft(ffty)

fn=np.arange(len(x))
ax2=plt.subplot(312)
ax2.plot(fn,np.abs(ffty))

ax3=plt.subplot(313)
ax3.plot(fn,np.angle(ffty))
plt.show()

图形如下：

fft输出如下：

0 0.0 ∠ -0.0     0.0+j*-0.0
1 0.068 ∠ -1.539     0.002+j*-0.068
2 5.938 ∠ -1.508     0.373+j*-5.927
3 0.226 ∠ 1.665     -0.021+j*0.225
4 3.02 ∠ 1.696     -0.379+j*2.996
5 0.17 ∠ -1.414     0.027+j*-0.168
6 1.978 ∠ -1.382     0.371+j*-1.943
7 0.195 ∠ 1.791     -0.043+j*0.191
8 1.489 ∠ 1.822     -0.37+j*1.442
9 0.195 ∠ -1.288     0.055+j*-0.188
10 1.208 ∠ -1.257     0.373+j*-1.149
 

显然，系统输出为 a- jb 的形式。b为正弦信号幅值；a为余弦信号幅值。

对 2hz 信号，fft结果余弦信号幅值为0.373；正弦信号幅值为 b = 5.927。显然此时相角应为正。即求相角时应将虚部乘以 -1 ，为照顾习惯，应实现 -3.1415926 ~ 3.1415926  到-180 ~ 180的映射工作。修改后的相位谱如下：