programing

Apple FFT 및 Accelerate Framework 사용

sourcejob 2023. 8. 12. 10:15
반응형

Apple FFT 및 Accelerate Framework 사용

누군가 사용한 적이 있습니까?Apple FFT아이폰 앱의 경우, 사용 방법에 대한 샘플 앱을 어디서 찾을 수 있는지 알고 계십니까?애플이 샘플 코드를 게시한 것은 알고 있지만 실제 프로젝트에 어떻게 구현해야 할지 잘 모르겠습니다.

방금 아이폰 프로젝트에서 FFT 코드를 받았습니다.

  • 새 프로젝트 생성
  • main.m 및 xxx_info.plist를 제외한 모든 파일을 삭제합니다.
  • 프로젝트 설정으로 이동하여 pch를 검색하고 .pch를 로드하지 않도록 합니다(우리가 방금 삭제한 것처럼 보기).
  • 코드 예제를 메인에 있는 모든 것 위에 복사하여 붙여넣습니다.
  • #의 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);

실제 플로트에서 FFT를 실행한 다음 역방향으로 전환하여 시작 위치로 돌아갑니다.ip는 in-place의 약자인데, 이것은 &A가 덮어쓴다는 것을 의미합니다. 이것이 바로 이 모든 특수 포장 말라키의 이유입니다. 그래서 우리는 반환 값을 송신 값과 동일한 공간에 밀어 넣을 수 있습니다.

예를 들어, 처음부터 이 기능을 사용하는 이유는 무엇일까요?와 같은 관점을 제공하기 위해 마이크 입력에 대해 피치 감지를 수행하고 마이크가 1024 플로트에 들어갈 때마다 콜백이 트리거되도록 설정했다고 가정합니다.마이크 샘플링 속도가 44.1kHz라고 가정하면 초당 최대 44프레임입니다.

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

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

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

이제 A에는 n/2개의 복소수가 포함됩니다.다음은 n/2 주파수 빈을 나타냅니다.

  • bin[1]idealFreq = 44Hz -- 즉, 신뢰성 있게 감지할 수 있는 가장 낮은 주파수는 해당 윈도우 내의 하나의 전체 파형, 즉 44Hz 파형입니다.

  • bin[2].이상적인 Freq = 2 * 44Hz

  • 기타.

  • bin[512].idealFreq = 512 * 44Hz -- 탐지할 수 있는 가장 높은 주파수(Nyquist 주파수라고 함)는 모든 포인트 쌍이 창 내의 512 전체 파형을 나타내는 곳입니다. 즉, 512 * 44Hz 또는: n/2 * bin[1].idealFreq

  • 실제로 흔히 'DC 오프셋'이라고 하는 추가적인 Bin[0]이 있습니다.Bin[0]과 Bin[n/2]은 항상 복합 성분 0을 갖게 되므로 A[0].realp는 Bin[0]과 A[0]을 저장하는 데 사용됩니다.imagp는 Bin[n/2]을 저장하는 데 사용됩니다.

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

보시다시피 거의 미세한 입자를 가지고 있지 않기 때문에 그다지 좋은 피치 검출기는 아닙니다.프레임 간 위상 변화를 사용하여 FFT Bin에서 정확한 주파수를 추출하여 주어진 Bin에 대한 정확한 주파수를 얻는 교묘한 속임수가 있습니다.

자, 이제 코드로 넘어가겠습니다.

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

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

실제 데이터는 분할 복합 형식으로 저장되며, 홀수 실제 데이터는 분할 복합 형식의 가상 측에 저장되고 실제 데이터도 실제 측에 저장됩니다.

이것이 아마도 가장 이해하기 어려운 것일 것입니다.우리는 공정 내내 같은 컨테이너(&A)를 사용하고 있습니다.그래서 처음에는 n개의 실수로 채우고자 합니다. FFT 이후에는 n/2개의 복소수를 보유하게 됩니다.그런 다음 우리는 그것을 역변환에 던지고, 바라건대 우리의 원래 n개의 실수를 꺼낼 수 있습니다.

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

그래서 우리는 먼저 실수를 생성합니다: 1, 2, ..., n.

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

다음으로 A에 n/2 complex #s로 포장합니다.

// 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));

다음에는 사전 계산을 합니다.


수학 본문에 대한 빠른 DSP 클래스: 푸리에 이론은 머리를 돌리는 데 오랜 시간이 걸립니다 (저는 몇 년 동안 그것을 계속 보고 있었습니다).

시소이드는 다음과 같습니다.

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

즉, 복소 평면에서 단위 원의 점.

복소수를 곱하면, 각도가 더해집니다.그래서 z^k는 단위 원 주위를 계속해서 깡충깡충 뛸 것입니다; z^k는 각도 k에서 찾을 수 있습니다.

  • z1 = 0+1i, 즉 실제 축에서 4분의 1 회전을 선택하고, z1^2 z1^3 z1^4가 각각 다른 4분의 1 회전을 주어 z1^4 = 1이 되도록 합니다.

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

  • 마찬가지로 z3 = '회전의 4분의 3' 지점, 즉 -i는 정확히 3 사이클을 완료하지만, 실제로는 매 번 3/4 앞으로 이동하는 것은 매 번 1/4 뒤로 이동하는 것과 같습니다.

즉, z3는 단지 z1이지만 반대 방향입니다. 그것은 에일리어스라고 불립니다.

z2는 최대 파형을 유지하기 위해 4개의 샘플을 선택했기 때문에 가장 의미 있는 주파수입니다.

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

임의의 4점 신호를 z0 z1과 z2의 선형 조합으로 표현할 수 있습니다. 즉, 이 기본 벡터에 투영하는 것입니다.

하지만 "시소이드에 신호를 투사한다는 것이 무슨 뜻입니까?"라고 묻는 것을 들었습니다.

다음과 같이 생각할 수 있습니다.바늘은 시소이드 주위를 회전하므로 샘플 k에서 바늘은 k.theta 방향을 가리키고 있으며 길이는 신호 [k]입니다.시소이드의 주파수와 정확히 일치하는 신호는 결과적인 모양을 어떤 방향으로 벌지시킬 것입니다.따라서 모든 기여를 합하면 강력한 결과 벡터를 얻을 수 있습니다.주파수가 거의 일치하면 팽대부보다 더 작고 원 주위에서 천천히 움직입니다.주파수와 일치하지 않는 신호의 경우 기여는 서로 취소됩니다.

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을 eg8로 설정하면 분해능 <= 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);

다음은 실제 예입니다. Accelerate의 vDSP fft 루틴을 사용하여 원격 IO 오디오 장치의 입력에 대해 자동 상관을 수행하는 c++의 일부입니다.이 프레임워크를 사용하는 것은 꽤 복잡하지만 문서화는 그리 나쁘지 않습니다.

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);
    }
}

애플의 FFT 프레임워크는 빠르지만,정확한 피치 검출(즉, 가장 지배적인 빈의 피치가 아닌 정확한 피치를 찾기 위해 각 연속 FFT에서 위상 차이를 계산)을 얻으려면 FFT의 작동 방식을 알아야 합니다.

도움이 될지 모르겠지만 튜너 앱(musicianskit.com/developer.php) 에서 Pitch Detector 개체를 업로드했습니다.다운로드를 위한 xCode 4 프로젝트의 예도 있습니다. 구현이 어떻게 작동하는지 확인할 수 있습니다.

FFT 구현 예제를 업로드하는 중입니다. 잠시 기다려 주십시오. 업데이트가 되면 업데이트하겠습니다.

해피 코딩!

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

반응형