iPhone FFT與Accelerate框架vDSP

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

感謝

Answer 1

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感興趣，則感興趣的分區i由i = N * f/Fs給出。

附加說明：您將需要申請一個合適的window function（例如Hann aka Hanning）到您的FFT輸入數據，之前計算FFT本身。

Answer 2

您可以檢查有關這方面的蘋果文件並妥善保管數據包裝。 這裏是我的例子，你需要小心的是計算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); 

} 

Answer 3

一件事。我將我的結果與fftw庫FFT進行了比較，並且使用vDSP庫計算的變換的虛部在索引0處（即0頻率，即DC）總是有不同的值。 我應用的另一種方法是將實部和虛部都除以2的倍數。我想這是由於函數中使用的算法。而且，這兩個問題都發生在FFT過程中，而不是IFFT過程中。

我用vDSP_fft_zrip。我希望這可以幫助。

Paolo