Utilisation des fonctions FFT de GNU Octave

Question

Je joue avec les fonctions fft d'Octave, et je n'arrive pas vraiment à trouver comment mettre à l'échelle leur sortie : J'utilise le code suivant (très court) pour approximer une fonction :

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));

En gros, ce qu'il fait, c'est prendre une fonction, calculer l'image d'un intervalle, faire un fft, puis garder quelques harmoniques, et faire un ifft du résultat. Pourtant, j'obtiens un tracé qui est compressé verticalement (l'échelle verticale de la sortie est fausse). Avez-vous une idée ?

Answer 1

Vous jetez la seconde moitié de la transformation. La transformation est symétrique hermitienne pour les entrées à valeur réelle et vous devez conserver ces lignes. Essaie ça :

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 transformée inverse aura invariablement une minuscule partie imaginaire due à l'erreur de calcul numérique, d'où l'utilisation de l'option real extraction.

Notez que input y size sont des mots-clés dans Octave ; les englober avec vos propres variables est un bon moyen d'obtenir des bogues vraiment bizarres par la suite !

Answer 2

Vous vous y prenez probablement mal. Vous supprimez toutes les fréquences "négatives" dans votre code. Vous devriez garder les basses fréquences positives et négatives. Voici un code en python et le résultat. Le graphique a la bonne échelle.

texte alternatif http://files.droplr.com/files/35740123/XUl90.fft.png

Le code :

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))