소나 에코를 찾기 위한 상호 상관

Question

나는 내 소리의 반향을 감지하려고 노력하고 있습니다. 짹짹 내 소리에 녹음 Android에서는 교차 상관이 두 신호의 FFT가 유사한 위치를 찾는 가장 적절한 방법인 것 같습니다. 거기에서 거리에 해당하는 교차 상관 배열의 피크를 식별할 수 있습니다.

내 이해를 바탕으로 다음과 같은 상호 상관 함수를 생각해 냈습니다.이 올바른지?처음에 0을 추가하고 몇 가지 요소를 다시 시작할지 잘 모르겠습니다.

public double[] xcorr1(double[] recording, double[] chirp) {        
    double[] recordingZeroPadded = new double[recording.length + chirp.length];

    for (int i = recording.length; i < recording.length + chirp.length; ++i)
            recordingZeroPadded[i] = 0;

    for (int i = 0; i < recording.length; ++i)
            recordingZeroPadded[i] = recording[i];

    double[] result = new double[recording.length + chirp.length - 1];

    for (int offset = 0; offset < recordingZeroPadded.length - chirp.length; ++offset)
        for (int i = 0; i < chirp.length; ++i)
            result[offset] += chirp[i] * recordingZeroPadded[offset + i];
    return result;
}

보조 질문:

에 따르면 이것 대답은 다음과 같이 계산할 수도 있습니다.

corr(a, b) = ifft(fft(a_and_zeros) * fft(b_and_zeros[reversed]))

전혀 이해하지 못하지만 구현하기에는 충분히 쉬운 것 같습니다.즉 나는 실패했다고 말했다(내 가정에서는 xcorr1 맞다).제가 완전히 오해하고 있는 것 같은데요?

public double[] xcorr2(double[] recording, double[] chirp) {
    // assume same length arguments for now
    DoubleFFT_1D fft = new DoubleFFT_1D(recording.length);
    fft.realForward(recording);
    reverse(chirp);
    fft.realForward(chirp);
    double[] result = new double[recording.length];

    for (int i = 0; i < result.length; ++i)
        result [i] = recording[i] * chirp[i];

    fft.realInverse(result, true);
    return result;
}

두 가지가 모두 작동한다고 가정하면 배열에 수천 개의 요소가 포함된다는 점을 고려하면 어떤 함수가 가장 적합할까요?

편집하다:그런데 FFT 버전에서는 두 배열의 양쪽 끝에 0을 추가해 보았습니다.

---

SleuthEye의 응답 후 편집:

내가 '실제' 데이터를 다루고 있기 때문에 실제 변환을 수행하여 계산의 절반(실제 부분)만 수행하면 된다는 것을 확인할 수 있습니까?

코드에서는 REAL 변환에서 반환된 배열의 홀수 요소가 가상인 것처럼 보입니다.여기서 무슨 일이 일어나고 있는 걸까요?

실수 배열에서 복소수 배열로 어떻게 가나요?아니면 이것이 변환의 목적입니까?실수를 복소수 영역으로 옮기려면?(그러나 실수는 복소수의 부분 집합일 뿐이므로 이미 이 영역에 있지 않습니까?)

realForward가 실제로 허수/복소수를 반환하는 경우 complexForward와 어떻게 다릅니까?그리고 결과를 어떻게 해석하나요?복소수의 크기는 무엇입니까?

변환에 대한 이해가 부족한 점에 대해 사과드립니다. 저는 지금까지 푸리에 급수만 연구했습니다.

코드를 보내주셔서 감사합니다.다음은 '내' 작업 구현입니다.

public double[] xcorr2(double[] recording, double[] chirp) {
    // pad to power of 2 for optimisation
    int y = 1;
    while (Math.pow(2,y) < recording.length + chirp.length)
        ++y;
    int paddedLength = (int)Math.pow(2,y);

    double[] paddedRecording = new double[paddedLength];
    double[] paddedChirp = new double[paddedLength];

    for (int i = 0; i < recording.length; ++i)
            paddedRecording[i] = recording[i];

    for (int i = recording.length; i < paddedLength; ++i)
            paddedRecording[i] = 0;

    for (int i = 0; i < chirp.length; ++i)
            paddedChirp[i] = chirp[i];

    for (int i = chirp.length; i < paddedLength; ++i)
            paddedChirp[i] = 0;

    reverse(chirp);
    DoubleFFT_1D fft = new DoubleFFT_1D(paddedLength);
    fft.realForward(paddedRecording);
    fft.realForward(paddedChirp);
    double[] result = new double[paddedLength];

    result[0] = paddedRecording[0] * paddedChirp[0]; // value at f=0Hz is real-valued
    result[1] = paddedRecording[1] * paddedChirp[1]; // value at f=fs/2 is real-valued and packed at index 1
    for (int i = 1; i < result.length / 2; ++i) {
        double a = paddedRecording[2*i];
        double b = paddedRecording[2*i + 1];
        double c = paddedChirp[2*i];
        double d = paddedChirp[2*i + 1];

        // (a+b*j)*(c-d*j) = (a*c+b*d) + (b*c-a*d)*j
        result[2*i]     = a*c + b*d;
        result[2*i + 1] = b*c - a*d;
    }

    fft.realInverse(result, true);

    // discard trailing zeros
    double[] result2 = new double[recording.length + chirp.length - 1];
    for (int i = 0; i < result2.length; ++i)
        result2[i] = result[i];

    return result2;
}

그러나 각각 약 5000개의 요소가 있을 때까지는 xcorr1이 더 빠른 것 같습니다.특별히 느린 작업을 수행하고 있나요(아마도 메모리를 지속적으로 '새로 만들기' - ArrayList로 캐스팅해야 할 수도 있음)?아니면 테스트하기 위해 배열을 생성한 임의의 방식인가요?아니면 반전하는 대신 공액체를 수행해야 합니까?즉, 성능은 실제로 문제가 되지 않으므로 뭔가 분명한 것이 없으면 최적화를 지적할 필요가 없습니다.

Answer 1

귀하의 구현 xcorr1 상호 상관의 표준 신호 처리 정의에 해당합니다.

처음에 0을 추가하는 것과 관련된 질문과 관련하여:첨가 chirp.length-1 0은 결과의 인덱스 0이 전송 시작에 해당하도록 만듭니다.그러나 상관관계 출력의 피크가 발생한다는 점에 유의하십시오. chirp.length-1 에코 시작 후의 샘플(처프는 수신된 전체 에코와 정렬되어야 함)에코 지연을 얻기 위해 피크 인덱스를 사용하면 지연을 빼거나 첫 번째 지연을 삭제하여 해당 상관기 지연을 조정해야 합니다. chirp.length-1 결과를 출력합니다.추가 0은 처음에 많은 추가 출력에 해당하므로 처음부터 해당 0을 추가하지 않는 것이 더 나을 것입니다.

을 위한 xcorr2 그러나 몇 가지 사항을 해결해야 합니다.첫째, 만약 recording 그리고 chirp 입력은 최소한 chirp+recording을 위해 아직 0으로 채워지지 않았습니다. 데이터 길이를 그렇게 해야 합니다(성능상의 이유로 2의 거듭제곱 길이가 바람직함).아시다시피 둘 다 동일한 길이로 패딩되어야 합니다.

둘째, 당신은 다음에 표시된 곱셈을 고려하지 않았습니다. 게시된 참조 답변, 실제로는 복잡한 곱셈에 해당합니다(반면 DoubleFFT_1D.realForward API는 double을 사용합니다).이제 처프의 FFT를 사용하여 복소수 곱셈과 같은 작업을 구현하려는 경우 처프의 FFT의 켤레 복소수를 사용하여 실제로 곱셈을 구현할 수도 있습니다(대체 구현은 다음 그림에 표시됨). 참고 답변), 시간 영역 값을 반대로 할 필요가 없습니다.

또한 회계 DoubleFFT_1D.realForward 짝수 길이 변환을 위한 포장 순서는 다음과 같습니다.

// [...]
fft.realForward(paddedRecording);
fft.realForward(paddedChirp);

result[0] = paddedRecording[0]*paddedChirp[0]; // value at f=0Hz is real-valued
result[1] = paddedRecording[1]*paddedChirp[1]; // value at f=fs/2 is real-valued and packed at index 1
for (int i = 1; i < result.length/2; ++i) {
    double a = paddedRecording[2*i];
    double b = paddedRecording[2*i+1];
    double c = paddedChirp[2*i];
    double d = paddedChirp[2*i+1];

    // (a+b*j)*(c-d*j) = (a*c+b*d) + (b*c-a*d)*j
    result[2*i]   = a*c + b*d;
    result[2*i+1] = b*c - a*d;
}
fft.realInverse(result, true);
// [...]

참고 result 배열의 크기는 다음과 같습니다. paddedRecording 그리고 paddedChirp, 하지만 첫 번째만 recording.length+chirp.length-1 유지되어야합니다.

마지막으로, 수천 개의 요소 배열에 가장 적합한 함수와 관련하여 FFT 버전 xcorr2 훨씬 더 빨라질 가능성이 높습니다(배열 길이를 2의 거듭제곱으로 제한하는 경우).

Answer 2

직접 버전은 처음으로 무늬가 Zero-Padding이 필요하지 않습니다. 길이 M와 길이의 Chirp의 기록을 녹음하고 길이 N의 결과를 계산합니다. 손으로 작은 예를 통해 작업을 수행하기 위해 손으로 일하십시오 :

recording = [1, 2, 3]
chirp = [4, 5]

  1 2 3
4 5

  1 2 3
  4 5

  1 2 3
    4 5

  1 2 3
      4 5


result = [1*5, 1*4 + 2*5, 2*4 + 3*5, 3*4] = [5, 14, 23, 4]

.

FFT 메소드는 긴 배열이있는 경우 훨씬 빠릅니다. 이 경우 M + N-1을 크기 M + N-1로 제로 패딩하여 두 입력 어레이가 모두 FFT를 사용하기 전에 과 동일한 크기입니다.

또한 FFT 출력은 복잡한 숫자이므로 복잡한 곱셈 를 사용해야합니다. ...에 (1 + 2J) * (3 + 4J)는 3 + 8J가 아닌 -5 + 10J입니다. 복잡한 숫자가 어떻게 정렬되거나 처리되는지 모르겠지만 이것이 맞는지 확인하십시오.

또는 이것은 이것이 변형의 목적입니다. 실수를 복잡한 도메인으로 이동하려면?

아니오, 푸리에 변환은 시간 영역에서 주파수 영역으로 변환합니다. 시간 영역 데이터는 실제 또는 복잡 할 수 있으며 주파수 도메인 데이터는 실제 또는 복잡 할 수 있습니다. 대부분의 경우 복잡한 스펙트럼이있는 실제 데이터가 있습니다. 푸리에 변환을 읽어야합니다.

실제로 상상의 / 복소수를 반환하는 것은 실제로, 복잡한 것과 어떻게 다른가요?

실제 FFT는 진짜 입력 을 취하고 복합 FFT는 복잡한 입력 을 취합니다. 두 변환 모두 복소수를 출력으로 생성합니다. 그것이 DFT가하는 일입니다. DFT가 실제 출력을 생성하는 유일한 시간은 입력 데이터가 대칭이면 (DCT를 사용하여 더 많은 시간을 절약 할 수 있음)