用python代码实现实验：应用快速傅里叶变换对信号进行频谱分析

如下为作者所见题目，代码在题目下方，具体参数读者可自调 。

1.下面是自己写的FFT算法程序：

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
def butterfly_computation(x):
    N = len(x)
    if N == 1:
        return x
    even = butterfly_computation(x[::2])
    odd = butterfly_computation(x[1::2])
    factor = np.exp(-2j * np.pi * np.arange(N) / N)
    x = np.concatenate([even + factor[:N//2] * odd, even + factor[N//2:] * odd])
    return x
def fft(x, M):
    N = len(x)
    if N > M:
        raise ValueError("M should >= N")
    x = np.concatenate((x, np.zeros(M-N)))  # 在输入信号后面补零
    return butterfly_computation(x)

def ifft(x,M):
    N = len(x)
    if N > M:
        raise ValueError("M should >= N")
    x = np.concatenate((x, np.zeros(M-N)))  # 在输入信号后面补零
    return np.conj(butterfly_computation(np.conj(x)))/N

2.下面的代码实现对于正弦信号的频谱分析 

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd
import fft
matplotlib.use ("Qt5Agg") # 修改配置的后端 backend
N=256
M=16384
f=0.245 #采样频率10Hz
sinx=[]
Fs=8000
for i in range(N):
    a=np.sin(2*np.pi*f*i)
    sinx.append(a)
Xn=np.arange(0,N,1)
sinx=np.array(sinx)
plt.subplot(2,1,1)
plt.plot(Xn,sinx)
plt.title(f'SinTime,f={f},N={N},M={M},Fs={Fs}')
plt.xlabel("M")
plt.ylabel("sinx")

# 幅频响应
fft_sinx=fft.fft(sinx,M)
abs_b=np.abs(fft_sinx)
normolization_b=[]
for i in range(M):
    a=abs_b[i]/(M/2)
    normolization_b.append(a)
normolization_b[0]=normolization_b[0]/2
Omiga=[]
for i in np.arange(0,M,1):
    a=2*np.pi*Fs*(i)/M
    Omiga.append(a)
plt.subplot(2,1,2)
plt.plot(Omiga,normolization_b)
plt.title(f"sinFuPin,f={f},N={N},M={M},Fs={Fs}")
plt.xlabel("Omiga")
plt.ylabel("FuPin")
plt.show()
# # 相频响应
# angle_b=np.angle(fft_sinx,deg=True)
# plt.plot(Omiga,angle_b)
# plt.title(f"sinXiangPin,f={f},N={N},M={M}")
# plt.xlabel("Omiga")
# plt.ylabel("XiangWei")
# plt.show()

3.下面代码实现高斯序列的频谱分析

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import fft


# x=np.linspace(0,1,128)
# y=7*np.sin(2*np.pi*200*x) + 5*np.sin(2*np.pi*400*x)+3*np.sin(2*np.pi*600*x)

# 生成高斯信号
M=1024
N=32
p=32
q=1
Xa=[]
for i in range(N):
    z=np.exp(-(i-p)**2/q)
    Xa.append(z)
Xn=np.arange(0,N,1)
Xa=np.array(Xa)
plt.plot(Xn,Xa)
plt.title(f'GaussTime,p={p}')
plt.show()
Fs=2000# 设采样频率为Fs=2kHz
# 幅频响应
Omiga=[]
for i in np.arange(0,M,1):
    a=2*np.pi*Fs*(i)/M
    Omiga.append(a)
Xa_fft=fft.fft(Xa,M)
abs_Xa=np.abs(Xa_fft)
normolization_Xa=[]
for i in range(M):
    a=abs_Xa[i]/(M/2)
    normolization_Xa.append(a)
normolization_Xa[0]=normolization_Xa[0]/2
plt.plot(normolization_Xa)
plt.title(f"GaussFUPIN,p={p}")
# plt.xlabel("Omiga")
plt.ylabel("FuDu")
plt.show()
# # 相频响应
# angle_Xa=np.angle(Xa_fft,deg=True)
# plt.plot(angle_Xa)
# # plt.plot(angle_Xa)
# plt.title(f"GaussXIANGPIN,q={q}")
# # plt.xlabel("Omiga")
# plt.ylabel("XiangWei")
# plt.show()

4.下面实现对于衰减正弦序列的频谱分析

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import fft
N=32
M=256    #a=0.1,f=0.5625，N=32,M=256
f=0.5625 #采样频率10Hz
s=0.1
easinx=[]
for i in range(N):
    a=np.sin(2*np.pi*f*i)*np.exp(-s*i)
    easinx.append(a)
Xn=np.arange(0,N,1)
easinx=np.array(easinx)
plt.plot(Xn,easinx)
plt.title(f'eaSinTime,f={f},N={N},M={M}')
plt.xlabel("N")
plt.ylabel("easinx")
plt.show()
# 幅频响应
fft_easinx=fft.fft(easinx,M)
abs_b=np.abs(fft_easinx)
normolization_b=[]
for i in range(M):
    a=abs_b[i]/(M/2)
    normolization_b.append(a)
normolization_b[0]=normolization_b[0]/2
Omiga=[]
for i in np.arange(0,M,1):
    a=2*np.pi*f*(i)/M
    Omiga.append(a)
plt.plot(Omiga,normolization_b)
plt.title(f"easinFuPin,f={f},N={N},M={M}")
plt.xlabel("Omiga")
plt.ylabel("FuPin")
plt.show()
# 相频响应
angle_b=np.angle(fft_easinx,deg=True)
# plt.plot(Omiga,angle_b)
plt.plot(angle_b)
plt.title(f"easinXiangPin,f={f},N={N},M={M}")
plt.xlabel("Omiga")
plt.ylabel("XiangWei")
plt.show()

5.下面实现三角序列的频谱分析

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import fft
N=8
M=256
f=100
d=[]
for n in range(N):
    if(n<4):
        d.append(n+1)
    if(4<=n<=7):
        d.append(8-n)
Xd=np.arange(0,8,1)
plt.plot(Xd,d)
plt.title(f"TrangularTime,f={f},N={N},M={M}")
plt.xlabel("n")
plt.ylabel("Xd")
plt.show()
# 幅频响应
fft_d=fft.fft(d,M)
abs_d=np.abs(fft_d)
normolization_d=[]
for i in range(M):
    a=abs_d[i]/(M/2)
    normolization_d.append(a)
normolization_d[0]=normolization_d[0]/2
Omiga=[]
for i in np.arange(0,M,1):
    a=2*np.pi*f*(i)/M
    Omiga.append(a)
plt.plot(normolization_d)
plt.title(f"TrangularFUPIN,,f={f},N={N},M={M}")
# plt.xlabel("Omiga")
plt.ylabel("FuDu")
plt.show()
# 相频响应
angle_d=np.angle(fft_d,deg=True)
plt.plot(angle_d)
plt.title("TrangularXiangPin")
# plt.xlabel("Omiga")
plt.ylabel("XiangWei")
plt.show()

6.下面实现反三角序列的频谱分析

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import fft
N=8
M=8
f=100
d=[]
for n in range(N):
    if(n<4):
        d.append(4-n)
    if(4<=n<=7):
        d.append(n-3)
Xd=np.arange(0,8,1)
plt.plot(Xd,d)
plt.title(f"5TrangularTime,,f={f},N={N},M={M}")
plt.xlabel("n")
plt.ylabel("Xd")
plt.show()
# 幅频响应
fft_d=fft.fft(d,M)
abs_d=np.abs(fft_d)
normolization_d=[]
for i in range(M):
    a=abs_d[i]/(M/2)
    normolization_d.append(a)
normolization_d[0]=normolization_d[0]/2
Omiga=[]
for i in np.arange(0,M,1):
    a=2*np.pi*f*(i)/M
    Omiga.append(a)
plt.plot(normolization_d)
plt.title(f"5TrangularFUPIN,,f={f},N={N},M={M}")
# plt.xlabel("Omiga")
plt.ylabel("FuDu")
plt.show()
# 相频响应
angle_d=np.angle(fft_d,deg=True)
plt.plot(angle_d)
plt.title("5TrangularXiangPin")
# plt.xlabel("Omiga")
plt.ylabel("XiangWei")
plt.show()

 如有不足，欢迎指正。