2010-05-08 25 views
6

estoy jugando con las funciones fft de octava, y no puedo encontrar la manera de escalar su salida: Yo uso el siguiente código (muy corto) para aproximar una función:El uso de GNU Octave FFT funciones

function y = f(x) 
    y = x .^ 2; 
endfunction; 

X=[-4096:4095]/64; 
Y = f(X); 
# plot(X, Y); 

F = fft(Y); 
S = [0:2047]/2048; 

function points = approximate(input, count) 
    size = size(input)(2); 
    fourier = [fft(input)(1:count) zeros(1, size-count)]; 
    points = ifft(fourier); 
endfunction; 

Y = f(X); plot(X, Y, X, approximate(Y, 10)); 

Básicamente, lo que hace es tomar una función, calcular la imagen de un intervalo, fft-it, luego mantener unos pocos armónicos, y obtener el resultado. Sin embargo, obtengo un gráfico que está verticalmente comprimido (la escala vertical de la salida es incorrecta). ¿Algunas ideas?

Respuesta

3

Probablemente lo estés haciendo mal. Eliminas todas las frecuencias "negativas" en tu código. Debes mantener bajas frecuencias tanto positivas como negativas. Aquí hay un código en python y el resultado. La trama tiene la escala correcta.

alt text http://files.droplr.com/files/35740123/XUl90.fft.png

El código:

from __future__ import division 

from scipy.signal import fft, ifft 
import numpy as np 

def approximate(signal, cutoff): 
    fourier = fft(signal) 
    size = len(signal) 
    # remove all frequencies except ground + offset positive, and offset negative: 
    fourier[1+cutoff:-cutoff] = 0 
    return ifft(fourier) 

def quad(x): 
    return x**2 

from pylab import plot 

X = np.arange(-4096,4096)/64 
Y = quad(X) 

plot(X,Y) 
plot(X,approximate(Y,3)) 
+0

Olivier, Rock :) Eso es exactamente lo que necesitaba, ¡Gracias! –

+0

Aunque, ¿qué hace el uso del subíndice de corte negativo? –

+0

CFP, ¡te alegra que te guste! '-cutoff' significa el índice de" corte al último ", es decir,' -1' significa el último índice. Por lo tanto, el corte '[size/2, -cutoff]' significa dejar todo a partir de la mitad, excepto el 'cutoff' last. Una forma más ordenada habría sido: 'fourier [cutoff + 1: -cutoff] = 0'. –

4

Estás tirando la segunda mitad de la transformación. La transformación es hermitiana simétrica para entradas con valores reales y tienes que mantener esas líneas. Prueba esto:

function points = approximate(inp, count) 
    fourier = fft(inp); 
    fourier((count+1):(length(fourier)-count+1)) = 0; 
    points = real(ifft(fourier)); %# max(imag(ifft(fourier))) should be around eps(real(...)) 
endfunction; 

La transformada inversa invariablemente tener alguna pequeña parte imaginaria debido al error de cálculo numérico, por lo tanto, la extracción real.

Tenga en cuenta que input y size son palabras clave en Octave; castigarlos con sus propias variables es una buena manera de obtener errores realmente extraños en el futuro.

+0

Muy bien, gracias! Lo entiendo ahora. ¿Conoces buenas fuentes de documentación sobre fft? –