Apple FFT 및 Accelerate Framework 사용

Apple FFT 및 Accelerate Framework 사용

---

아무도 Apple FFTiPhone 앱에를 사용 했거나 사용 방법에 대한 샘플 응용 프로그램을 어디서 찾을 수 있는지 알고 있습니까? Apple에 샘플 코드가 게시되어 있다는 것을 알고 있지만 실제 프로젝트에 구현하는 방법을 잘 모르겠습니다.

---

방금 iPhone 프로젝트에서 작동하는 FFT 코드가 있습니다.

• 새 프로젝트 생성
• main.m 및 xxx_info.plist를 제외한 모든 파일 삭제
• 프로젝트 설정으로 이동하여 pch를 검색하고 .pch로드 시도를 중지합니다 (방금 삭제 한 것처럼 보임).
• 복사하여 main.m에있는 코드 예제를 붙여 넣으십시오.
• #include의 Carbon 줄을 제거하십시오. Carbon은 OSX 용입니다.
• 모든 프레임 워크를 삭제하고 가속 프레임 워크를 추가합니다.

프로젝트에 xib를로드하라고 알려주는 info.plist에서 항목을 제거해야 할 수도 있지만 90 %는 신경 쓸 필요가 없다고 확신합니다.

참고 : 프로그램 출력은 콘솔로, 결과는 오류가 아닌 0.000으로 나옵니다 – 매우 빠릅니다.

이 코드는 정말 어리석게 모호합니다. 관대하게 주석을 달았지만 실제로 삶을 더 쉽게 만들지는 않습니다.

기본적으로 핵심은 다음과 같습니다.

vDSP_fft_zrip(setupReal, &A, stride, log2n, FFT_FORWARD);
vDSP_fft_zrip(setupReal, &A, stride, log2n, FFT_INVERSE);

n 개의 실제 수레에 FFT를 적용한 다음 역순으로 시작하여 다시 시작합니다. ip는 in-place를 의미합니다. 이는 & A가 덮어 쓰여짐을 의미합니다. 이것이이 모든 특별한 패킹 malarkey에 대한 이유입니다. 그래서 우리는 반환 값을 send 값과 같은 공간에 스쿼시 할 수 있습니다.

몇 가지 관점을 제공하기 위해 (예 : 처음에이 기능을 사용하는 이유는 무엇입니까?) 마이크 입력에서 피치 감지를 수행하고 매번 콜백이 트리거되도록 설정했다고 가정 해 보겠습니다. 마이크는 1024 개의 플로트에 들어갑니다. 마이크 샘플링 속도가 44.1kHz라고 가정하면 ~ 44 프레임 / 초입니다.

따라서 우리의 시간 창은 1024 개 샘플의 시간 지속 시간, 즉 1/44 초입니다.

따라서 마이크에서 1024 개의 플로트를 A에 넣고 log2n = 10 (2 ^ 10 = 1024)을 설정하고 일부 보빈 (setupReal)을 미리 계산하고 다음을 수행합니다.

vDSP_fft_zrip(setupReal, &A, stride, log2n, FFT_FORWARD);

이제 A는 n / 2 개의 복소수를 포함합니다. 이는 n / 2 개의 주파수 빈을 나타냅니다.

• bin [1] .idealFreq = 44Hz-즉 우리가 안정적으로 감지 할 수있는 가장 낮은 주파수는 해당 창 내에서 하나의 완전한 파, 즉 44Hz 파입니다.

• bin [2] .idealFreq = 2 * 44Hz

• 기타

• bin [512] .idealFreq = 512 * 44Hz-감지 할 수있는 가장 높은 주파수 (나이 퀴 스트 주파수라고 함)는 모든 점 쌍이 파동을 나타내는 곳입니다. 즉, 창 내 512 개의 완전한 파동, 즉 512 * 44Hz 또는 다음과 같습니다. n / 2 * bin [1] .idealFreq

• 실제로 'DC 오프셋'이라고하는 추가 Bin, Bin [0]이 있습니다. 따라서 Bin [0] 및 Bin [n / 2]는 항상 복잡한 구성 요소 0을 가지므로 A [0] .realp는 Bin [0]을 저장하는 데 사용되고 A [0] .imagp는 Bin [을 저장하는 데 사용됩니다. n / 2]

그리고 각 복소수의 크기는 그 주파수 주변에서 진동하는 에너지의 양입니다.

따라서 보시다시피 세분성이 거의 없기 때문에 훌륭한 피치 감지기는 아닙니다. 주어진 빈에 대한 정확한 주파수를 얻기 위해 프레임 간 위상 변화를 사용하여 FFT 빈에서 정확한 주파수를 추출 하는 교활한 트릭이 있습니다 .

좋아, 이제 코드에 :

vDSP_fft_zrip의 'ip', = 'in place'즉, 출력이 A를 덮어 씁니다 ( 'r'는 실제 입력이 필요함을 의미 함)

vDSP_fft_zrip에 대한 설명서를 참조하십시오.

실수 데이터는 분할 복합 형식으로 저장되며, 홀수 실수는 분할 복합 형식의 허수쪽에 저장되고 실수는 실수쪽에 저장됩니다.

이것은 아마도 가장 이해하기 어려운 것입니다. 프로세스 내내 동일한 컨테이너 (& A)를 사용하고 있습니다. 그래서 처음에는 n 개의 실수로 채워야합니다. FFT 후에는 n / 2 개의 복소수를 보유하게됩니다. 그런 다음 그것을 역변환에 던져서 원래 n 개의 실수를 얻길 바랍니다.

이제 A의 구조는 복잡한 값에 대한 설정입니다. 따라서 vDSP는 실수를 압축하는 방법을 표준화해야합니다.

그래서 먼저 n 개의 실수를 생성합니다 : 1, 2, ..., n

for (i = 0; i < n; i++)
    originalReal[i] = (float) (i + 1);

다음으로 우리는 그것들을 n / 2 복잡한 #으로 A에 압축합니다.

// 1. masquerades n real #s as n/2 complex #s = {1+2i, 3+4i, ...}
// 2. splits to 
//   A.realP = {1,3,...} (n/2 elts)
//   A.compP = {2,4,...} (n/2 elts)
//
vDSP_ctoz(
          (COMPLEX *) originalReal, 
          2,                            // stride 2, as each complex # is 2 floats
          &A, 
          1,                            // stride 1 in A.realP & .compP
          nOver2);                      // n/2 elts

이것을 얻기 위해 A가 어떻게 할당되는지 살펴 봐야 할 것입니다. 문서에서 COMPLEX_SPLIT를 찾아보십시오.

A.realp = (float *) malloc(nOver2 * sizeof(float));
A.imagp = (float *) malloc(nOver2 * sizeof(float));

다음으로 사전 계산을합니다.

---

수학 bods를위한 빠른 DSP 수업 : 푸리에 이론은 머리를 돌리는 데 오랜 시간이 걸립니다 (지금까지 몇 년 동안 켜고 끌었습니다).

cisoid는 다음과 같습니다.

z = exp(i.theta) = cos(theta) + i.sin(theta)

즉, 복잡한 평면에서 단위 원의 한 점입니다.

복소수를 곱하면 각도가 더해집니다. 따라서 z ^ k는 단위 원 주위를 계속 홉핑합니다. z ^ k는 각도 k.theta에서 찾을 수 있습니다.

• z1 = 0 + 1i, 즉 실제 축에서 1/4 회전을 선택하고 z1 ^ 2 z1 ^ 3 z1 ^ 4가 각각 다른 1/4 회전을 제공하므로 z1 ^ 4 = 1이됩니다.

• z2 = -1, 즉 반 회전을 선택합니다. 또한 z2 ^ 4 = 1이지만 z2는이 시점에서 2주기를 완료했습니다 (z2 ^ 2도 = 1 임). 따라서 z1을 기본 주파수로, z2를 첫 번째 고조파로 생각할 수 있습니다.

• 유사하게 z3 = '3/4 회전'지점, 즉 -i는 정확히 3주기를 완료하지만 실제로 매번 3/4 씩 앞으로가는 것은 매번 1/4로 뒤로 이동하는 것과 같습니다.

즉, z3은 z1이지만 반대 방향입니다. 앨리어싱이라고합니다.

z2는 전체 파동을 유지하기 위해 4 개의 샘플을 선택 했으므로 가장 의미있는 주파수입니다.

• z0 = 1 + 0i, z0 ^ (anything) = 1, 이것은 DC 오프셋입니다.

모든 4 점 신호를 z0 z1과 z2의 선형 조합으로 표현할 수 있습니다. 즉, 이러한 기저 벡터에 투영합니다.

하지만 "시소 이드에 신호를 투사하는 것이 무엇을 의미합니까?"라고 묻는 것을 들었습니다.

당신은 이것을 이렇게 생각할 수 있습니다 : 바늘은 cisoid 주위를 회전합니다. 그래서 샘플 k에서 바늘은 k.theta 방향을 가리키고 있고 길이는 signal [k]입니다. cisoid의 주파수와 정확히 일치하는 신호는 결과 모양을 어떤 방향으로 팽창시킵니다. 따라서 모든 기여를 더하면 강력한 결과 벡터를 얻을 수 있습니다. 주파수가 거의 일치하면 벌지보다 작아지고 원 주위에서 천천히 움직입니다. 주파수와 일치하지 않는 신호의 경우 기여는 서로를 취소합니다.

http://complextoreal.com/tutorials/tutorial-4-fourier-analysis-made-easy-part-1/ 은 직관적 인 이해를 돕습니다.

그러나 요점은 다음과 같습니다. 1024 개의 샘플을 {z0, ..., z512}에 투영하도록 선택했다면 z0부터 z512까지 사전 계산을했을 것 입니다 . 이것이 바로이 사전 계산 단계입니다 .

---

실제 코드에서이 작업을 수행하는 경우 앱이로드 될 때이 작업을 한 번 수행하고 종료 할 때 보완 릴리스 함수를 한 번 호출하고 싶을 것입니다. 여러 번하지 마십시오. 비용이 많이 듭니다.

// let's say log2n = 8, so n=2^8=256 samples, or 'harmonics' or 'terms'
// if we pre-calculate the 256th roots of unity (of which there are 256) 
// that will save us time later.
//
// Note that this call creates an array which will need to be released 
// later to avoid leaking
setupReal = vDSP_create_fftsetup(log2n, FFT_RADIX2);

log2n을 예를 들어 8로 설정하면 미리 계산 된 값을 해상도 <= 2 ^ 8을 사용하는 모든 fft 함수에 던질 수 있다는 점에 주목할 가치가 있습니다. 따라서 (궁극적 인 메모리 최적화를 원하지 않는 한) 필요한 가장 높은 해상도를 위해 하나의 세트를 만들고 모든 것에 사용하십시오.

이제 우리가 방금 미리 계산 한 내용을 사용하는 실제 변환 :

vDSP_fft_zrip(setupReal, &A, stride, log2n, FFT_FORWARD);

이 시점에서 A는 n / 2 개의 복소수를 포함하고 첫 번째 것은 실제로 복소수로 가장하는 두 개의 실수 (DC 오프셋, Nyquist #)입니다. 문서 개요는이 패킹을 설명합니다. 매우 깔끔합니다. 기본적으로 변환의 (복잡한) 결과를 (실제이지만 이상하게 패키지 된) 입력과 동일한 메모리 공간에 압축 할 수 있습니다.

vDSP_fft_zrip(setupReal, &A, stride, log2n, FFT_INVERSE);

그리고 다시 ... 우리는 여전히 A에서 원래 배열을 풀어야 할 것입니다. 그리고 우리는 우리가 시작했던 것을 정확히 되찾았는지 확인하기 위해 비교하고 미리 계산 된 보빈을 해제하고 완료했습니다!

하지만 기다려! 포장을 풀기 전에 마지막으로해야 할 일이 있습니다.

// Need to see the documentation for this one...
// in order to optimise, different routines return values 
// that need to be scaled by different amounts in order to 
// be correct as per the math
// In this case...
scale = (float) 1.0 / (2 * n);

vDSP_vsmul(A.realp, 1, &scale, A.realp, 1, nOver2);
vDSP_vsmul(A.imagp, 1, &scale, A.imagp, 1, nOver2);

---

Here's a real-world example: A snippet of c++ that uses Accelerate's vDSP fft routines to do autocorrelation on the Remote IO audio unit's input. Using this framework is pretty complicated, but the documentation isn't too bad.

OSStatus DSPCore::initialize (double _sampleRate, uint16_t _bufferSize) {
    sampleRate = _sampleRate;
    bufferSize = _bufferSize;
    peakIndex = 0;
    frequency = 0.f;
    uint32_t maxFrames = getMaxFramesPerSlice();
    displayData = (float*)malloc(maxFrames*sizeof(float));
    bzero(displayData, maxFrames*sizeof(float));
    log2n = log2f(maxFrames);
    n = 1 << log2n;
    assert(n == maxFrames);
    nOver2 = maxFrames/2;
    A.realp = (float*)malloc(nOver2 * sizeof(float));
    A.imagp = (float*)malloc(nOver2 * sizeof(float));
    FFTSetup fftSetup = vDSP_create_fftsetup(log2n, FFT_RADIX2);

    return noErr;
}

void DSPCore::Render(uint32_t numFrames, AudioBufferList *ioData) {

    bufferSize = numFrames;
    float ln = log2f(numFrames);

    //vDSP autocorrelation

    //convert real input to even-odd
    vDSP_ctoz((COMPLEX*)ioData->mBuffers[0].mData, 2, &A, 1, numFrames/2);
    memset(ioData->mBuffers[0].mData, 0, ioData->mBuffers[0].mDataByteSize);
    //fft
    vDSP_fft_zrip(fftSetup, &A, 1, ln, FFT_FORWARD);

    // Absolute square (equivalent to mag^2)
    vDSP_zvmags(&A, 1, A.realp, 1, numFrames/2);
    bzero(A.imagp, (numFrames/2) * sizeof(float));    

    // Inverse FFT
    vDSP_fft_zrip(fftSetup, &A, 1, ln, FFT_INVERSE);

    //convert complex split to real
    vDSP_ztoc(&A, 1, (COMPLEX*)displayData, 2, numFrames/2);

    // Normalize
    float scale = 1.f/displayData[0];
    vDSP_vsmul(displayData, 1, &scale, displayData, 1, numFrames);

    // Naive peak-pick: find the first local maximum
    peakIndex = 0;
    for (size_t ii=1; ii < numFrames-1; ++ii) {
        if ((displayData[ii] > displayData[ii-1]) && (displayData[ii] > displayData[ii+1])) {
            peakIndex = ii;
            break;
        }
    }

    // Calculate frequency
    frequency = sampleRate / peakIndex + quadInterpolate(&displayData[peakIndex-1]);

    bufferSize = numFrames;

    for (int ii=0; ii<ioData->mNumberBuffers; ++ii) {
        bzero(ioData->mBuffers[ii].mData, ioData->mBuffers[ii].mDataByteSize);
    }
}

---

While I will say Apple's FFT Framework is fast... You need to know how an FFT works in order to get accurate pitch detection (i.e. calculating the phase difference on each successive FFT in order to find the exact pitch, not the pitch of the most dominate bin).

I don't know if it's of any help, but I uploaded my Pitch Detector object from my tuner app (musicianskit.com/developer.php). There is an example xCode 4 project for download also (so you can see how the implementation works).

I'm working on uploading an example FFT implementation -- so stay tuned and I'll update this once that happens.

Happy coding!

---

Here is an another real-world example: https://github.com/krafter/DetectingAudioFrequency

참고URL : https://stackoverflow.com/questions/3398753/using-the-apple-fft-and-accelerate-framework