代码之家 › 专栏 › 技术社区 › Clément

使用GNU八度FFT函数

gnuplot octave fft

Clément · 技术社区 · 15 年前

我在玩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));

基本上,它所做的就是取一个函数,计算一个区间的图像,对其进行快速傅立叶变换,然后保留一些谐波,并对结果进行模糊处理。但是我得到了一个垂直压缩的图(输出的垂直比例是错误的)。有什么想法吗?

2 回复 | 直到 15 年前

Olivier Verdier 15 年前

你可能做错了。删除代码中的所有“负”频率。你应该保持正负低频。下面是一个用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))

mtrw 15 年前

你正在扔掉转换的后半部分。对于实值输入,转换是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 是八度音阶的关键词;用你自己的变数来打击它们是一个很好的方法,让真正奇怪的虫子上路!

推荐文章

AmericanMade · 倍频程中的deal函数不工作

7 年前

Li Feiyi · octave软件包在windows上安装不成功

7 年前

Dimitris · Matlab/Octave中所有内置符号的列表

7 年前

daiana05 · 在倍频程中实现Sobel梯度

7 年前

marolafm · 倍频程/Matlab-使用现有绘图对图形进行子绘图

7 年前

sajran · Octave handle类方法返回object而不是number

7 年前

Pavlin · 在3d矩阵中查找值

7 年前

Eduardo Pignatelli · 如何以八度音阶解析json文件?

7 年前

Levi · 在倍频程中绘制判定边界线

7 年前

Rochak Agrawal · 以下八度音阶代码段的矢量化代码是什么

7 年前