带有加速框架vDSP的iPhone FFT

我很难用vDSP来实现FFT。 我理解这个理论，但是我正在寻找一个具体的代码示例。

我有一个wav文件的数据如下：

问题1.如何将audio数据放入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. 将audio采样数据放入input的实部，并将虚部归零。

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 。

附加说明：在计算FFT本身之前，您需要将适当的窗口函数 （例如Hann aka Hanning ）应用于FFTinput数据。

Answer 2

你可以检查有关的苹果文件，并妥善保pipe数据包装。 这是我的例子

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

有一点你需要注意的是计算FFT的直stream分量。 我将结果与fftw库FFT进行了比较，并且使用vDSP库计算的变换的虚部总是在索引0（即0频率，即DC）处具有不同的值。 我应用的另一个方法是将实部和虚部都除以2。我想这是由于函数中使用的algorithm。 而且，这两个问题都发生在FFT过程中，而不是在IFFT过程中。

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

保罗