Saltar al contenido

cómo extraer la frecuencia asociada con los valores fft en python

Solución:

np.fft.fftfreq le dice las frecuencias asociadas con los coeficientes:

import numpy as np

x = np.array([1,2,1,0,1,2,1,0])
w = np.fft.fft(x)
freqs = np.fft.fftfreq(len(x))

for coef,freq in zip(w,freqs):
    if coef:
        print('{c:>6} * exp(2 pi i t * {f})'.format(c=coef,f=freq))

# (8+0j) * exp(2 pi i t * 0.0)
#    -4j * exp(2 pi i t * 0.25)
#     4j * exp(2 pi i t * -0.25)

El OP pregunta cómo encontrar la frecuencia en Hertz. Creo que la fórmula es frequency (Hz) = abs(fft_freq * frame_rate).

Aquí hay un código que lo demuestra.

Primero, creamos un archivo de onda a 440 Hz:

import math
import wave
import struct

if __name__ == '__main__':
    # http://stackoverflow.com/questions/3637350/how-to-write-stereo-wav-files-in-python
    # http://www.sonicspot.com/guide/wavefiles.html
    freq = 440.0
    data_size = 40000
    fname = "test.wav"
    frate = 11025.0
    amp = 64000.0
    nchannels = 1
    sampwidth = 2
    framerate = int(frate)
    nframes = data_size
    comptype = "NONE"
    compname = "not compressed"
    data = [math.sin(2 * math.pi * freq * (x / frate))
            for x in range(data_size)]
    wav_file = wave.open(fname, 'w')
    wav_file.setparams(
        (nchannels, sampwidth, framerate, nframes, comptype, compname))
    for v in data:
        wav_file.writeframes(struct.pack('h', int(v * amp / 2)))
    wav_file.close()

Esto crea el archivo test.wav. Ahora leemos en los datos, FFT, encontramos el coeficiente con la potencia máxima y encontramos la frecuencia fft correspondiente, y luego convertimos a Hertz:

import wave
import struct
import numpy as np

if __name__ == '__main__':
    data_size = 40000
    fname = "test.wav"
    frate = 11025.0
    wav_file = wave.open(fname, 'r')
    data = wav_file.readframes(data_size)
    wav_file.close()
    data = struct.unpack('{n}h'.format(n=data_size), data)
    data = np.array(data)

    w = np.fft.fft(data)
    freqs = np.fft.fftfreq(len(w))
    print(freqs.min(), freqs.max())
    # (-0.5, 0.499975)

    # Find the peak in the coefficients
    idx = np.argmax(np.abs(w))
    freq = freqs[idx]
    freq_in_hertz = abs(freq * frate)
    print(freq_in_hertz)
    # 439.8975

Frecuencias asociadas con los valores DFT (en python)

Por fft, Transformada Rápida de Fourier, entendemos a un miembro de una gran familia de algoritmos que permiten la rápido cálculo de la DFT, Transformada de Fourier discreta, de una señal equimuestreada.

Una DFT convierte una lista de norte números complejos a una lista de norte números complejos, en el entendido de que ambas listas son periódicas con período norte.

Aquí nos ocupamos de la numpy implementacion de fft.

En muchos casos, piensas en

  • una señal X definido en el dominio del tiempo de la longitud norte, muestreado en un intervalo constante dt,
  • su DFT X (aquí específicamente X = np.fft.fft(x)), cuyos elementos se muestrean en el eje de frecuencia con una frecuencia de muestreo dw.

Alguna definición

  • el período (también conocido como duración) de la señal x, muestreado en dt con N muestras es es

    T = dt*N
    
  • las frecuencias fundamentales (en Hz y en rad / s) de X, su DFT es

    df = 1/T
    dw = 2*pi/T # =df*2*pi
    
  • la frecuencia superior es la frecuencia de Nyquist

    ny = dw*N/2
    

    (y no es dw*N)

Las frecuencias asociadas con un elemento particular en la DFT.

Las frecuencias correspondientes a los elementos en X = np.fft.fft(x) para un índice dado 0<=n<N se puede calcular de la siguiente manera:

def rad_on_s(n, N, dw):
    return dw*n if n<N/2 else dw*(n-N)

o en un solo barrido

w = np.array([dw*n if n<N/2 else dw*(n-N) for n in range(N)])

si prefiere considerar las frecuencias en Hz, s/w/f/

f = np.array([df*n if n<N/2 else df*(n-N) for n in range(N)])

Usando esas frecuencias

Si desea modificar la señal original x -> y aplicando un operador en el dominio de la frecuencia en la forma de una función de frecuencia solamente, el camino a seguir es calcular el w‘arena

Y = X*f(w)
y = ifft(Y)

Introduciendo np.fft.fftfreq

Por supuesto numpy tiene una función de conveniencia np.fft.fftfreq que vuelve frecuencias adimensionales en vez de los dimensionales pero es tan fácil como

f = np.fft.fftfreq(N)*N*df
w = np.fft.fftfreq(N)*N*dw

Porque df = 1/T y T = N/sps (sps siendo el número de muestras por segundo) también se puede escribir

f = np.fft.fftfreq(N)*sps

La frecuencia es solo el índice de la matriz. En el índice norte, la frecuencia es 2πn / la longitud de la matriz (radianes por unidad). Considerar:

>>> numpy.fft.fft([1,2,1,0,1,2,1,0])
array([ 8.+0.j,  0.+0.j,  0.-4.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+4.j,
        0.+0.j])

el resultado tiene valores distintos de cero en los índices 0, 2 y 6. Hay 8 elementos. Esto significa

       2πit/8 × 0       2πit/8 × 2       2πit/8 × 6
    8 e           - 4i e           + 4i e
y ~ ———————————————————————————————————————————————
                          8
¡Haz clic para puntuar esta entrada!
(Votos: 0 Promedio: 0)



Utiliza Nuestro Buscador

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *