非2的幂，离散傅里叶变化DFT的快速实现

本文是接上一篇结尾的问题，“当数据长度非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_lenfor n in range(data_len):if n % 1000 == 0:print('=====正在第{:}个点的DFT结果====='.format(n))data_fft_temp = data * np.exp(W_n)**ndata_fft_i = np.sum(data_fft_temp)data_fft[n] = data_fft_ireturn 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_dotdef 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 = 0for q in range(Q):  # 取每块里面的每个值# 内存循环data_dft_inter = 0for 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 % Ptwiddlle_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 % Ntwiddle_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 % Qtwiddle_fac_outer = np.exp(-1j * 2 * np.pi * index_outer / Q)  # 考虑提前存储N个值data_dft_outer += data_dft_middle * twiddle_fac_outerreturn data_dft_outerdef 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_facdef 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 = 0for q in range(Q):  # 取块里面的每个值# 内存循环data_dft_inter = 0for p in range(P):  # 取每块里的第q个值data_pi = data[p * Q + q]index_inter = p * s % Ptwiddlle_fac_inter = twiddlle_fac_inter_all[index_inter]  # 读取提前存储P个值data_dft_inter += data_pi * twiddlle_fac_inter# 乘以单独的一个旋转因子index_middle = q * s % Ntwiddle_fac_middle = twiddle_fac_middle_all[index_middle]  # 读取提前存储N个值data_dft_middle = data_dft_inter * twiddle_fac_middle# 外层循环index_outer = q * r % Qtwiddle_fac_outer = twiddle_fac_outer_all[index_outer]  # 读取提前存储N个值data_dft_outer += data_dft_middle * twiddle_fac_outerreturn data_dft_outerdef 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 % Ntwiddle_fac_index[q, s] = index_ijtwiddle_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_dotdef 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) * 1jflt_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：问题永远是一个接一个，当你提出一个问题，然后在解决它的时候，又会提出不止一个问题，这是一个递归，甚至是一个向下的递归，最终也不知道会怎么样~~~~~~ 就先让它迭着吧！！

预告：下一期想实现区分男女生声音的功能~敬请期待，谢谢思密达！！

非2的幂，离散傅里叶变化DFT的快速实现相关推荐

如何用C语言做离散傅里叶变化
1 离散傅里叶变换此图大家都再熟悉不过了,任何周期信号都可以看作是一系列不同频率正弦信号叠加的结果. 信号从时域变换到频域的过程,可以借助傅里叶变换(Fourier transform)完成.我们知 ...
各种FIFO硬件设计（FIFO概念、异步、同步、非2次幂深度FIFO）
文章目录一.FIFO概述二.FIFO分类三.FIFO重要信号与参数 3.1 信号 3.2 参数 3.2.1 data_depth的确定四.FIFO存储原理五.同步FIFO 5.1 空满信号判 ...
任意深度同步FIFO设计总结（非2次幂）
同步FIFO属于是最基本的模块之一了,但是对于任意深度(非2次幂与2次幂均可以)的同步FIFO可能有些难度,但是只要了解格雷码的性质就能解决.接下来将介绍整个设计的思考流程,不单单是给出解决方法. 对 ...
小试牛刀--我的快速离散傅里叶变化matlab函数（FFT）
小试牛刀--我的快速离散傅里叶变化matlab函数(FFT) 想法来源例子分析函数说明代码展示想法来源我在实验室里使用示波器观察实验数据的波形时,有时需要对实验数据进行傅里叶变换,观察其中的 ...
蒙哥马利幂模算法（二分快速幂）
蒙哥马利幂模算法(二分快速幂) 用于加快幂次取模运算的速度. 形式为这种ab%c 可以等于 (ab1%c)*(ab2%c)%c,其中b=b1+b2:如果b1=b2就更简便一些,(ab/2%c)2%c, ...
离散傅里叶变换 (DFT)、快速傅里叶变换 (FFT)
目录离散傅里叶变换 (DFT) 离散傅里叶变换的基离散傅里叶变换快速傅里叶变换 (FFT) 卷积线性时不变系统傅里叶级数参考文献离散傅里叶变换 (DFT) 离散傅里叶变换的基对于周期为 ...
傅里叶变换、离散傅里叶变换(DFT)、快速傅里叶变换(FFT)详解
前置知识以下内容参考<复变函数与积分变换>,如果对积分变换有所了解,完全可以跳过忽略复数的三角表达式如下 Z=r(cosθ+isinθ)Z=r(cos\theta+isin\theta ...
【DSP数字信号处理学习笔记】—— 详细推导DFT的快速实现算法：FFT 基于库利-图基算法的实现
引言:尽管离散傅里叶变换(DFT)让频谱分析技术在计算机上的实现成为可能,但是受限于DFT算法庞大的计算量 O(N2)O(N^2)O(N2),使得DFT在一开始并没有被广泛使用,直到快速傅里叶变换算法 ...
mysql中如何幂次方的函数_幂次方的四种快速取法（不使用pow函数）
Pow(x, n) 方法一:暴力法方法二:递归快速幂算法方法三:迭代快速幂算法方法四:位运算法方法一:暴力法思路只需模拟将 x 相乘 n 次的过程. 如果 \(n < 0\),我们可 ...

非2的幂，离散傅里叶变化DFT的快速实现

非2的幂，离散傅里叶变化DFT的快速实现相关推荐

最新文章

热门文章