dget函数的批量使用方法,dget函数能否下拉填充

首页 > 经验 > 作者:YD1662022-11-07 23:08:34

由于W(N,n,k)与W(N,n,K N/2)的关系,当n为欧时,两者相等。当n为奇数时,两者相反。

又由于:F(k) 的关于 W(N,n,k)的for循环,而F(k N/2)是关于W(N,n,k N/2)的for 循环,所以,我们用n的奇偶行将将F(k)的实现拆为两部分:

dget函数的批量使用方法,dget函数能否下拉填充(5)

拆成n为偶数,n为奇数两部分后,关键就是求出每个部分的值了,我么将F(k)表示为:

F(k)=E(k) O(k),其中E(k)为偶数部分,O(k)为奇数部分。

E(k):

E(k)=0 for n in 0:N/2: E(k) =f(2*n)*W(N,2*n,k)

O(k):

O(k)=0 for n in 0:N/2: O(k) =f(2*n 1)*W(N,2*n 1,k)

只要分别求出E(k),O(k)的值,F(k)=E(k) O(k),F(k N/2)=E(k)-O(k).

其实这里写代码也就可以了,写for循环就能求出相应的值。

但是,这个因为还是有for循环,所以虽然加速了,但是仍然慢,我们要想办法加速,看能不能把for循环消灭掉。

如何求E(k)的值?

将f(n)中偶数部分单独拿出来作为一个信号:fe(m), 则E(k)的公式可改写为:

E(k):

E(k)=0 for m in 0:N/2: E(k) =fe(m)*W(N,2*m,k)

也可以写为:

E(k)=0 for m in 0:N/2: E(k) =fe(m)*W(N/2,m,k)

我们再接着把N/2表示为:M,则也可以写为:

E(k)=0 for m in 0:M: E(k) =f1(m)*W(M,m,k)

注意到这个写法,与我们求F(k)的写法一样:

F(k)=0 for n in range(0:N): F(k) =f(n)*W(N,n,k)

观察E(k)与F(k),发现结构一摸一样。

回想一下,我们求F(k),k=1,2,3,4,....N时,只用求前半部分F(1),F(2),,,,,F(N/2)即可,F(N/2 1)...F(N)的值,都可以由前半部分给出。

对于F(k),将其分解成了F(k) = E(k) O(k),k=1,2,3,4,...N/2

这列E(k)的写法与F(k)一摸一样了,F(k)的性质, E(k)同样具有。

也就是要求E(k),k=1,2,3,...N/2,只用求E(1),E(2),...E(N/4)即可,从E(N/4 1),..到E(N/2)都可以用前半部分来求。

对E(k)做奇偶分解:E(k)=E_E(k) E_O(k)

总的来说,就是将F(k)做奇偶分解得到的E和O ,其中E 又可以做奇偶分解。

如果O能够做奇偶分解的话,就完美了,此时,我们只用实现一个奇偶分解的函数,然后递归调用就能计算出F(k)的值了。

如何求O(k)的值?

观察O:

O(k):

O(k)=0 for n in 0:N/2: O(k) =f(2*n 1)*W(N,2*n 1,k)

f(n)中的奇数部分表示为fo(m), 然后令:M=N/2,得:

O(k):

O(k)=0 for m in 0:M: O(k) =fo(m)*W(2M,2*m 1,k)

而W(2M,2m 1,k) = W(2M,2m,k)*W(2M,1,k)

又因为:W(N,2n,k) = W(N/2,n,k),

所以上式W(2M,2m,k)*W(2M,1,k)=W(2M/2,m,k)*W(2M,1,k)=W(M,m,k)*W(2M,1,k)

️:

O(k):

O(k)=0 for m in 0:M: O(k) =fo(m)*W(M,m,k)*W(2M,1,k)

注意到,在for循环中,m是变化的,所以W(2M,1,k)是个常量,我们可以将其提取出来,for循环结束后再乘,即:

O(k):

O(k)=0 for m in 0:M: O(k) =fo(m)*W(M,m,k) O(k)=O(k)*W(2M,1,k)

我们只看前三行,是不是形式与F(k)求取过程又一样了。这样形式是可以进行奇偶分解的:O(k) = [OE(k) OO(k) ]* W(2M,1,k)

我们总结一下:

F(k)=E(k) O(k)

E(k)=EE(k) EO(k)

O(k)=W(2M,1,k)*[ OE(k) OO(k) ]

也就是说,F(k)进行奇偶分解后得到的E和O,又都可以进行奇偶分解。

显然,EE(k)也是可以进行奇偶分解的,EO(k),OE(k),OO(k)都是可以进行奇偶分解的。

直到M=1时便不能再分解了。

所以,我们的伪代码为:

设输入信号f(n),n=1,2,3,4...N 要求其FFT为F(k),k=1,2,3,...N

则:

def getF[ f(n) ] E,O=奇偶分解[f(n)] F[k]=E O,k=1,2,3,4....N/2 F[k]=E-O,k=N/2 1,N/2 2,...N. return F def 奇偶分解[f(n)]: E = getF( f[0:2:n] ),偶数部分可以用得到F的方法得到,因为他们公式相同, 输入的信号变为f[0:2:n]了,只取偶数部分 O = getF( f[1:2:n] ),偶数部分可以用得到F的方法得到,因为他们公式相同输入的信号变为f[0:2:n]了,只取偶数部分 O=O*W(2M,1,k) retunr E,O

可以看到,我们写了两个函数: getF 和奇偶分解,他们组成递归调用,就可以递归的进行奇偶分解,递归调用的过程中,每进入一次奇偶分解函数,输入到getF 函数的序列长度就减半,当这个序列的长度为2时,就可以不再递归了:

N=len(f(n))

如果 N=2:

那么 E= f(0)*W(N,0,0)

O=f(0)*W(N,0,0)*W(2*M,1,0)

此时M=N/2=1

我们把这段写进奇偶分解函数里:

def 奇偶分解[f(n)]: if N==2: E= f(0)*W(N,0,0) M=N/2 O=f(0)*W(N,0,0)*W(2*M,1,0) else: E = getF( f[0:2:n] ),偶数部分可以用得到F的方法得到,因为他们公式相同, 输入的信号变为f[0:2:n]了,只取偶数部分 O = getF( f[1:2:n] ),偶数部分可以用得到F的方法得到,因为他们公式相同输入的信号变为f[0:2:n]了,只取偶数部分 O=O*W(2M,1,k) retunr E,O

到这里,完成了连个函数getF,和奇偶分解函数,利用这两个函数递归的调用,完成了对奇偶分解后的值递归的进行奇偶分解的作用。

当 N==2时,就不需要再用for循环了,所以,就没有必要再调用getF函数了,直接写出E 和 O值即可。

python实现

奇偶分解函数的实现:

dget函数的批量使用方法,dget函数能否下拉填充(6)

注意到,这里fft_add 及时getF函数

MultiWnk是一个函数,实现了复数W(N,n,k) 与复数 W(N,1,k)乘法的功能。

getF函数的实现:

dget函数的批量使用方法,dget函数能否下拉填充(7)

使用这两个函数完成FFT的蝶形算法,并与scipy的fft进行了对比,代码为:

#coding:utf-8 import numpy as np import matplotlib.pyplot as plt from scipy.fftpack import fft,ifft # input : ys # output: F def Wnk(N,n,k): N=float(N) n=float(n) k=float(k) return np.array([np.cos(2*np.pi/N*n*k),np.sin(2*np.pi/N*n*k)]) #(a bi)(c di)=(ac-bd) (bc ad)i def MultiWnk(wnk1,wnk2): a=wnk1[0] b=wnk1[1] c=wnk2[0] d=wnk2[1] return np.array([a*c-b*d,b*c a*d]) def EO(xs): N=len(xs) if N==2: E = xs[0]*Wnk(N,0,0) #O = xs[1]*Wnk(N,0,0)*Wnk(N,1,0) O = xs[1]*Wnk(N,0 1,0) #O = xs[1]*MultiWnk(Wnk(N,0,0),Wnk(N,1,0)) #O=[O[k,:]*Wnk(N,1,k) for k in range(0,N/2)] #O = O*Wnk(N,1,0) else: E=fft_add(xs[range(0,N,2)]) O=fft_add(xs[range(1,N,2)]) #O=[O[k,:]*Wnk(N,1,k) for k in range(0,N/2)] # 使用复数乘法 O=[MultiWnk(O[k,:],Wnk(N,1,k)) for k in range(0,N/2)] return E,O def EO1(xs): N = len(xs) E=np.zeros((N/2,2)) O=np.zeros((N/2,2)) for k in range(0,N/2): for i in range(0,N/2): E[k,:] =xs[2*i]*Wnk(N,2*i,k) #O[k,:] =xs[2*i 1]*Wnk(N,2*i 1,k) O[k,:] =xs[2*i 1]*MultiWnk(Wnk(N,2*i,k),Wnk(N,1,k)) return E,O def fft_add(xs): N = len(xs) F=np.zeros((N,2)) E,O=EO(xs) F[range(N/2),:]=E O; F[range(N/2,N),:]=E-O; print("F:%s"%(F)) return F def absF(F): Fabs=np.zeros((F.shape[0])) for i in range(len(F)): Fabs[i]=(F[i,0]**2 F[i,1]**2)**0.5 return Fabs if __name__=="__main__": print("hello") xs=np.linspace(0,1,32) print(xs) ys=5.3*np.sin(2*np.pi*5*xs) print(ys) plt.subplot(221);plt.plot(np.arange(len(xs)),ys);plt.title("ys") F=fft_add(ys) print(absF(F)) plt.subplot(222);plt.plot(np.arange(len(xs)),absF(F)/len(xs));plt.title("ys's fft by zmm") plt.subplot(224);plt.plot(np.arange(len(xs)),abs(fft(ys))/len(xs));plt.title("ys' fft by scipy") plt.show()

dget函数的批量使用方法,dget函数能否下拉填充(8)

上一页123下一页

栏目热文

文档排行

本站推荐

Copyright © 2018 - 2021 www.yd166.com., All Rights Reserved.