Transformed problem of Fourier

Asked

Viewed 505 times

2

I’m having a problem using the fft function of numpy. 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
  • 1

    I’m not sure, but in function transformada_direta(), the sum divided by N is not the transformation inverse? It shouldn’t just be c[k]=soma? (instead of c[k]=soma/(N))

  • Is not x=np.linspace(0.0,T*(N-1),N) that should be x=np.linspace(0.0, T*N, N)? In the above form, each interval of x will be 2/3 pi, while in the suggested would be 1/2 pi.

1 answer

1

Your FFT is different from Wikipedia, see clicking here.

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

Your sample is too small, which is why you are not being able to properly verify the functioning of your code. The larger the sample, the more faithful and less inaccurate FFT’s 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))

Browser other questions tagged python numpy

You are not signed in. Login or sign up in order to post.