python 数字信号产生并求其FFT

python 数字信号产生并求其FFT

本文写了一篇关于数字信号发生和做FFT处理的文章。

信号产生模块

这个模块用于产生单位阶跃信号、单位脉冲信号、简单的锯齿波信号。

也就是常用的数字信号处理中常用的序列。
话不多说，直接上代码：

signal_module.py:

import numpy as np

def saw_tooth(n0,n1,n2):
	# Generates x(n)=saw_tooth
	# -----------------------------------------
	# n,x = saw_tooth(n0, n1, n2)  n1 <= n0, n0 <= n2
	num=0
	if (n0<n1 or n2<n1 or n2<n0):
		raise Exception('arguments must satisfy n1<=n0<=n2')
	else:
		n=np.arange(n1,n2+1)
		x=np.zeros(len(n))
		no = np.arange(len(n))
		for a in no:
			if a < n0-n1:
				x[a]=0
			else:
				x[a]=num
				num+=1
	return n,x

def step_seq(n0,n1,n2):
	# Generates
	# x(n) = u(n - n0);
	# n1 <= n0, n0 <= n2
	if (n0<n1) or (n2<n1) or (n2<n0):
		raise Exception('arguments must satisfy n1<=n0<=n2')
	else:
		n=np.arange(n1,n2+1)
		x=(n>=n0)*1.
	return n,x

def imp_seq(n0,n1,n2):
	# Generates x(n)=delta(n-n0);n1<=n0,n0<=n2
	# -----------------------------------------
	# [x, n] = impseq(n0, n1, n2)
	if (n0<n1) or (n2<n1) or (n2<n0):
		raise Exception('arguments must satisfy n1<=n0<=n2')
	else:
		n=np.arange(n1,n2+1)
		x=(n==n0)*1.
	return n,x

if __name__ =="__main__":
	from matplotlib import pyplot as plt
	n1,x1=saw_tooth(5,1,10)
	n2,x2=imp_seq(5,1,10)
	n3,x3=step_seq(5,1,10)
	plt.subplot(3, 1, 1)
	plt.stem(n1, x1, use_line_collection=True)
	plt.title("sawtooth")
	plt.xlabel("n")
	plt.subplot(3, 1, 2)
	plt.stem(n2, x2, use_line_collection=True)
	plt.title("impseq")
	plt.xlabel("n")
	plt.subplot(3, 1, 3)
	plt.stem(n3, x3, use_line_collection=True)
	plt.title("stepseq")
	plt.xlabel("n")
	plt.show()

运行结果：

快速傅里叶变换

这里参考了：https://blog.csdn.net/qq_27825451/article/details/88553441

import numpy as np
from scipy.fftpack import fft, ifft
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl
from signal_module import *

mpl.rcParams['font.sans-serif'] = ['SimHei']  # 显示中文
mpl.rcParams['axes.unicode_minus'] = False  # 显示负号

n, y1 = imp_seq(40, 0, 99)
n, y2 = imp_seq(50, 0, 99)
n, y3 = imp_seq(60, 0, 99)
n, y4 = imp_seq(70, 0, 99)
y=y1+y2+y3+y4                 #4个序列相加
fft_y = fft(y)  # 快速傅里叶变换

N = 100
x = np.arange(N)  # 频率个数
half_x = x[range(int(N / 2))]  # 取一半区间

abs_y = np.abs(fft_y)  # 取复数的绝对值，即复数的模(双边频谱)
angle_y = np.angle(fft_y)  # 取复数的角度
normalization_y = abs_y / N  # 归一化处理（双边频谱）
normalization_half_y = normalization_y[range(int(N / 2))]  # 由于对称性，只取一半区间（单边频谱）

plt.subplot(231)
plt.plot(n, y)
plt.title('原始波形')

plt.subplot(232)
plt.plot(x, fft_y, 'black')
plt.title('双边振幅谱(未求振幅绝对值)', fontsize=9, color='black')

plt.subplot(233)
plt.plot(x, abs_y, 'r')
plt.title('双边振幅谱(未归一化)', fontsize=9, color='red')

plt.subplot(234)
plt.plot(x, angle_y, 'violet')
plt.title('双边相位谱(未归一化)', fontsize=9, color='violet')

plt.subplot(235)
plt.plot(x, normalization_y, 'g')
plt.title('双边振幅谱(归一化)', fontsize=9, color='green')

plt.subplot(236)
plt.plot(half_x, normalization_half_y, 'blue')
plt.title('单边振幅谱(归一化)', fontsize=9, color='blue')

plt.show()

运行结果：