Usando o Apple FFT e o Accelerate Framework

Alguém já usou o Apple FFT para um aplicativo do iPhone ou sabe onde posso encontrar um exemplo de aplicativo para usá-lo? Eu sei que a Apple tem algum código de exemplo postado, mas não tenho certeza de como implementá-lo em um projeto real.

    Acabei de obter o código da FFT funcionando para um projeto do iPhone:

    • criar um novo projeto
    • exclua todos os arquivos, exceto main.m e xxx_info.plist
    • indo para as configurações do projeto e procurar por pch e pará-lo de tentar carregar um .pch (visto que acabamos de apagá-lo)
    • copie e cole o exemplo de código sobre o que você tem no main.m
    • remova a linha que # inclui Carbono. O carbono é para o OSX.
    • exclua todas as estruturas e adicione a aceleração da estrutura

    Você também pode precisar remover uma input do info.plist que informa ao projeto para carregar um xib, mas tenho 90% de certeza de que você não precisa se preocupar com isso.

    NOTA: Programe as saídas para o console, os resultados saem como 0.000 que não é um erro – é muito rápido

    Este código é realmente estupidamente obscuro; é generosamente comentado, mas os comentários não facilitam a vida.

    Basicamente no coração disso é:

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

    FFT em n flutuações reais e, em seguida, reverter para voltar para onde começamos. ip significa in-loco, o que significa que & A é sobrescrito. Essa é a razão de toda essa malarkey de empacotamento especial – para que possamos esmagar o valor de retorno no mesmo espaço que o valor de envio.

    Para dar alguma perspectiva (como, por exemplo: por que estaríamos usando essa function em primeiro lugar?), Digamos que queremos executar detecção de pitch na input de microfone, e configuramos para que algum callback seja acionado toda vez o microfone entra em 1024 flutuadores. Supondo que a taxa de amostragem do microfone fosse de 44.1kHz, isto é ~ 44 frameworks / seg.

    Portanto, nossa janela de tempo é qualquer que seja a duração de 1024 amostras, ou seja, 1/44 s.

    Então nós empacotamos A com 1024 flutuadores do microfone, setamos log2n = 10 (2 ^ 10 = 1024), pré-calculamos algumas bobinas (setupReal) e:

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

    Agora A conterá n / 2 números complexos. Estes representam bins de frequência n / 2:

    • bin [1] .idealFreq = 44Hz – ie A frequência mais baixa que podemos detectar com segurança é UMA onda completa dentro dessa janela, ou seja, uma onda de 44Hz.

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

    • etc.

    • bin [512] .idealFreq = 512 * 44Hz – A maior freqüência que podemos detectar (conhecida como freqüência de Nyquist) é quando cada par de pontos representa uma onda, ou seja, 512 ondas completas dentro da janela, ou seja, 512 * 44Hz, ou: n / 2 * bin [1] .idealFreq

    • Na verdade, há um Bin extra, Bin [0], que é frequentemente chamado de ‘DC Offset’. Acontece que Bin [0] e Bin [n / 2] sempre terão componente complexo 0, então A [0] .realp é usado para armazenar Bin [0] e A [0] .imagp é usado para armazenar Bin [ n / 2]

    E a magnitude de cada número complexo é a quantidade de energia que vibra em torno dessa frequência.

    Então, como você pode ver, não seria um detector de pitch muito grande, já que não tem granularidade suficiente. Existe um truque astuto Extraindo frequências precisas de FFT Bins usando a mudança de fase entre os frameworks para obter a frequência precisa para um determinado bin.

    Ok, agora no código:

    Observe o ‘ip’ em vDSP_fft_zrip, = ‘in place’, ou seja, a saída sobrescreve A (‘r’ significa que recebe inputs reais)

    Veja a documentação em vDSP_fft_zrip,

    Os dados reais são armazenados na forma complexa dividida, com reais ímpares armazenados no lado imaginário da forma complexa dividida e até reais armazenados no lado real.

    esta é provavelmente a coisa mais difícil de entender. Estamos usando o mesmo contêiner (& A) durante todo o processo. Então, no começo, queremos preenchê-lo com n números reais. depois da FFT, ela vai conter números complexos n / 2. nós então jogamos isso na transformada inversa, e esperamos obter nossos n números reais reais.

    agora a estrutura de A é sua configuração para valores complexos. Então, o vDSP precisa padronizar como compactar números reais nele.

    Então, primeiro nós geramos n números reais: 1, 2, …, n

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

    Em seguida, nós os empacotamos em A como n / 2 #s complexos:

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

    Você realmente precisaria ver como A está alocado para obter isso, talvez procure COMPLEX_SPLIT na documentação.

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

    Em seguida, fazemos um pré-cálculo.


    Aula rápida de DSP para bods de matemática: a teoria de Fourier leva muito tempo para se entender (eu tenho visto isso por vários anos)

    Um cisiso é:

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

    isto é, um ponto no círculo unitário no plano complexo.

    Quando você multiplica números complexos, os ângulos são adicionados. Então z ^ k continuará pulando ao redor do círculo unitário; z ^ k pode ser encontrado em um ângulo k.theta

    • Escolha z1 = 0 + 1i, ou seja, um quarto de volta do eixo real, e observe que z1 ^ 2 z1 ^ 3 z1 ^ 4 dão outro quarto de volta de forma que z1 ^ 4 = 1

    • Escolha z2 = -1, ou seja, meia volta. também z2 ^ 4 = 1, mas z2 completou 2 ciclos neste ponto (z2 ^ 2 também é = 1). Então você poderia pensar em z1 como a frequência fundamental e z2 como o primeiro harmônico

    • Da mesma forma, z3 = o ponto 'três quartos de uma revolução', ou seja, -i completa exatamente 3 ciclos, mas, na verdade, avançando 3/4 de cada vez é o mesmo que retroceder 1/4 de cada vez

    ou seja, z3 é apenas z1, mas na direção oposta - é chamado de aliasing

    z2 é a maior frequência significativa, pois escolhemos 4 amostras para manter uma onda completa.

    • z0 = 1 + 0i, z0 ^ (qualquer coisa) = 1, isso é offset CC

    Você pode expressar qualquer sinal de 4 pontos como uma combinação linear de z0 z1 e z2, ou seja, você está projetando-o nesses vetores de base

    mas eu ouço você perguntando "o que significa projetar um sinal em um cisóide?"

    Você pode pensar desta forma: A agulha gira em torno do cisóide, então na amostra k, a agulha está apontando na direção k.theta, e o comprimento é sinal [k]. Um sinal que corresponda exatamente à frequência do cissóide irá inchar a forma resultante em alguma direção. Então, se você adicionar todas as contribuições, terá um vetor resultante forte. Se a frequência quase coincidir, a protuberância será menor e se moverá lentamente ao redor do círculo. Para um sinal que não corresponde à frequência, as contribuições serão canceladas.

    http://complextoreal.com/tutorials/tutorial-4-fourier-analysis-made-easy-part-1/ irá ajudá-lo a obter uma compreensão intuitiva.

    Mas a essência é; se tivermos escolhido projetar 1024 amostras em {z0, ..., z512}, teríamos o pré-cálculo de z0 a z512, e é isso que essa etapa de pré-cálculo é .


    Observe que, se você estiver fazendo isso em código real, provavelmente desejará fazer isso uma vez quando o aplicativo for carregado e chamar a function de liberação complementar uma vez quando ele for encerrado. Não faça isso muitas vezes - é caro.

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

    Vale a pena notar que, se definirmos log2n como, por exemplo, 8, você pode lançar esses valores pré-calculados em qualquer function fft que use resolução < = 2 ^ 8. Então (a menos que você queira otimização de memória final) basta criar um conjunto para a maior resolução que você vai precisar e usá-lo para tudo.

    Agora as transformações reais, fazendo uso das coisas que acabamos de calcular:

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

    Neste ponto A estará contendo n / 2 números complexos, apenas o primeiro é na verdade dois números reais (offset CC, Nyquist #) mascarados como um número complexo. A visão geral da documentação explica essa embalagem. É bastante interessante - basicamente permite que os resultados (complexos) da transformação sejam empacotados na mesma memory que as inputs (reais, mas estranhamente empacotadas).

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

    e de volta ... ainda precisaremos descompactar nossa matriz original de A. então comparamos apenas para verificar se recuperamos exatamente com o que começamos, liberamos nossas bobinas pré-calculadas e concluímos!

    Mas espere! antes de você descompactar, há uma última coisa que precisa ser feita:

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

    Aqui está um exemplo do mundo real: Um snippet de c ++ que usa as rotinas vftp fft do Accelerate para fazer autocorrelação na input da unidade de áudio Remote IO. Usar esta estrutura é bem complicado, mas a documentação não é tão ruim.

     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; iimNumberBuffers; ++ii) { bzero(ioData->mBuffers[ii].mData, ioData->mBuffers[ii].mDataByteSize); } } 

    Embora eu diga que o FFT Framework da Apple é rápido … Você precisa saber como uma FFT funciona para obter uma detecção de pitch precisa (ou seja, calcular a diferença de fase em cada FFT sucessiva para encontrar o pitch exato, e não o pitch do a maioria domina bin).

    Não sei se é de alguma ajuda, mas enviei meu object Pitch Detector do meu aplicativo de afinação (musicianskit.com/developer.php). Existe também um exemplo do projeto xCode 4 para download (assim você pode ver como a implementação funciona).

    Estou trabalhando para fazer o upload de um exemplo de implementação de FFT. Fique ligado e atualizarei isso assim que isso acontecer.

    Codificação feliz!

    Aqui está outro exemplo do mundo real: https://github.com/krafter/DetectingAudioFrequency