Python实现DFT

1.利用公式实现DFT:

参阅书上的公式定义可知，可利用Wn展开矩阵的形式实现；

# coding=GBK
import numpy as np
from numpy import arange, sin, pi, cos
import matplotlib.pyplot as plt

fs=8192
Ts=1/fs 
N=64    
n = np.arange(N)  
xn1=cos(2*pi*512.*n*Ts)+0.75*cos(2*pi*1024.*n*Ts)
#print (type(xn1))
#画出时域采样图
plt.subplot(311)
plt.stem(n,xn1,use_line_collection=True)
plt.xlabel('n')
plt.ylabel('xn')
plt.title('Time domain sampling')
my_x_ticks = np.arange(0, 70, 5)
my_y_ticks = np.arange(0, 3, 1)
plt.xticks(my_x_ticks)
plt.yticks(my_y_ticks)

#DFT
n = np.expand_dims(n, axis=1)         
k = n  
m = n.T * k / N                     
S = np.exp(1j * 2 * pi *m)
S=S.conjugate(S)
X2 = np.dot(xn1, S)

mX = np.abs(X2)
pX = np.angle(X2)

# plot the magnitude and phase
plt.subplot(312)
plt.stem(mX,use_line_collection=True)
plt.xlabel('k');
plt.ylabel('|X(K)|')
plt.title('Amplitude frequency characteristic')
my_x_ticks = np.arange(0, 70, 5)
my_y_ticks = np.arange(0, 50, 10)
plt.xticks(my_x_ticks)
plt.yticks(my_y_ticks)

plt.subplot(313)
plt.stem(n,pX,use_line_collection=True)
plt.xlabel('k');
plt.ylabel('φ(w)/rad');
plt.title('Phase frequency characteristics')
my_x_ticks = np.arange(0, 70, 5)
my_y_ticks = np.arange(-4, 5, 1)

该代码，将xn1先抽样成序列，之后进行DFT变换。