本文是接上一篇结尾的问题，“当数据长度非2时怎么办？”
废话不多说，先上得出的结论（仅限个人观点！）
各种尝试办法得出的结果如图：
我分别用numpy自带的fft.fft作为参考，然后用原始DFT和CT优化的DFT两种方法与参考进行对比，最终得到的就是上图；
图看起来只有一条线，好像我是拿一个结论重复画了好多次一样。
然而并不是，恰恰说明了各个办法得到的结果是有多么的接近，perfect~~
误差：1e-19
介绍：
因为工作原因，每天都在np.fft.fft、plt.psd，天天在敲这俩函数，敲得越多心就越虚，咋搞？ 一个字拆！！！
一开始是拆psd，搞了一堆代码，参考fft与psd的关系【傅里叶变化求功率谱】；
然后是拆 fft，对于理论，基本是个理工科生，大部分都略知一二，至于从离散傅里叶级数 ----嘟嘟— 的推导出离散傅里叶变化，也不是本文重点，哪天突然热血沸腾了再慢慢把这个过程敲出来吧（注意：我写了两个嘟嘟，不是三个也不是四个，品一品）
回归本文：
（1）第一步：直接用离散傅里叶变化的公式实现，以np.fft.fft做对比，从二者误差可以看出，np.fft的结果就是对这个过程进行优化~
# 手写DFT
def Fun_DFT(data):
data = np.array(data)
data_len = data.shape[0]
data_fft = np.zeros_like(data)
# 循环的公共部分
W_n = -1j * 2 * np.pi * np.arange(data_len) / data_len
for n in range(data_len):
if n % 1000 == 0:
print('=====正在第{:}个点的DFT结果====='.format(n))
data_fft_temp = data * np.exp(W_n)**n
data_fft_i = np.sum(data_fft_temp)
data_fft[n] = data_fft_i
return data_fft
（2）第二步：第一步的结果是用for循环写的，旋转因子的计算尽量提到循环外了，但计算时间挺长的，受不了！八千多个点计算了这么长时间，就考虑用numpy的矩阵实现乘法，结果时间上确实节省了接近10s；（计算结果懒得贴了，和上图一模一样）
# 优化手写DFT：耗时的原因是因为for，改成矩阵计算提升运算速度
# 对旋转因子矩阵的计算进行优化
def Fun_Cac_Twiddle_Fac_dot_plus(N,P,Q):
val_temp = np.arange((P-1)*(Q-1)+1) % N # 没有必要对N取余，以为P*Q的最大值也就是N，这里加上N是为了和之前代码对称
twiddle_fac_index = np.zeros([P,Q])
min_size, max_size = min(P,Q), max(P,Q)
for i in range(1, min_size):
twiddle_fac_index[i, :] = val_temp[::i][:max_size]
twiddle_fac_dot = np.exp(-1j * 2*np.pi * twiddle_fac_index / N)
return twiddle_fac_dot
def Fun_DFT_plus(data):
data = np.array(data)
N = data.shape[0]
W_n = Fun_Cac_Twiddle_Fac_dot_plus(N, N, N)
data_dft = np.dot(data, W_n)
return data_dft
（3）第三步：虽然第二步用矩阵实现的方式，确认节省了很多时间，但计算8000多个点的dft需要用5s，还是受不了，然后就引出了本文的重点——CT
理论推导见下图：
（写的比较细，懒得转成LaTeX公式了，也许哪天无聊，就突然把这一堆转成在线公式了）
ps：为了大家看的清楚，专门把乱七八糟的草稿整理了一下，小编不易给个好评谢谢~~~
另外感谢公众号“88A写字的地方”博主解答的疑惑，感兴趣的童鞋可以关注
代码内有详细注释
# 按照公式计算第r*s点的fft结果（用时118s）
def DFT_CT_main(data, Q, P, r, s):
N = data.shape[0]
data_dft_outer = 0
for q in range(Q): # 取每块里面的每个值
# 内存循环
data_dft_inter = 0
for p in range(P):
data_pi = data[p*Q+q]
# q*s就是旋转项，对N取余是因为：对于旋转因子e**(2pi/N * index)，
# 当index=(N+x)和(N-x)时，旋转因子的结果是一致的，公式推导见图片
twiddlle_fac_inter_index = p*s % P
twiddlle_fac_interx = np.exp(-1j * 2 * np.pi * twiddlle_fac_inter_index / P) # 考虑提前存储P个值
data_dft_inter += data_pi * twiddlle_fac_interx
# 乘以单独的一个旋转因子
index_middle = q * s % N
twiddle_fac_middle = np.exp(-1j * 2 * np.pi * index_middle / N) # 考虑提前存储N个值
data_dft_middle = data_dft_inter * twiddle_fac_middle
# 外层循环
index_outer = q * r % Q
twiddle_fac_outer = np.exp(-1j * 2 * np.pi * index_outer / Q) # 考虑提前存储N个值
data_dft_outer += data_dft_middle * twiddle_fac_outer
return data_dft_outer
def DFT_CT(data, P, Q):
CT_data = []
for r in range(P):
if r % 50 == 0:
print('========当前计算位置：{:}========'.format(r))
for s in range(Q):
CT_data.append(DFT_CT_main(data, P, Q, r, s))
return CT_data
（4）第四步：第三个用推导的公式优化后，代码的用时竟然118s，这就更受不了了，进一步优化。耗时的原因无非就是for用的太多，但是这一步先暂时不优化它，因为在写第三步代码的时候，看到旋转因子的计算是重复计算了一堆，所以考虑先优化旋转因子，提到循环外面。
结论也是相对可喜可贺的，用了45s，节省的时间确实挺多的~~
# 优化1：每次dft循环的时候，计算旋转因子是可以提前存储的，因为无论何时，twiddle_fac的分子index永远是小于分母N的，
# 所以提前便利存储i=[0:N]的e**(-1j*2pi*i) （用时45s）
# 旋转因子矩阵计算
def Fun_Cac_Twiddle_Fac(N):
# twiddle_fac = np.zeros(N, 'complex')
# for index in range(N):
# twiddle_fac[index] = np.exp(-1j * 2*np.pi * index / N)
twiddle_fac_index = np.arange(N)
twiddle_fac = np.exp(-1j * 2*np.pi * twiddle_fac_index / N)
return twiddle_fac
def DFT_CT_plus1_main(data, Q, P, r, s, twiddlle_fac_inter_all, twiddle_fac_middle_all, twiddle_fac_outer_all):
N = data.shape[0]
data_dft_outer = 0
for q in range(Q): # 取块里面的每个值
# 内存循环
data_dft_inter = 0
for p in range(P): # 取每块里的第q个值
data_pi = data[p * Q + q]
index_inter = p * s % P
twiddlle_fac_inter = twiddlle_fac_inter_all[index_inter] # 读取提前存储P个值
data_dft_inter += data_pi * twiddlle_fac_inter
# 乘以单独的一个旋转因子
index_middle = q * s % N
twiddle_fac_middle = twiddle_fac_middle_all[index_middle] # 读取提前存储N个值
data_dft_middle = data_dft_inter * twiddle_fac_middle
# 外层循环
index_outer = q * r % Q
twiddle_fac_outer = twiddle_fac_outer_all[index_outer] # 读取提前存储N个值
data_dft_outer += data_dft_middle * twiddle_fac_outer
return data_dft_outer
def DFT_CT_plus1(data, Q, P):
res = []
N = data.shape[0]
# 导入提前计算好的旋转因子
twiddlle_fac_inter_all = Fun_Cac_Twiddle_Fac(P)
twiddle_fac_middle_all = Fun_Cac_Twiddle_Fac(N)
twiddle_fac_outer_all = Fun_Cac_Twiddle_Fac(Q)
for r in range(Q):
if r % 50 == 0:
print('========当前计算位置：{:}========'.format(r))
for s in range(P):
res.append(DFT_CT_plus1_main(data, Q, P, r, s, twiddlle_fac_inter_all,
twiddle_fac_middle_all, twiddle_fac_outer_all))
return res
（5）第五步：这一步就考虑对for循环进行优化，用矩阵乘的方式节省for循环耗费的地址跳转时间。
怎么转换成矩阵实现的方式？
按我比较蠢得办法是先手动写两次循环，把这两次循环转成矩阵乘，就可以得到最终的办法，上图~
然后站在上帝视角看这个转换成矩阵的过程，原始的公式里一共有p、q、r、s四个变量，两个求和公式，固定一个变量，两个求和公式其实就是两个矩阵相乘；
比如当固定r时，内存循环就是矩阵的某一行(s行) * 另一个矩阵的不同列(q取不同值)，得到的是一个矩阵，然后固定p，外层循环就是矩阵的某一行(r行) * 另一个矩阵的不同列(s取不同值)，最终得到一个r*s大小的矩阵，矩阵的每个位置分别表示第rs位置的dft结果；
再站的高一点：四个变量，两个矩阵相乘会消耗掉一个变量，再乘以一个矩阵就又会消耗掉一个变量，最终得到两个变量~
代码：
# 优化2：优化1速度慢的原因是因为调用了一堆循环，所以优化的方法是把循环变成矩阵
# 首先对于PQ的循环，其实是在固定rs的基础下的循环，所以rs固定后，PQ的循环是一个矩阵乘；
# 其次对于RS的循环，相当于在PQ循环后再乘以一个矩阵（或者可以理解为固定PQ的基础下的循环）
# 旋转因子矩阵计算
def Fun_Cac_Twiddle_Fac_dot(N,P,Q):
twiddle_fac_index = np.zeros([P,Q], 'complex')
for q in range(Q):
for s in range(P):
# q*s就是旋转项，对N取余是因为：对于旋转因子e**(2pi/N * index)，
# 当index=(N+x)和(N-x)时，旋转因子的结果是一致的，公式推导见图片
index_ij = q * s % N
twiddle_fac_index[q, s] = index_ij
twiddle_fac = np.exp(-1j * 2*np.pi * twiddle_fac_index / N)
return twiddle_fac
#
# twiddle_fact_test_data = Fun_Cac_Twiddle_Fac_dot(4,4,4)
# 对旋转因子矩阵的计算进行优化
def Fun_Cac_Twiddle_Fac_dot_plus(N,P,Q):
val_temp = np.arange((P-1)*(Q-1)+1) % N # 没有必要对N取余，以为P*Q的最大值也就是N，这里加上N是为了和之前代码对称
twiddle_fac_index = np.zeros([P,Q])
min_size, max_size = min(P,Q), max(P,Q)
for i in range(1, min_size):
twiddle_fac_index[i, :] = val_temp[::i][:max_size]
twiddle_fac_dot = np.exp(-1j * 2*np.pi * twiddle_fac_index / N)
return twiddle_fac_dot
def DFT_CT_plus2(data,P,Q):
data = np.array(data)
N = data.shape[0]
data_x_dot = data.reshape([P, Q])
W_fac_dot = Fun_Cac_Twiddle_Fac_dot_plus(P,P,P)
M_fac_dot = Fun_Cac_Twiddle_Fac_dot_plus(N,Q,P)
N_fac_dot = Fun_Cac_Twiddle_Fac_dot_plus(Q,Q,Q)
data_dft_inter = np.dot(W_fac_dot.T, data_x_dot)
data_y_dot = data_dft_inter * M_fac_dot.T # 注意这里是点乘
data_dft_outer = np.dot(N_fac_dot, data_y_dot.T)
data_dft = data_dft_outer.reshape(-1)
return data_dft
最终：
最后就是把所有方法同时调用一遍，把结果画到一个图里，就得到一开始的图；
但是，该说不说，用CT矩阵实现的方式竟然比np.fft.fft的时间还要短很多，至于原因怀疑是np.fft.fft兼容了更多的功能，毕竟括号里面有一堆关键字；
def Fun_Cac_Time(func,data, org_method=True, P=1, Q=1):
time_start = time.time()
if org_method:
res = func(data)
else:
res = func(data, P, Q)
data_fft = np.roll(res, len(res) // 2)
data_db = cac_db(data_fft)
plt.plot(data_db, '*-')
plt.grid(True)
time_end = time.time()
return time_end - time_start, data_db
# %%
if __name__ == '__main__':
# 构造数据
noise_data = np.random.uniform(low=-1, high=1, size=8192) + np.random.uniform(low=-1, high=1, size=8192) * 1j
flt_coef = sig.remez(150, [0, 50, 100, 491.52], [1, 0], fs=983.04)
data = np.convolve(noise_data, flt_coef)
# python第三方库
time_numpy, data_numpy = Fun_Cac_Time(np.fft.fft, data)
print('numpy自带FFT用时：{:}'.format(time_numpy))
# DFT公式for实现
time_dft_for, data_dft_for = Fun_Cac_Time(die.Fun_DFT, data)
print('原始DFT for循环实现用时：{:}'.format(time_dft_for))
# DFT公式矩阵实现
time_dft_dot, data_dft_dot = Fun_Cac_Time(die.Fun_DFT_plus, data)
print('原始DFT 矩阵实现用时：{:}'.format(time_dft_dot))
# CT优化DFT for实现
P, Q = Fun_Cac_CT_PQ(len(data))
time_ct_for, data_ct_for = Fun_Cac_Time(DFT_CT, data, org_method=False, P=P, Q=Q)
print('CT优化DFT for循环实现用时：{:}'.format(time_ct_for))
# CT优化DFT:旋转因子提前计算
P, Q = Fun_Cac_CT_PQ(len(data))
time_ct_for1, data_ct_for1 = Fun_Cac_Time(DFT_CT_plus1, data, org_method=False, P=P, Q=Q)
print('CT优化DFT for循环旋转因子提前计算实现用时：{:}'.format(time_ct_for1))
# CT优化DFT：矩阵实现
P, Q = Fun_Cac_CT_PQ(len(data))
time_ct_dot, data_ct_dot = Fun_Cac_Time(DFT_CT_plus2, data, org_method=False, P=P, Q=Q)
print('CT优化DFT 矩阵实现用时：{:}'.format(time_ct_dot))
# # FFT
# time_fft, data_fft = Fun_Cac_Time(die.Fun_FFT, data)
# print('fft 实现用时：{:}'.format(time_fft))
plt.legend(['numpy自带FFT-{:.2f}s'.format(time_numpy),
'原始DFT for循环实现-{:.2f}s'.format(time_dft_for),
'原始DFT 矩阵实现-{:.2f}s'.format(time_dft_dot),
'CT优化DFT for循环实现-{:.2f}s'.format(time_ct_for),
'CT优化DFT for循环旋转因子提前-{:.2f}s'.format(time_ct_for1),
'CT优化DFT 矩阵实现-{:.2f}s'.format(time_ct_dot),
# 'FFT 实现-{:.2f}s'.format(time_ct_dot)
])
疑问解答：
问题1：DFT原始公式用矩阵实现和 DFT用CT优化后的矩阵实现对比，后者为什么计算会快那么多？
答：观察CT矩阵实现最后一个矩乘的计算过程，可以发现矩阵W的一行与矩阵Y的每一列都进行了相乘，分别得到最终结果的前Q个结果，说明前Q个计算结果共用了W一行的结果；而原始公式矩阵实现计算前Q个值时并没有公用，直接拿N个点相乘，这就是省计算量的地方；
问题2：直接拿原始8314长度数据画出来的频谱和截取8192画出来的差距为什么这么大？
答：截取8192相当于在时域上加了一个带内8192总长8314长度的矩形窗，
时域上进行乘积操作，相当于在频域上用矩形窗的频谱和原数据的频谱进行了卷积，就会把原信号带内部分的能量带到带外，所以带外信号会抬升,仿真代码见下：
# 参考结果
data_org_fft = np.fft.fft(data[:8192])
data_org_fft = np.roll(data_org_fft, len(data_org_fft)//2)
plt.plot(cac_db(data_org_fft))
plt.grid(True)
# 生成时域上的矩形窗
win_time = np.zeros_like(data)
win_time[:8192] = 1
data_win_fft = np.fft.fft(win_time)
data_win_fft = np.roll(data_win_fft, len(data_win_fft)//2)
# plt.plot(data_win_fft/max(data_win_fft), '*-')
# 原始信号
data_org_fft = np.fft.fft(data)
data_org_fft = np.roll(data_org_fft, len(data_org_fft)//2)
plt.plot(cac_db(data_org_fft))
plt.grid(True)
# 频域上进行卷积(窗系数除sum是为了不把滤波器的增益带到结果中)
con_res = np.convolve(data_org_fft, data_win_fft/sum(data_win_fft), 'same')
plt.plot(cac_db(con_res))
plt.grid(True)
眼看千遍不如手过一遍，只用眼睛看，永远是别人的~
ps：问题永远是一个接一个，当你提出一个问题，然后在解决它的时候，又会提出不止一个问题，这是一个递归，甚至是一个向下的递归，最终也不知道会怎么样~~~~~~ 就先让它迭着吧！！
预告：下一期想实现区分男女生声音的功能~敬请期待，谢谢思密达！！
1782
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
