問題描述
我正在嘗試使用以下代碼計算時間序列中示例窗口的自相關性.我將 FFT 應用于該窗口,然后計算實部和虛部的大小并將虛部設置為零,最后對其進行逆變換以獲得自相關:
I'm trying to calculate autocorrelation of sample windows in a time series using the code below. I'm applying FFT to that window, then computing magnitudes of real and imaginary parts and setting imaginary part to zero, lastly taking inverse transform of it to obtain autocorrelation:
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
數組中的值?
The first question is: It is seen that this code maps autocorrelation result to [0,1]
range, although correlation is supposed to be between -1 and 1. Of course it is easy to map the results to [-1,1]
range, but I'm not sure if this mapping is correct. How can we interpret the values in the resulting autocorr
array?
其次,使用此代碼,我在某些周期性序列中得到了很好的結果,也就是說,我根據信號的周期為特定的自相關指數獲得了更高的值.但是,當我將其應用于非周期性信號時,結果變得很奇怪:autocorr
數組中的所有值似乎都非常接近 1.這是什么原因?
Secondly, with this code I'm geting good results for some periodic series, that is I get higher values for specific autocorrelation indices in accordance with the period of signal. However the result go weird when I apply it to non-periodic signals: all the values in autocorr
array appear to be very close to 1. What is the reason for that?
推薦答案
要使基于 FFT 的算法正常工作,您必須特別注意定義,包括您正在使用的庫的約定.您似乎混淆了 AC 的信號處理"約定和統計"約定.然后是 FFT 換行和零填充.
For FFT-based algorithms to work, you must pay careful attention to the definitions, including conventions of the library you are using. You seem to be confusing the "signal processing" convention for AC and the "statistical" one. And then there is FFT wrapping and zero padding.
這是一個適用于偶數 N 情況的代碼,即信號處理約定.它針對蠻力包裝的自相關進行了測試.注釋顯示了如何將其轉換為信號處理約定.對于統計 ac,減去數據的平均值.這可以通過將 FFT 的0Hz"分量歸零來完成.那么 ac 的第零個元素就是方差,你可以通過除以這個量來歸一化.如您所說,結果值將落在-1..1.
Here is a code that's working for the even N case, signal processing convention. It's tested against a brute force wrapped autocorrelation. The comments show how to convert it to the signal processing convention. For statistical ac, the mean of the data is subtracted out. This can be done merely by zeroing out the "0Hz" component of the FFT. Then the zero'th element of the ac is the variance, and you can normalize by dividing through by this quantity. The resulting values will fall in -1..1 as you say.
您的代碼似乎在進行除法,但并未忽略數據的 0 Hz 分量.所以它正在計算某種約定的混搭.
Your code seems to be doing the dividing through, but not ignoring the 0 Hz component of the data. So it's computing some kind of mashup of the conventions.
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();
}
}
這篇關于使用 JTransforms 庫通過 FFT 計算自相關的文章就介紹到這了,希望我們推薦的答案對大家有所幫助,也希望大家多多支持html5模板網!