Fourier transform problem

Question

Fourier transform problem

Navigation

#1 by (1 votes)

2

I'm having a problem using numpy's fft function. When rotating it, the answer is giving different from the Fourier transform that I implemented manually and by the tests I did the manual seems correct. I think I'm using fft the wrong way.

import numpy as np
import math as m

def transformada_direta(F,x,N):
    c=[]
    for k in range (N):
        c.append(0)
    for k in range (N):
        soma=0
        for j in range (N):
            soma=soma+(F[j]*(np.exp((0-1j)*k*x[j])))
        c[k]=soma/(N)
    return c
def main():
    N=4
    T=(2*(np.pi))/N
    x=np.linspace(0.0,T*(N-1),N)
    F=[]
    for j in range (N):
        F.append(0)
    F[0]=5
    F[1]=-1
    F[2]=3
    F[3]=1
    c=transformada_direta(F,x,N)
    print(c)
    yf=np.fft.fft(F)
    print(yf)
main()

python numpy

asked by anonymous 23.06.2018 / 04:33

1 answer

Error "Transaction is required to perform this operation" when inserting, editing or deleting - Wildfly Save checked checkboxes and load them when the page is reloaded

score 1 · Answer 1

Your FFT is different from Wikipedia, see here .

Below, I rewrote the C ++ algorithm, in the wiki above, for Python. Make a comparison.

Your sample is too small, so you are not able to properly check your code. The larger the sample, the more faithful and less inaccurate the FFT response will be.

Note that, by Nyquist's theorem, the FFT will only have valid response for frequencies below 32 Hz, in the example below. Therefore, values above 32 Hz should be disregarded.

import numpy
import math

def separate(X):
    return X[::2] + X[1::2]

def fft(F):
    N = len(F)

    if N >= 2:
        S = separate(F)
        FirstHalf = fft(S[:N/2])
        SecondHalf = fft(S[N/2:])
        S = FirstHalf + SecondHalf
        for i in range(N/2):
            e = S[i]
            o = S[N/2+i]
            w = numpy.exp((0-2j)*numpy.pi*i)
            S[i] = e+ o*w
            S[N/2+i] = e- o*w
        return S
    else:
        return F

nsamples = 64
timespan = 1.0
samplerate = nsamples / timespan
freqresolution = samplerate / nsamples
freqs = [2,5,11,17,29]
X = [ sum(map(lambda f : math.sin(2*numpy.pi*t/nsamples*f), freqs)) 
      for t in range(nsamples) ]
F = fft(X)
maxF = max(map(lambda f: numpy.abs(f), F))

print("Smp\t  Measured\t       FFT\t Histogram\tFRs")
for i in range(nsamples):
    print("{0:3d}\t{1:10.2f}\t{2:10.2f}\t{3:10}\t{4:3.0f}"
        .format(i, X[i], numpy.abs(F[i]), '*'*(int(numpy.abs(F[i])/maxF*10.0)), i*freqresolution))