使用JTransforms库计算自相关性的FFT

10

我正在尝试使用以下代码计算时间序列中样本窗口的自相关性。 我将FFT应用于该窗口,然后计算实部和虚部的幅度,并将虚部设置为零,最后对其进行逆变换以获得自相关性:

DoubleFFT_1D fft = new DoubleFFT_1D(magCnt);
fft.realForward(magFFT);

magFFT[0] = (magFFT[0] * magFFT[0]);
for (int i = 1; i < (magCnt - (magCnt%2)) / 2; i++) {
    magFFT[2*i] = magFFT[2*i] * magFFT[2*i] + magFFT[2*i + 1] * magFFT[2*i + 1];
    magFFT[2*i + 1] = 0.0;
}

if (magCnt % 2 == 0) {
    magFFT[1] = (magFFT[1] * magFFT[1]);
} else {
    magFFT[magCnt/2] = (magFFT[magCnt-1] * magFFT[magCnt-1] + magFFT[1] * magFFT[1]);
}

autocorr = new double[magCnt];
System.arraycopy(magFFT, 0, autocorr, 0, magCnt);
DoubleFFT_1D ifft = new DoubleFFT_1D(magCnt);
ifft.realInverse(autocorr, false);

for (int i = 1; i < autocorr.length; i++)
    autocorr[i] /= autocorr[0];
autocorr[0] = 1.0;

首先,问题是:这段代码将自相关结果映射到[0,1]范围内,虽然相关性应该在-1和1之间。当然,可以将结果映射到[-1,1]范围内,但我不确定这种映射是否正确。我们如何解释结果中autocorr数组中的值?

其次,使用此代码对某些周期性序列获得良好结果,即根据信号的周期性获得特定自相关指标的较高值。但是,当我将其应用于非周期性信号时,结果会变得奇怪:所有autocorr数组中的值都非常接近于1。这是什么原因?

1个回答

4
要使基于FFT的算法正常工作,您必须仔细研究定义,包括您正在使用的库的约定。您似乎混淆了“信号处理”约定和“统计”约定。还有FFT包装和零填充。
以下是适用于偶数N情况、信号处理约定的代码。已针对Brute Force Wrapped Autocorrelation进行了测试。注释显示如何将其转换为信号处理约定。对于统计ac,需要减去数据的平均值。这只需将FFT的“0Hz”分量清零即可完成。然后,自相关的零元素就是方差,可以通过除以此数量来归一化。结果值将落在-1..1之间,正如您所说。
您的代码似乎正在执行除法操作,但未忽略数据的0 Hz分量。因此,它计算出了某种约定混合体。
import edu.emory.mathcs.jtransforms.fft.DoubleFFT_1D;
import java.util.Arrays;

public class TestFFT {

    void print(String msg, double [] x) {
        System.out.println(msg);
        for (double d : x) System.out.println(d);
    }

    /**
     * This is a "wrapped" signal processing-style autocorrelation. 
     * For "true" autocorrelation, the data must be zero padded.  
     */
    public void bruteForceAutoCorrelation(double [] x, double [] ac) {
        Arrays.fill(ac, 0);
        int n = x.length;
        for (int j = 0; j < n; j++) {
            for (int i = 0; i < n; i++) {
                ac[j] += x[i] * x[(n + i - j) % n];
            }
        }
    }

    private double sqr(double x) {
        return x * x;
    }

    public void fftAutoCorrelation(double [] x, double [] ac) {
        int n = x.length;
        // Assumes n is even.
        DoubleFFT_1D fft = new DoubleFFT_1D(n);
        fft.realForward(x);
        ac[0] = sqr(x[0]);
        // ac[0] = 0;  // For statistical convention, zero out the mean 
        ac[1] = sqr(x[1]);
        for (int i = 2; i < n; i += 2) {
            ac[i] = sqr(x[i]) + sqr(x[i+1]);
            ac[i+1] = 0;
        }
        DoubleFFT_1D ifft = new DoubleFFT_1D(n); 
        ifft.realInverse(ac, true);
        // For statistical convention, normalize by dividing through with variance
        //for (int i = 1; i < n; i++)
        //    ac[i] /= ac[0];
        //ac[0] = 1;
    }

    void test() {
        double [] data = { 1, -81, 2, -15, 8, 2, -9, 0};
        double [] ac1 = new double [data.length];
        double [] ac2 = new double [data.length];
        bruteForceAutoCorrelation(data, ac1);
        fftAutoCorrelation(data, ac2);
        print("bf", ac1);
        print("fft", ac2);
        double err = 0;
        for (int i = 0; i < ac1.length; i++)
            err += sqr(ac1[i] - ac2[i]);
        System.out.println("err = " + err);
    }

    public static void main(String[] args) {
        new TestFFT().test();
    }
}

在得到另一个问题的答案后(http://dsp.stackexchange.com/questions/3337/finding-peaks-in-an-autocorrelation-function),我也有同样的想法。需要消除直流分量,但我不确定如何消除它。所以,你的意思是将第一个FFT系数赋值为零将消除直流分量吗?我现在没有机会尝试它。 - mostar
是的。如果您查看fft的定义,第0个分量是输入数据的平均值。因此,将此分量清零等同于在进行FFT之前从所有数据元素中减去平均值。上面注释掉的代码经过轻微测试。 - Gene
在执行 ifft.realInverse(ac, true); 之前,@Gene 在做类似计算功率谱的事情,但是公式不应该像这样:ac[i] = Math.sqrt(sqr(x[i]) + sqr(x[i+1])) 而不是没有 Math.sqrt - Yuriy Kravets
@YuriyKravets 正如我在答案中所说,这段代码对于通常的带环自相关信号处理定义是正确的。如果不是这样,它就不会产生与暴力计算相同的结果,后者实际上实现了该定义。如果您有不同的定义,当然需要不同的代码。 - Gene

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接