使用加速框架 vDSP 进行 iPhone FFT

Question

我正在尝试使用vDSP实现FFT，但遇到了困难。我理解其原理，但需要一个具体的代码示例。

我有来自wav文件的数据如下：

问题1：如何将音频数据放入FFT中？

问题2：如何从FFT中获取输出数据？

问题3：最终目标是检测低频声音。我应该如何做到这一点？

-(OSStatus)open:(CFURLRef)inputURL{
OSStatus result = -1;

result = AudioFileOpenURL (inputURL, kAudioFileReadPermission, 0, &mAudioFile);
if (result == noErr) {
    //get  format info
    UInt32 size = sizeof(mASBD);

    result = AudioFileGetProperty(mAudioFile, kAudioFilePropertyDataFormat, &size, &mASBD);

    UInt32 dataSize = sizeof packetCount;
    result = AudioFileGetProperty(mAudioFile, kAudioFilePropertyAudioDataPacketCount, &dataSize, &packetCount);
    NSLog([NSString stringWithFormat:@"File Opened, packet Count: %d", packetCount]);

    UInt32 packetsRead = packetCount;
    UInt32 numBytesRead = -1;
    if (packetCount > 0) { 
        //allocate  buffer
        audioData = (SInt16*)malloc( 2 *packetCount);
        //read the packets
        result = AudioFileReadPackets (mAudioFile, false, &numBytesRead, NULL, 0, &packetsRead,  audioData); 
        NSLog([NSString stringWithFormat:@"Read %d  bytes,  %d packets", numBytesRead, packetsRead]);
    }
}
return result;
}

以下是FFT代码：

log2n = N;
n = 1 << log2n;
stride = 1;
nOver2 = n / 2;

printf("1D real FFT of length log2 ( %d ) = %d\n\n", n, log2n);

/* Allocate memory for the input operands and check its availability,
 * use the vector version to get 16-byte alignment. */

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

if (originalReal == NULL || A.realp == NULL || A.imagp == NULL) {
printf("\nmalloc failed to allocate memory for  the real FFT"
"section of the sample.\n");
exit(0);
}

/* Generate an input signal in the real domain. */
for (i = 0; i < n; i++)

    originalReal[i] = (float) (i + 1);

/* Look at the real signal as an interleaved complex vector  by
 * casting it.  Then call the transformation function vDSP_ctoz to
 * get a split complex vector, which for a real signal, divides into
 * an even-odd configuration. */

vDSP_ctoz((COMPLEX *) originalReal, 2, &A, 1, nOver2);

/* Set up the required memory for the FFT routines and check  its
 * availability. */

setupReal = vDSP_create_fftsetup(log2n, FFT_RADIX2);

if (setupReal == NULL) {

printf("\nFFT_Setup failed to allocate enough memory  for"
"the real FFT.\n");

exit(0);
}

/* Carry out a Forward and Inverse FFT transform. */
vDSP_fft_zrip(setupReal, &A, stride, log2n, FFT_FORWARD);
vDSP_fft_zrip(setupReal, &A, stride, log2n, FFT_INVERSE);

/* Verify correctness of the results, but first scale it by  2n. */
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);

/* The output signal is now in a split real form.  Use the  function
 * vDSP_ztoc to get a split real vector. */
vDSP_ztoc(&A, 1, (COMPLEX *) obtainedReal, 2, nOver2);

/* Check for accuracy by looking at the inverse transform  results. */
Compare(originalReal, obtainedReal, n);

感谢

- Simon

---

如果你只是想检测低频声音，那么使用FFT可能会过度。你要寻找什么特定的频率，以及需要多少分辨率？ - Paul R

我正在寻找包含鼓或低音声音的任何频率，以便我可以响应节奏。谢谢。 - Simon

在这种情况下，你最好使用低通滤波器+包络检测器——实现起来更简单，而且由于计算成本要少得多，所以对电池寿命应该更有利。 - Paul R

Answer 1

1. 将音频样本数据放入输入的实部，虚部清零。

2. 如果你只对每个频域中的幅度感兴趣，则可以计算每个输出频率的sqrt(re*re + im*im)。如果你只对相对幅度感兴趣，则可以省略sqrt并计算平方幅度，(re*re + im*im)。

3. 你需要查看与你感兴趣的频率对应的频率或频率（参见2）。如果你的采样率为Fs，FFT大小为N，则输出bin i对应的频率为f = i * Fs / N。反之，如果你对特定频率f感兴趣，则感兴趣的bin i由i = N * f / Fs给出。

附注：在计算FFT本身之前，您需要对FFT输入数据应用适当的窗口函数（例如汉宁窗）。

- Paul R

---

你能举个例子，在使用vDSP方法进行fft之前应用窗口函数的方法吗？ - some_id

Answer 2

你可以查看苹果的文档并注意数据打包。以下是我的示例：

//  main.cpp
//  FFTTest
//
//  Created by Harry-Chris Stamatopoulos on 11/23/12.
//  

/* 
 This is an example of a hilbert transformer using 
 Apple's VDSP fft/ifft & other VDSP calls.
 Output signal has a PI/2 phase shift.
 COMPLEX_SPLIT vector "B" was used to cross-check
 real and imaginary parts coherence with the original vector "A"
 that is obtained straight from the fft.
 Tested and working. 
 Cheers!
*/

#include <iostream>
#include <Accelerate/Accelerate.h>
#define PI 3.14159265
#define DEBUG_PRINT 1

int main(int argc, const char * argv[])
{
    float fs = 44100;           //sample rate
    float f0 = 440;             //sine frequency
    uint32_t i = 0;

    uint32_t L = 1024;

    /* vector allocations*/
    float *input = new float [L];
    float *output = new float[L];
    float *mag = new float[L/2];
    float *phase = new float[L/2];

    for (i = 0 ; i < L; i++)
    {
        input[i] = cos(2*PI*f0*i/fs);
    }

    uint32_t log2n = log2f((float)L);
    uint32_t n = 1 << log2n;
    //printf("FFT LENGTH = %lu\n", n);

    FFTSetup fftSetup;
    COMPLEX_SPLIT A;
    COMPLEX_SPLIT B;
    A.realp = (float*) malloc(sizeof(float) * L/2);
    A.imagp = (float*) malloc(sizeof(float) * L/2);

    B.realp = (float*) malloc(sizeof(float) * L/2);
    B.imagp = (float*) malloc(sizeof(float) * L/2);

    fftSetup = vDSP_create_fftsetup(log2n, FFT_RADIX2);

    /* Carry out a Forward and Inverse FFT transform. */
    vDSP_ctoz((COMPLEX *) input, 2, &A, 1, L/2);
    vDSP_fft_zrip(fftSetup, &A, 1, log2n, FFT_FORWARD);

    mag[0] = sqrtf(A.realp[0]*A.realp[0]);

    //get phase
    vDSP_zvphas (&A, 1, phase, 1, L/2);
    phase[0] = 0;

    //get magnitude;
    for(i = 1; i < L/2; i++){
        mag[i] = sqrtf(A.realp[i]*A.realp[i] + A.imagp[i] * A.imagp[i]);
    }

    //after done with possible phase and mag processing re-pack the vectors in VDSP format
    B.realp[0] = mag[0];
    B.imagp[0] = mag[L/2 - 1];;

    //unwrap, process & re-wrap phase
    for(i = 1; i < L/2; i++){
        phase[i] -= 2*PI*i * fs/L;
        phase[i] -= PI / 2 ;
        phase[i] += 2*PI*i * fs/L;
    }

    //construct real & imaginary part of the output packed vector (input to ifft)
    for(i = 1; i < L/2; i++){
        B.realp[i] = mag[i] * cosf(phase[i]);
        B.imagp[i] = mag[i] * sinf(phase[i]);
    }

#if DEBUG_PRINT
    for (i = 0 ; i < L/2; i++)
    {
       printf("A REAL = %f \t A IMAG = %f \n", A.realp[i], A.imagp[i]);
       printf("B REAL = %f \t B IMAG = %f \n", B.realp[i], B.imagp[i]);
    }
#endif
    //ifft
    vDSP_fft_zrip(fftSetup, &B, 1, log2n, FFT_INVERSE);

    //scale factor
    float scale = (float) 1.0 / (2*L);

    //scale values
    vDSP_vsmul(B.realp, 1, &scale, B.realp, 1, L/2);
    vDSP_vsmul(B.imagp, 1, &scale, B.imagp, 1, L/2);

    //unpack B to real interleaved output
    vDSP_ztoc(&B, 1, (COMPLEX *) output, 2, L/2);

    // print output signal values to console
    printf("Shifted signal x = \n");
    for (i = 0 ; i < L/2; i++)
        printf("%f\n", output[i]);

    //release resources
    free(input);
    free(output);
    free(A.realp);
    free(A.imagp);
    free(B.imagp);
    free(B.realp);
    free(mag);
    free(phase);
}

- Harry Chris Stamatopoulos

Answer 3

需要注意的一点是计算FFT时的直流分量。我将我的结果与fftw库FFT进行了比较，使用vDSP库计算的变换虚部在索引0处（即0频率，也就是直流）总是有不同的值。

另一个措施是将实部和虚部都除以2。我猜测这是由于函数中使用的算法导致的。此外，这两个问题都发生在FFT过程中而非IFFT过程中。

我使用了vDSP_fft_zrip。

- Paolo Turati