如何解决具有Accelerate框架vDSP的iPhone FFT
| 我在使用vDSP实现FFT时遇到困难。我了解理论,但正在寻找特定的代码示例。 我有来自wav文件的数据,如下所示: 问题1.如何将音频数据放入FFT? 问题2.如何从FFT中获取输出数据? 问题3.最终目标是检查低频声音。我该怎么做?-(Osstatus)open:(CFURLRef)inputURL{
Osstatus result = -1;
result = AudioFileOpenURL (inputURL,kAudioFileReadPermission,&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,&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,stride,log2n,FFT_FORWARD);
vDSP_fft_zrip(setupReal,FFT_INVERSE);
/* Verify correctness of the results,but first scale it by 2n. */
scale = (float) 1.0 / (2 * n);
vDSP_vsmul(A.realp,&scale,A.realp,nOver2);
vDSP_vsmul(A.imagp,A.imagp,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,(COMPLEX *) obtainedReal,nOver2);
/* Check for accuracy by looking at the inverse transform results. */
Compare(originalReal,obtainedReal,n);
谢谢
解决方法
您将音频样本数据放入输入的实部,并将虚部置零。
如果您只对频域中每个bin的大小感兴趣,则可以为每个输出bin计算
sqrt(re*re + im*im)
。如果您只对相对幅度感兴趣,则可以降低sqrt并仅计算平方幅度(re*re + im*im)
。
您将查看与您感兴趣的一个或多个频率相对应的一个或多个bin的大小(请参阅(2))。如果采样率为Fs,FFT大小为N,则输出仓bin4ѭ的对应频率由f = i * Fs / N
给出。相反,如果您对特定频率f感兴趣,则感兴趣的区间i
由i = N * f / Fs
给出。
补充说明:在计算FFT本身之前,您需要对FFT输入数据应用合适的窗口函数(例如Hann aka Hanning)。
, 您可以查看Apple的文档并妥善保管数据打包。
这是我的示例:
// 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,log2n,FFT_FORWARD);
mag[0] = sqrtf(A.realp[0]*A.realp[0]);
//get phase
vDSP_zvphas (&A,phase,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,FFT_INVERSE);
//scale factor
float scale = (float) 1.0 / (2*L);
//scale values
vDSP_vsmul(B.realp,&scale,B.realp,L/2);
vDSP_vsmul(B.imagp,B.imagp,L/2);
//unpack B to real interleaved output
vDSP_ztoc(&B,(COMPLEX *) output,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);
}
, 您需要注意的一件事是计算出的FFT的直流分量。我将结果与fftw库FFT进行了比较,并且使用vDSP库计算的变换的虚部始终在索引0处具有不同的值(这意味着频率为0,因此为DC)。
我采用的另一项措施是将实部和虚部均除以2。我猜这是由于函数中使用的算法所致。同样,这两个问题都在FFT过程中发生,但在IFFT过程中没有发生。
我使用了vDSP_fft_zrip。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。