| DSP |
Autokorrelation
Moin,
zur Bestimmung der Grundfrequenz eines Audiosignals benötige ich die Autokorrelation. Diese kann ja bekanntermaßen über mehrere Wege berechnet werden: einmal über die Faltung des originalen Signals mit sich selbst um [tau] verschoben und außerdem über die inverse Fouriertransformation des Autoleistungsspektrums (Wiener-Khinchin-Theorem). Ich habe mal beide Methoden implementiert (aus Effizienzgründen soll in der finalen Anwendung letztere genutzt werden), leider unterscheiden sich jedoch die Ergebnisse. 
Zunächst einmal der Code der ersten Methode:
| code: |
1:
2:
3:
4:
5:
6:
7:
8:
9:
|
static float* calculateAutocorrelationUsingBruteForce(float *m_samples) {
float *autocorr = (float *)malloc(sizeof(float)*N);
for (int tau = 0; tau < N; tau++) {
float sum = 0;
for (int t = 0; t < N - tau; t++) sum += m_samples[t] * m_samples[t + tau];
autocorr[tau] = sum;
}
return autocorr;
} |
|
Eigentlich relativ straightforward. 
Nun zur zweiten Methode: Für die FFTs und Vektoroperationen habe ich Apples Accelerate Framework genutzt, da dieses besonders effizient und ziemlich schnell rechnet. Hier der Code:
| code: |
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
|
static float* calculateAutocorrelationUsingFFT(float *m_samples) {
// calculate autocorrelation according to the Wiener–Khinchin theorem
// prepare calculation
FFTSetup m_fftSetup = vDSP_create_fftsetup(LOG_2N, kFFTRadix2);
DSPSplitComplex m_splitComplex;
m_splitComplex.realp = (float *)malloc(sizeof(float)*N/2);
m_splitComplex.imagp = (float *)malloc(sizeof(float)*N/2);
vDSP_ctoz((DSPComplex *)m_samples, 2, &m_splitComplex, 1, N/2);
// perform FFT
vDSP_fft_zrip(m_fftSetup, &m_splitComplex, 1, LOG_2N, kFFTDirection_Forward);
// calculate the absolute squares respectively the squared magnitudes
vDSP_zvmags(&m_splitComplex, 1, m_splitComplex.realp, 1, N/2);
memset(m_splitComplex.imagp, 0, sizeof(float)*N/2);
// perform IFFT
vDSP_fft_zrip(m_fftSetup, &m_splitComplex, 1, LOG_2N, kFFTDirection_Inverse);
// store autocorrelation result vector
vDSP_ztoc(&m_splitComplex, 1, (DSPComplex *)m_samples, 2, N/2);
return m_samples;
} |
|
Leider kann ich in keiner der angeführten Funktionen einen Fehler feststellen, dennoch sind die Ergebnisse grundverschieden. Kann mir da vielleicht jemand von Euch weiterhelfen? 
Ich bin für jede Hilfe sehr dankbar!!
Liebe Grüße,
Alex |
