python实现DSP中相关函数及应用

如下为作者所见题目： 

 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.pyplot as plt
import pandas as pd
# 线性卷积函数
def convolve(signal1, signal2):
    conv_result = []
    len1 = len(signal1)
    len2 = len(signal2)
    for i in range(len1 + len2 - 1):
        conv_sum = 0
        for j in range(len1):
            if i - j >= 0 and i - j < len2:
                conv_sum += signal1[j] * signal2[i - j]
        conv_result.append(conv_sum)
    return conv_result

3.对高斯信号生成自相关函数

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

import conv
import fft
import matplotlib
matplotlib.use ("Qt5Agg") # 修改配置的后端 backend
# 生成高斯信号
M=512
N=32
p=16
q=30
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.figure(1)
plt.subplot(2,1,1)
plt.plot(Xn,Xa)
plt.title(f'GaussTime,p={p},N={N},q={q}')
# 归一化处理：求高斯信号的DFT，求其对应的绝对值abs_a,
# 将abs_a[0]除以N，将其余的abs_a除以N/2
dft_a=fft.fft(Xa,M)
abs_a=np.abs(dft_a)
norm_a=[]
for i in range(M):
    a=abs_a[i]/(N/2)
    norm_a.append(a)
norm_a[0]=norm_a[0]/2
Xm=np.arange(M)
plt.subplot(2,1,2)
plt.plot(Xm,norm_a)
plt.title("GaussFUPIN")
plt.show()
plt.figure(2)
# 时域求自相关
# 利用numpy自带的correlate函数求解自相关
# 利用FFT算法求其功率谱，并归一化
plt.subplot(2,1,1)
# cor_a=np.correlate(Xa,Xa,"full")
cor_a=conv.convolve(Xa,Xa)#调用写的线性卷积函数conv()
X_a=np.arange(-N+1,N,1)
plt.plot(X_a,cor_a)
plt.title(f"GaussCorrelate,p={p},N={N},q={q}")

dft_cor=fft.fft(cor_a,M)
abs_cor=np.abs(dft_cor)
norm_cor=[]
for i in range(len(abs_cor)):
    a=abs_cor[i]/(len(abs_cor)/2)
    norm_cor.append(a)
norm_cor[0]=norm_cor[0]/2
plt.subplot(2,1,2)
plt.plot(norm_cor)
plt.title("corFUPIN")
plt.show()

4.对正弦信号生成自相关函数

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import fft
import matplotlib
matplotlib.use ("Qt5Agg") # 修改配置的后端 backend
# 生成sin信号
M=256
N=128
f=0.0625
Fs=8000
Xa=[]
for i in range(N):
    z=np.sin(2*np.pi*f*i)
    Xa.append(z)
Xn=np.arange(0,N,1)
Xa=np.array(Xa)
plt.figure(1)
plt.subplot(2,1,1)
plt.plot(Xn,Xa)
plt.title(f'SinTime,f={f},N={N},M={M}')
# 归一化处理：求高斯信号的DFT，求其对应的绝对值abs_a,
# 将abs_a[0]除以N，将其余的abs_a除以N/2
dft_a=fft.fft(Xa,M)
abs_a=np.abs(dft_a)
norm_a=[]
for i in range(M):
    a=abs_a[i]/(N/2)
    norm_a.append(a)
norm_a[0]=norm_a[0]/2
# Xm=np.arange(M)
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,norm_a)
plt.xlabel("Omiga")
plt.title(f"SinFUPIN,f={f},N={N},M={M},Fs={Fs}")
plt.show()




# 频域求自相关
# 利用FFT算法求其功率谱，并归一化
plt.figure(2)
fft_a=fft.fft(Xa,M)
fft_aa=[]
for i in range(M):
    a=fft_a[i]**2
    fft_aa.append(a)
ifft_a=fft.ifft(fft_aa,M)



# 频域求自相关
plt.subplot(2,1,1)
X_a=np.arange(len(ifft_a))

plt.plot(X_a,ifft_a)
plt.title(f"SinCorrelate,f={f},N={N},M={M}'")

dft_cor=fft.fft(ifft_a,M)
abs_cor=np.abs(dft_cor)
norm_cor=[]
for i in range(len(abs_cor)):
    a=abs_cor[i]/(len(abs_cor)/2)
    norm_cor.append(a)
norm_cor[0]=norm_cor[0]/2
plt.subplot(2,1,2)
plt.plot(norm_cor)
# plt.xlabel("corFUPIN")
plt.title("corFUPIN")
plt.show()

 欢迎在评论区留言交流。