Jamoma API  0.6.0.a19
CASpectralProcessor.cpp
1 /* Copyright 2007 Apple Inc. All Rights Reserved.
2 
3  Disclaimer: IMPORTANT: This Apple software is supplied to you by
4  Apple Inc. ("Apple") in consideration of your agreement to the
5  following terms, and your use, installation, modification or
6  redistribution of this Apple software constitutes acceptance of these
7  terms. If you do not agree with these terms, please do not use,
8  install, modify or redistribute this Apple software.
9 
10  In consideration of your agreement to abide by the following terms, and
11  subject to these terms, Apple grants you a personal, non-exclusive
12  license, under Apple's copyrights in this original Apple software (the
13  "Apple Software"), to use, reproduce, modify and redistribute the Apple
14  Software, with or without modifications, in source and/or binary forms;
15  provided that if you redistribute the Apple Software in its entirety and
16  without modifications, you must retain this notice and the following
17  text and disclaimers in all such redistributions of the Apple Software.
18  Neither the name, trademarks, service marks or logos of Apple Inc.
19  may be used to endorse or promote products derived from the Apple
20  Software without specific prior written permission from Apple. Except
21  as expressly stated in this notice, no other rights or licenses, express
22  or implied, are granted by Apple herein, including but not limited to
23  any patent rights that may be infringed by your derivative works or by
24  other works in which the Apple Software may be incorporated.
25 
26  The Apple Software is provided by Apple on an "AS IS" basis. APPLE
27  MAKES NO WARRANTIES, EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION
28  THE IMPLIED WARRANTIES OF NON-INFRINGEMENT, MERCHANTABILITY AND FITNESS
29  FOR A PARTICULAR PURPOSE, REGARDING THE APPLE SOFTWARE OR ITS USE AND
30  OPERATION ALONE OR IN COMBINATION WITH YOUR PRODUCTS.
31 
32  IN NO EVENT SHALL APPLE BE LIABLE FOR ANY SPECIAL, INDIRECT, INCIDENTAL
33  OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
34  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
35  INTERRUPTION) ARISING IN ANY WAY OUT OF THE USE, REPRODUCTION,
36  MODIFICATION AND/OR DISTRIBUTION OF THE APPLE SOFTWARE, HOWEVER CAUSED
37  AND WHETHER UNDER THEORY OF CONTRACT, TORT (INCLUDING NEGLIGENCE),
38  STRICT LIABILITY OR OTHERWISE, EVEN IF APPLE HAS BEEN ADVISED OF THE
39  POSSIBILITY OF SUCH DAMAGE.
40 */
41 
42 //#include "AudioFormulas.h"
43 #include "CASpectralProcessor.h"
44 #include "CABitOperations.h"
45 
46 
47 #include <vecLib/vectorOps.h>
48 
49 
50 #define OFFSETOF(class, field)((size_t)&((class*)0)->field)
51 
52 CASpectralProcessor::CASpectralProcessor(UInt32 inFFTSize, UInt32 inHopSize, UInt32 inNumChannels, UInt32 inMaxFrames)
53  : mFFTSize(inFFTSize), mHopSize(inHopSize), mNumChannels(inNumChannels), mMaxFrames(inMaxFrames),
54  mLog2FFTSize(Log2Ceil(mFFTSize)),
55  mFFTMask(mFFTSize - 1),
56  mFFTByteSize(mFFTSize * sizeof(Float32)),
57  mIOBufSize(NextPowerOfTwo(mFFTSize + mMaxFrames)),
58  mIOMask(mIOBufSize - 1),
59  mInputSize(0),
60  mInputPos(0), mOutputPos(-mFFTSize & mIOMask),
61  mInFFTPos(0), mOutFFTPos(0),
62  mSpectralFunction(0), mUserData(0)
63 {
64  mWindow.alloc(mFFTSize, false);
65  SineWindow(); // set default window.
66 
67  mChannels.alloc(mNumChannels);
68  mSpectralBufferList.allocBytes(OFFSETOF(SpectralBufferList, mDSPSplitComplex[mNumChannels]), true);
69  mSpectralBufferList->mNumberSpectra = mNumChannels;
70  for (UInt32 i = 0; i < mNumChannels; ++i)
71  {
72  mChannels[i].mInputBuf.alloc(mIOBufSize, true);
73  mChannels[i].mOutputBuf.alloc(mIOBufSize, true);
74  mChannels[i].mFFTBuf.alloc(mFFTSize, true);
75  mChannels[i].mSplitFFTBuf.alloc(mFFTSize, true);
76  mSpectralBufferList->mDSPSplitComplex[i].realp = mChannels[i].mSplitFFTBuf();
77  mSpectralBufferList->mDSPSplitComplex[i].imagp = mChannels[i].mSplitFFTBuf() + (mFFTSize >> 1);
78  }
79 
80  mFFTSetup = vDSP_create_fftsetup (mLog2FFTSize, FFT_RADIX2);
81 
82 }
83 
84 CASpectralProcessor::~CASpectralProcessor()
85 {
86  mWindow.free();
87  mChannels.free();
88  mSpectralBufferList.free();
89  vDSP_destroy_fftsetup(mFFTSetup);
90 }
91 
92 void CASpectralProcessor::Reset()
93 {
94  mInputPos = 0;
95  mOutputPos = -mFFTSize & mIOMask;
96  mInFFTPos = 0;
97  mOutFFTPos = 0;
98 
99  for (UInt32 i = 0; i < mNumChannels; ++i)
100  {
101  memset(mChannels[i].mInputBuf(), 0, mIOBufSize * sizeof(Float32));
102  memset(mChannels[i].mOutputBuf(), 0, mIOBufSize * sizeof(Float32));
103  memset(mChannels[i].mFFTBuf(), 0, mFFTSize * sizeof(Float32));
104  }
105 }
106 
107 const double two_pi = 2. * M_PI;
108 
109 void CASpectralProcessor::HanningWindow()
110 {
111  // this is also vector optimized
112 
113  double w = two_pi / (double)(mFFTSize - 1);
114  for (UInt32 i = 0; i < mFFTSize; ++i)
115  {
116  mWindow[i] = (0.5 - 0.5 * cos(w * (double)i));
117  }
118 }
119 
120 void CASpectralProcessor::SineWindow()
121 {
122  double w = M_PI / (double)(mFFTSize - 1);
123  for (UInt32 i = 0; i < mFFTSize; ++i)
124  {
125  mWindow[i] = sin(w * (double)i);
126  }
127 }
128 
129 void CASpectralProcessor::Process(UInt32 inNumFrames, AudioBufferList* inInput, AudioBufferList* outOutput)
130 {
131  // copy from buffer list to input buffer
132  CopyInput(inNumFrames, inInput);
133 
134  // if enough input to process, then process.
135  while (mInputSize >= mFFTSize)
136  {
137  CopyInputToFFT(); // copy from input buffer to fft buffer
138  DoWindowing();
139  DoFwdFFT();
140  ProcessSpectrum(mFFTSize, mSpectralBufferList());
141  DoInvFFT();
142  DoWindowing();
143  OverlapAddOutput();
144  }
145 
146  // copy from output buffer to buffer list
147  CopyOutput(inNumFrames, outOutput);
148 }
149 
150 void CASpectralProcessor::DoWindowing()
151 {
152  Float32 *win = mWindow();
153  if (!win) return;
154  for (UInt32 i=0; i<mNumChannels; ++i) {
155  Float32 *x = mChannels[i].mFFTBuf();
156  vDSP_vmul(x, 1, win, 1, x, 1, mFFTSize);
157  }
158  //printf("DoWindowing %g %g\n", mChannels[0].mFFTBuf()[0], mChannels[0].mFFTBuf()[200]);
159 }
160 
161 
162 
163 void CASpectralProcessor::CopyInput(UInt32 inNumFrames, AudioBufferList* inInput)
164 {
165  UInt32 numBytes = inNumFrames * sizeof(Float32);
166  UInt32 firstPart = mIOBufSize - mInputPos;
167 
168 
169  if (firstPart < inNumFrames) {
170  UInt32 firstPartBytes = firstPart * sizeof(Float32);
171  UInt32 secondPartBytes = numBytes - firstPartBytes;
172  for (UInt32 i=0; i<mNumChannels; ++i) {
173  memcpy(mChannels[i].mInputBuf + mInputPos, inInput->mBuffers[i].mData, firstPartBytes);
174  memcpy(mChannels[i].mInputBuf, (UInt8*)inInput->mBuffers[i].mData + firstPartBytes, secondPartBytes);
175  }
176  } else {
177  UInt32 numBytes = inNumFrames * sizeof(Float32);
178  for (UInt32 i=0; i<mNumChannels; ++i) {
179  memcpy(mChannels[i].mInputBuf + mInputPos, inInput->mBuffers[i].mData, numBytes);
180  }
181  }
182  //printf("CopyInput %g %g\n", mChannels[0].mInputBuf[mInputPos], mChannels[0].mInputBuf[(mInputPos + 200) & mIOMask]);
183  //printf("CopyInput mInputPos %u mIOBufSize %u\n", (unsigned)mInputPos, (unsigned)mIOBufSize);
184  mInputSize += inNumFrames;
185  mInputPos = (mInputPos + inNumFrames) & mIOMask;
186 }
187 
188 void CASpectralProcessor::CopyOutput(UInt32 inNumFrames, AudioBufferList* outOutput)
189 {
190  //printf("->CopyOutput %g %g\n", mChannels[0].mOutputBuf[mOutputPos], mChannels[0].mOutputBuf[(mOutputPos + 200) & mIOMask]);
191  //printf("CopyOutput mOutputPos %u\n", (unsigned)mOutputPos);
192  UInt32 numBytes = inNumFrames * sizeof(Float32);
193  UInt32 firstPart = mIOBufSize - mOutputPos;
194  if (firstPart < inNumFrames) {
195  UInt32 firstPartBytes = firstPart * sizeof(Float32);
196  UInt32 secondPartBytes = numBytes - firstPartBytes;
197  for (UInt32 i=0; i<mNumChannels; ++i) {
198  memcpy(outOutput->mBuffers[i].mData, mChannels[i].mOutputBuf + mOutputPos, firstPartBytes);
199  memcpy((UInt8*)outOutput->mBuffers[i].mData + firstPartBytes, mChannels[i].mOutputBuf, secondPartBytes);
200  memset(mChannels[i].mOutputBuf + mOutputPos, 0, firstPartBytes);
201  memset(mChannels[i].mOutputBuf, 0, secondPartBytes);
202  }
203  } else {
204  for (UInt32 i=0; i<mNumChannels; ++i) {
205  memcpy(outOutput->mBuffers[i].mData, mChannels[i].mOutputBuf + mOutputPos, numBytes);
206  memset(mChannels[i].mOutputBuf + mOutputPos, 0, numBytes);
207  }
208  }
209  //printf("<-CopyOutput %g %g\n", ((Float32*)outOutput->mBuffers[0].mData)[0], ((Float32*)outOutput->mBuffers[0].mData)[200]);
210  mOutputPos = (mOutputPos + inNumFrames) & mIOMask;
211 }
212 
213 void CASpectralProcessor::PrintSpectralBufferList()
214 {
215  UInt32 half = mFFTSize >> 1;
216  for (UInt32 i=0; i<mNumChannels; ++i) {
217  DSPSplitComplex &freqData = mSpectralBufferList->mDSPSplitComplex[i];
218 
219  for (UInt32 j=0; j<half; j++){
220  printf(" bin[%d]: %lf + %lfi\n", (int) j, freqData.realp[j], freqData.imagp[j]);
221  }
222  }
223 }
224 
225 
226 void CASpectralProcessor::CopyInputToFFT()
227 {
228  //printf("CopyInputToFFT mInFFTPos %u\n", (unsigned)mInFFTPos);
229  UInt32 firstPart = mIOBufSize - mInFFTPos;
230  UInt32 firstPartBytes = firstPart * sizeof(Float32);
231  if (firstPartBytes < mFFTByteSize) {
232  UInt32 secondPartBytes = mFFTByteSize - firstPartBytes;
233  for (UInt32 i=0; i<mNumChannels; ++i) {
234  memcpy(mChannels[i].mFFTBuf(), mChannels[i].mInputBuf() + mInFFTPos, firstPartBytes);
235  memcpy((UInt8*)mChannels[i].mFFTBuf() + firstPartBytes, mChannels[i].mInputBuf(), secondPartBytes);
236  }
237  } else {
238  for (UInt32 i=0; i<mNumChannels; ++i) {
239  memcpy(mChannels[i].mFFTBuf(), mChannels[i].mInputBuf() + mInFFTPos, mFFTByteSize);
240  }
241  }
242  mInputSize -= mHopSize;
243  mInFFTPos = (mInFFTPos + mHopSize) & mIOMask;
244  //printf("CopyInputToFFT %g %g\n", mChannels[0].mFFTBuf()[0], mChannels[0].mFFTBuf()[200]);
245 }
246 
247 void CASpectralProcessor::OverlapAddOutput()
248 {
249  //printf("OverlapAddOutput mOutFFTPos %u\n", (unsigned)mOutFFTPos);
250  UInt32 firstPart = mIOBufSize - mOutFFTPos;
251  if (firstPart < mFFTSize) {
252  UInt32 secondPart = mFFTSize - firstPart;
253  for (UInt32 i=0; i<mNumChannels; ++i) {
254  float* out1 = mChannels[i].mOutputBuf() + mOutFFTPos;
255  vDSP_vadd(out1, 1, mChannels[i].mFFTBuf(), 1, out1, 1, firstPart);
256  float* out2 = mChannels[i].mOutputBuf();
257  vDSP_vadd(out2, 1, mChannels[i].mFFTBuf() + firstPart, 1, out2, 1, secondPart);
258  }
259  } else {
260  for (UInt32 i=0; i<mNumChannels; ++i) {
261  float* out1 = mChannels[i].mOutputBuf() + mOutFFTPos;
262  vDSP_vadd(out1, 1, mChannels[i].mFFTBuf(), 1, out1, 1, mFFTSize);
263  }
264  }
265  //printf("OverlapAddOutput %g %g\n", mChannels[0].mOutputBuf[mOutFFTPos], mChannels[0].mOutputBuf[(mOutFFTPos + 200) & mIOMask]);
266  mOutFFTPos = (mOutFFTPos + mHopSize) & mIOMask;
267 }
268 
269 
270 void CASpectralProcessor::DoFwdFFT()
271 {
272  //printf("->DoFwdFFT %g %g\n", mChannels[0].mFFTBuf()[0], mChannels[0].mFFTBuf()[200]);
273  UInt32 half = mFFTSize >> 1;
274  for (UInt32 i=0; i<mNumChannels; ++i)
275  {
276  vDSP_ctoz((DSPComplex*)mChannels[i].mFFTBuf(), 2, &mSpectralBufferList->mDSPSplitComplex[i], 1, half);
277  vDSP_fft_zrip(mFFTSetup, &mSpectralBufferList->mDSPSplitComplex[i], 1, mLog2FFTSize, FFT_FORWARD);
278  }
279  //printf("<-DoFwdFFT %g %g\n", direction, mChannels[0].mFFTBuf()[0], mChannels[0].mFFTBuf()[200]);
280 }
281 
282 void CASpectralProcessor::DoInvFFT()
283 {
284  //printf("->DoInvFFT %g %g\n", mChannels[0].mFFTBuf()[0], mChannels[0].mFFTBuf()[200]);
285  UInt32 half = mFFTSize >> 1;
286  for (UInt32 i=0; i<mNumChannels; ++i)
287  {
288  vDSP_fft_zrip(mFFTSetup, &mSpectralBufferList->mDSPSplitComplex[i], 1, mLog2FFTSize, FFT_INVERSE);
289  vDSP_ztoc(&mSpectralBufferList->mDSPSplitComplex[i], 1, (DSPComplex*)mChannels[i].mFFTBuf(), 2, half);
290  float scale = 0.5 / mFFTSize;
291  vDSP_vsmul(mChannels[i].mFFTBuf(), 1, &scale, mChannels[i].mFFTBuf(), 1, mFFTSize );
292  }
293  //printf("<-DoInvFFT %g %g\n", direction, mChannels[0].mFFTBuf()[0], mChannels[0].mFFTBuf()[200]);
294 }
295 
296 void CASpectralProcessor::SetSpectralFunction(SpectralFunction inFunction, void* inUserData)
297 {
298  mSpectralFunction = inFunction;
299  mUserData = inUserData;
300 }
301 
302 void CASpectralProcessor::ProcessSpectrum(UInt32 inFFTSize, SpectralBufferList* inSpectra)
303 {
304  if (mSpectralFunction)
305  (mSpectralFunction)(inSpectra, mUserData);
306 }
307 
308 #pragma mark ___Utility___
309 
310 void CASpectralProcessor::GetMagnitude(AudioBufferList* list, Float32* min, Float32* max)
311 {
312  UInt32 half = mFFTSize >> 1;
313  for (UInt32 i=0; i<mNumChannels; ++i) {
314  DSPSplitComplex &freqData = mSpectralBufferList->mDSPSplitComplex[i];
315 
316  Float32* b = (Float32*) list->mBuffers[i].mData;
317 
318  vDSP_zvabs(&freqData,1,b,1,half);
319 
320  vDSP_maxmgv(b, 1, &max[i], half);
321  vDSP_minmgv(b, 1, &min[i], half);
322 
323  }
324 }
325 
326 
327 void CASpectralProcessor::GetFrequencies(Float32* freqs, Float32 sampleRate)
328 {
329  UInt32 half = mFFTSize >> 1;
330 
331  for (UInt32 i=0; i< half; i++){
332  freqs[i] = ((Float32)(i))*sampleRate/((Float32)mFFTSize);
333  }
334 }
335 
336 
337 bool CASpectralProcessor::ProcessForwards(UInt32 inNumFrames, AudioBufferList* inInput)
338 {
339  // copy from buffer list to input buffer
340  CopyInput(inNumFrames, inInput);
341 
342  bool processed = false;
343  // if enough input to process, then process.
344  while (mInputSize >= mFFTSize)
345  {
346  CopyInputToFFT(); // copy from input buffer to fft buffer
347  DoWindowing();
348  DoFwdFFT();
349  ProcessSpectrum(mFFTSize, mSpectralBufferList()); // here you would copy the fft results out to a buffer indicated in mUserData, say for sonogram drawing
350  processed = true;
351  }
352 
353  return processed;
354 }
355 
356 bool CASpectralProcessor::ProcessBackwards(UInt32 inNumFrames, AudioBufferList* outOutput)
357 {
358 
359  ProcessSpectrum(mFFTSize, mSpectralBufferList());
360  DoInvFFT();
361  DoWindowing();
362  OverlapAddOutput();
363 
364  // copy from output buffer to buffer list
365  CopyOutput(inNumFrames, outOutput);
366 
367  return true;
368 }
369 
370