2011-06-15 93 views
8

我很难实现使用vDSP的FFT。我理解这个理论,但是我正在寻找一个具体的代码示例。iPhone FFT与Accelerate框架vDSP

我有如下一个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); 

感谢

+0

如果您只想检测低频声音,那么使用FFT可能会矫枉过正。你在寻找什么特定的频率/频率,以及多少分辨率? – 2011-06-15 14:05:18

+0

我正在寻找任何包含鼓或贝司声音的频率,以便我可以响应节拍​​。谢谢 – Simon 2011-06-15 15:16:37

+2

在这种情况下,使用低通滤波器+包络检波器可能会更好 - 实施起来更简单,电池寿命应该更容易,因为它的计算成本更低。 – 2011-06-15 15:19:06

回答

9
  1. 你把你的音频采样数据到输入的实部和零虚部。

  2. 如果您只是对频域中每个bin的大小感兴趣,那么您可以为每个输出bin计算sqrt(re*re + im*im)。如果您只对相对幅度感兴趣,那么您可以删除sqrt并只计算平方幅度(re*re + im*im)

  3. 您可以查看与您感兴趣的频率或频率相对应的箱或箱的大小(请参阅(2))。如果您的采样率是Fs,并且您的FFT大小是N,则输出仓i的相应频率由f = i * Fs/N给出。相反,如果您对特定频率f感兴趣,则感兴趣的分区ii = N * f/Fs给出。

附加说明:您将需要申请一个合适的window function(例如Hann aka Hanning)到您的FFT输入数据,之前计算FFT本身。

+2

您可以举一个使用vDSP方法在fft之前应用窗口函数的例子吗? – jarryd 2013-05-31 10:45:18

5

您可以检查有关这方面的苹果文件并妥善保管数据包装。 这里是我的例子,你需要小心的是计算FFT的直流分量

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

} 
0

一件事。我将我的结果与fftw库FFT进行了比较,并且使用vDSP库计算的变换的虚部在索引0处(即0频率,即DC)总是有不同的值。 我应用的另一种方法是将实部和虚部都除以2的倍数。我想这是由于函数中使用的算法。而且,这两个问题都发生在FFT过程中,而不是IFFT过程中。

我用vDSP_fft_zrip。我希望这可以帮助。

Paolo