使用 GNU Octave FFT 函数

Question

我正在使用 octave 的 fft 函数，但我真的不知道如何缩放它们的输出：我使用以下（非常短的）代码来近似函数：

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

基本上，它所做的就是获取一个函数，计算一个区间的图像，fft-it，然后保留一些谐波，然后对结果进行 ifft。然而，我得到了一个垂直压缩的图（输出的垂直比例是错误的）。有什么想法吗？

Answer 1

您正在丢弃变换的后半部分。对于实值输入，变换是 Hermitian 对称的，您必须保留这些线。试试这个：

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;

由于数值计算错误，逆变换总是会有一些微小的虚部，因此real提取。

注意input 和size 是Octave 中的关键字；用你自己的变量来破坏它们是一个很好的方法来解决非常奇怪的错误！

Answer 2

你可能做错了。您删除了代码中的所有“负”频率。您应该同时保留正负低频。这是python中的代码和结果。该情节具有正确的比例。

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

代码：

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