Advertisement
Guest User

Untitled

a guest
Jul 7th, 2015
745
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # modulo fft
  2. from math import
  3.  
  4. def fft(m, f=[]): #genera fft
  5.     result = []  
  6.         n=2<<(m-1) 
  7.     nesp=0.5
  8.     for k in range(1,m+1):
  9.         d=n/(2<<(k-1)
  10.         nesp=nesp*2
  11.         s=0
  12.         for i in range(1,nesp+1):
  13.             r=s
  14.             z1=r>>(m-k)
  15.             z=0
  16.             for q in range(k):
  17.                 z=(((z1>>q) & 01) | (z<<1)
  18.            
  19.             w=2*pi*z/(2<<(k-1))
  20.             e=cos(w)-sin(w)j
  21.             for l in range (d):
  22.                 r=l+s
  23.                 r1=r+d
  24.                 t=e*f[r1]  
  25.                 f[r1]=f[r]-t
  26.                 f[r]=f[r]+t
  27.             s=s+2*d
  28.     #bitReversal(n, m, f=[])       
  29.     for x in range(1, n-1):
  30.         y=0
  31.         for i in range(m):
  32.             y=(((x>>i) & 01) | (y<<1))
  33.         if x<y :
  34.             t=f[x]
  35.             f[x]=f[y]
  36.             f[y]=t
  37.     #Ampiezze fasi 
  38.     nmezzi=n/2 
  39.     for i in range(n):
  40.         f[i]=f[i]/nmezzi
  41.     A[0]=f[0].real
  42.     phi[0]=0
  43.     for k in range(1, nmezzi+1):
  44.         A[k]=sqrt(f[k].real*f[k].real+f[k].imag*f[k].imag])
  45.         result.append(A[k])
  46.     if f[k].real == 0:
  47.         if f[k].imag>0:
  48.             phi[k]=pi
  49.         if f[k].imag==0:
  50.             phi[k]=0
  51.         if f[k].imag<0:
  52.             phi[k]=2*pi-pi/2
  53.     if f[k].real>0:
  54.         if f[k]>=0:
  55.             phi[k]=atan(f[k].imag, f[k].real)
  56.         if f[k]<0:
  57.             phi[k]=2*pi + atan(f[k].imag, f[k].real)
  58.     if f[k].real<0:
  59.         phi[k]=atan(f[k].imag, f[k].real)
  60.      return result
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement