Retired Document
Important: This sample code may not represent best practices for current development. The project may use deprecated symbols and illustrate technologies and techniques that are no longer recommended.
Test vBigDSP.c
/* |
Disclaimer: IMPORTANT: This Apple software is supplied to you by Apple Computer, Inc. |
("Apple") in consideration of your agreement to the following terms, and your |
use, installation, modification or redistribution of this Apple software |
constitutes acceptance of these terms. If you do not agree with these terms, |
please do not use, install, modify or redistribute this Apple software. |
In consideration of your agreement to abide by the following terms, and subject |
to these terms, Apple grants you a personal, non-exclusive license, under AppleÕs |
copyrights in this original Apple software (the "Apple Software"), to use, |
reproduce, modify and redistribute the Apple Software, with or without |
modifications, in source and/or binary forms; provided that if you redistribute |
the Apple Software in its entirety and without modifications, you must retain |
this notice and the following text and disclaimers in all such redistributions of |
the Apple Software. Neither the name, trademarks, service marks or logos of |
Apple Computer, Inc. may be used to endorse or promote products derived from the |
Apple Software without specific prior written permission from Apple. Except as |
expressly stated in this notice, no other rights or licenses, express or implied, |
are granted by Apple herein, including but not limited to any patent rights that |
may be infringed by your derivative works or by other works in which the Apple |
Software may be incorporated. |
The Apple Software is provided by Apple on an "AS IS" basis. APPLE MAKES NO |
WARRANTIES, EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION THE IMPLIED |
WARRANTIES OF NON-INFRINGEMENT, MERCHANTABILITY AND FITNESS FOR A PARTICULAR |
PURPOSE, REGARDING THE APPLE SOFTWARE OR ITS USE AND OPERATION ALONE OR IN |
COMBINATION WITH YOUR PRODUCTS. |
IN NO EVENT SHALL APPLE BE LIABLE FOR ANY SPECIAL, INDIRECT, INCIDENTAL OR |
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE |
GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) |
ARISING IN ANY WAY OUT OF THE USE, REPRODUCTION, MODIFICATION AND/OR DISTRIBUTION |
OF THE APPLE SOFTWARE, HOWEVER CAUSED AND WHETHER UNDER THEORY OF CONTRACT, TORT |
(INCLUDING NEGLIGENCE), STRICT LIABILITY OR OTHERWISE, EVEN IF APPLE HAS BEEN |
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
*/ |
#include <MacMemory.h> |
#include <math.h> |
#include <stdio.h> |
#include "vBigDSP.h" |
#ifndef PI |
#define PI (3.1415926535897932384626433832795) |
#endif |
#pragma mark - |
static double |
RMSE_real( |
float * x, |
float * y, |
int n |
) |
/* Calculates the root-mean-squared error between two real signals |
* x[n] and y[n]. |
*/ |
{ |
double t, err=0.0; |
int i; |
for (i=0; i<n; i++) |
{ |
t = x[i] - y[i]; |
err += (t * t); |
} |
return(sqrt(err/n)); |
} |
static double |
RMSE_complex( |
float * x, |
float * y, |
int n |
) |
/* Calculates the root-mean-squared error between two complex signals |
* x[n] and y[n]. |
*/ |
{ |
double t, err=0.0; |
int i; |
for (i=0; i<n; i++) |
{ |
t = x[2*i] - y[2*i]; |
err += t * t; |
t = x[2*i+1] - y[2*i+1]; |
err += t * t; |
} |
return(sqrt(err/n)); |
} |
#pragma mark - |
static void |
ConvolveRealLiteral( |
float * x, |
float * y, |
float * z, |
int n |
) |
/* Performs a literal cyclic convolution on x and y, placing the |
* result in z. This is an O(n^2) algorithm useful |
* for accuracy testing. |
*/ |
{ |
int s, i, q; |
for (s = 0; s < n; s++) |
{ |
z[s] = 0; |
for (q = 0; q < n; q++) |
{ |
i = (s-q)%n; |
if(i<0) i+= n; |
z[s] += y[i] * x[q]; |
} |
} |
} |
static void ConvolveComplexLiteral( |
float * x, |
float * y, |
float * z, |
int n |
) |
/* Performs a literal cyclic convolution on x and y, placing the |
* result in z. This is an O(n^2) algorithm useful |
* for accuracy testing. |
*/ |
{ |
int s, i, q; |
for (s = 0; s < n; s++) |
{ |
z[2*s] = z[2*s+1] = 0; |
for (q = 0; q < n; q++) |
{ |
i = (s-q)%n; |
if(i<0) i+= n; |
z[2*s] += y[2*i] * x[2*q] - y[2*i+1] * x[2*q+1]; |
z[2*s+1] += y[2*i+1] * x[2*q] + y[2*i] * x[2*q+1]; |
} |
} |
} |
static void |
Convolve2DRealLiteral( |
float * x, |
float * y, |
float * z, |
int w, |
int h |
) |
/* Literal cyclic convolution in two dimensions, for testing purposes. The |
* convolution of x[] and y[] is output in z[]. |
*/ |
{ |
int k, j, q, p, kk, jj; |
for (k=0; k<h; k++) |
{ |
for (j=0; j<w; j++) |
{ |
z[j + k*w] = 0; |
for (q = 0; q < h; q++) |
{ |
kk = (k-q<0)?k-q+h:k-q; |
for (p = 0; p < w; p++) |
{ |
jj = (j-p<0)?j-p+w:j-p; |
z[j + k*w] += x[p + q*w] * y[jj + kk*w]; |
} |
} |
} |
} |
} |
#pragma mark - |
#define COMPLEX_TEST_MIN_POWER 0 |
#define COMPLEX_TEST_MAX_POWER 18 |
static void TestComplexFFT() |
{ |
float *real_data; |
float *vector_data; |
int i, j; |
long currentLength; |
OSErr result; |
double rmse_vector; |
for (i = COMPLEX_TEST_MIN_POWER; i<=COMPLEX_TEST_MAX_POWER; i++) { |
currentLength = 1 << i; |
real_data = (float*)NewPtr(2*currentLength*sizeof(float)); |
vector_data = (float*)NewPtr(2*currentLength*sizeof(float)); |
if (real_data == nil || vector_data == nil) { |
printf("\nerror allocating space for complex test"); |
} else { |
///////////////////////////////////////////////// |
// initialize signal |
///////////////////////////////////////////////// |
for (j=0; j<currentLength; j++) { |
real_data[2*j] = vector_data[2*j] = |
cos(PI*j*j/currentLength); |
real_data[2*j+1] = vector_data[2*j+1] = |
sin(PI*j*j/currentLength); |
} |
///////////////////////////////////////////////// |
// perform forward and inverse FFT and compare |
// result with original signal |
///////////////////////////////////////////////// |
result = FFTComplex(vector_data, currentLength, -1); |
if (result != noErr) { |
printf("\nerror in forward complex FFT"); |
} |
result = FFTComplex(vector_data, currentLength, 1); |
if (result != noErr) { |
printf("\nerror in inverse complex FFT"); |
} |
rmse_vector = RMSE_complex(real_data, vector_data, currentLength); |
printf("\n complex length %d FFT rmse: %g", currentLength, rmse_vector); |
} |
if (real_data) DisposePtr((Ptr)real_data); |
if (vector_data) DisposePtr((Ptr)vector_data); |
} |
} |
#define REAL_TEST_MIN_POWER 0 |
#define REAL_TEST_MAX_POWER 18 |
static void TestRealFFT() |
{ |
float *real_data; |
float *vector_data; |
int i, j; |
double negativeOneToI = -1; |
double rmse_vector; |
long currentLength; |
OSErr result; |
for (i=REAL_TEST_MIN_POWER; i<= REAL_TEST_MAX_POWER; i++) { |
currentLength = 1 << i; |
real_data = (float*)NewPtr(currentLength*sizeof(float)); |
vector_data = (float*)NewPtr(currentLength*sizeof(float)); |
if (real_data == nil || vector_data == nil) { |
printf("\nerror allocating space for complex test"); |
} else { |
///////////////////////////////////////////////// |
// initialize signal |
///////////////////////////////////////////////// |
for (j=0; j<currentLength; j++) { |
negativeOneToI = -negativeOneToI; |
real_data[j] = vector_data[j] = |
negativeOneToI + sin(3*PI*j/currentLength); |
} |
///////////////////////////////////////////////// |
// perform forward and inverse FFT and compare |
// result with original signal |
///////////////////////////////////////////////// |
result = FFTRealForward(vector_data, currentLength); |
if (result != noErr) { |
printf("\nerror in forward real FFT"); |
} |
result = FFTRealInverse(vector_data, currentLength); |
if (result != noErr) { |
printf("\nerror in forward real FFT"); |
} |
rmse_vector = RMSE_real(real_data, vector_data, currentLength); |
printf("\n real length %d FFT rmse: %g", currentLength, rmse_vector); |
} |
if (real_data) DisposePtr((Ptr)real_data); |
if (vector_data) DisposePtr((Ptr)vector_data); |
} |
} |
#define COMPLEX_CONVOLVE_MIN_POWER 2 |
#define COMPLEX_CONVOLVE_MAX_POWER 13 |
static void TestComplexConvolve() |
{ |
float *data1; |
float *data2; |
float *literalConvolveData; |
int i, j; |
double negativeOneToI = -1; |
double rmse_vector; |
long currentLength; |
OSErr result; |
for (i=COMPLEX_CONVOLVE_MIN_POWER; i <= COMPLEX_CONVOLVE_MAX_POWER; i++) { |
currentLength = 1 << i; |
data1 = (float*)NewPtr(2*currentLength*sizeof(float)); |
data2 = (float*)NewPtr(2*currentLength*sizeof(float)); |
literalConvolveData = (float*)NewPtr(2*currentLength*sizeof(float)); |
if (data1 == nil || data2 == nil || literalConvolveData == nil) { |
printf("\nerror allocating space for complex convolve test"); |
} else { |
///////////////////////////////////////////////// |
// initialize signals |
///////////////////////////////////////////////// |
for (j=0; j<currentLength; j++) { |
data1[2*j] = cos(PI*j*j/currentLength); |
data1[2*j+1] = sin(PI*j*j/currentLength); |
data2[2*j] = 2*sin(3.3*PI*j*j/currentLength); |
data2[2*j+1] = -cos(1.4*PI*j*j/currentLength); |
} |
///////////////////////////////////////////////// |
// perform literal convolution, and compare |
// with convolution performed through FFT |
///////////////////////////////////////////////// |
ConvolveComplexLiteral(data1, data2, literalConvolveData, currentLength); |
result = ConvolveComplexAltivec(data1, data2, currentLength); |
if (result != noErr) { |
printf("\nerror in convolve complex"); |
} |
rmse_vector = RMSE_complex(data2, literalConvolveData, currentLength); |
printf("\n complex length %d convolve rmse: %g", currentLength, rmse_vector); |
} |
if (data1) DisposePtr((Ptr)data1); |
if (data2) DisposePtr((Ptr)data2); |
if (literalConvolveData) DisposePtr((Ptr)literalConvolveData); |
} |
} |
#define REAL_CONVOLVE_MIN_POWER 0 |
#define REAL_CONVOLVE_MAX_POWER 14 |
static void TestRealConvolve() |
{ |
float *data1; |
float *data2; |
float *literalConvolveData; |
int i, j; |
double negativeOneToI = -1; |
double rmse_vector; |
long currentLength; |
OSErr result; |
for (i=REAL_CONVOLVE_MIN_POWER; i <= REAL_CONVOLVE_MAX_POWER; i++) { |
currentLength = 1 << i; |
data1 = (float*)NewPtr(currentLength*sizeof(float)); |
data2 = (float*)NewPtr(currentLength*sizeof(float)); |
literalConvolveData = (float*)NewPtr(currentLength*sizeof(float)); |
if (data1 == nil || data2 == nil || literalConvolveData == nil) { |
printf("\nerror allocating space for real convolve test"); |
} else { |
///////////////////////////////////////////////// |
// initialize signals |
///////////////////////////////////////////////// |
for (j=0; j<currentLength; j++) { |
negativeOneToI = -negativeOneToI; |
data1[j] = |
negativeOneToI + sin(3*PI*j/currentLength); |
data2[j] = |
3 * negativeOneToI + sin(4.4*PI*j/currentLength); |
} |
///////////////////////////////////////////////// |
// perform literal convolution, and compare |
// with convolution performed through FFT |
///////////////////////////////////////////////// |
ConvolveRealLiteral(data1, data2, literalConvolveData, currentLength); |
result = ConvolveRealAltivec(data1, data2, currentLength); |
if (result != noErr) { |
printf("\nerror in convolve real"); |
} |
rmse_vector = RMSE_real(data2, literalConvolveData, currentLength); |
printf("\n real length %d convolve rmse: %g", currentLength, rmse_vector); |
} |
if (data1) DisposePtr((Ptr)data1); |
if (data2) DisposePtr((Ptr)data2); |
if (literalConvolveData) DisposePtr((Ptr)literalConvolveData); |
} |
} |
#define COMPLEX_FFT2D_MIN_POWER_X 1 |
#define COMPLEX_FFT2D_MAX_POWER_X 8 |
#define COMPLEX_FFT2D_MIN_POWER_Y 1 |
#define COMPLEX_FFT2D_MAX_POWER_Y 8 |
static void TestFFT2DComplex() |
{ |
float *data1; |
float *data2; |
int xsize, ysize, j; |
double negativeOneToI = -1; |
double rmse_vector; |
long currentLength; |
OSErr result; |
for (xsize=COMPLEX_FFT2D_MIN_POWER_X; xsize <= COMPLEX_FFT2D_MAX_POWER_X; xsize++) { |
for (ysize = COMPLEX_FFT2D_MIN_POWER_Y; ysize <= COMPLEX_FFT2D_MAX_POWER_Y; ysize++) { |
currentLength = (1 << xsize) * (1 << ysize); |
data1 = (float*)NewPtr(2*currentLength*sizeof(float)); |
data2 = (float*)NewPtr(2*currentLength*sizeof(float)); |
if (data1 == nil || data2 == nil ) { |
printf("\nerror allocating space for FFT 2D test"); |
} else { |
///////////////////////////////////////////////// |
// initialize signal |
///////////////////////////////////////////////// |
for (j=0; j<currentLength; j++) { |
data1[2*j] = data2[2*j] = |
cos(PI*j*j/currentLength); |
data1[2*j+1] = data2[2*j+1] = |
sin(PI*j*j/currentLength); |
} |
///////////////////////////////////////////////// |
// perform forward and inverse FFT and compare |
// result with original signal |
///////////////////////////////////////////////// |
result = FFT2DComplex(data1, 1 << xsize, 1 << ysize, -1); |
if (result != noErr) { |
printf("\nerror in fft 2d complex"); |
} |
result = FFT2DComplex(data1, 1 << xsize, 1 << ysize, 1); |
if (result != noErr) { |
printf("\nerror in fft 2d complex"); |
} |
rmse_vector = RMSE_complex(data1, data2, currentLength); |
printf("\n complex 2D %d X %d FFT rmse: %g", 1 << xsize, 1 << ysize, rmse_vector); |
} |
if (data1) DisposePtr((Ptr)data1); |
if (data2) DisposePtr((Ptr)data2); |
} |
} |
} |
#define REAL_FFT2D_MIN_POWER_X 3 |
#define REAL_FFT2D_MAX_POWER_X 9 |
#define REAL_FFT2D_MIN_POWER_Y 2 |
#define REAL_FFT2D_MAX_POWER_Y 9 |
static void TestFFT2DReal() |
{ |
float *data1; |
float *data2; |
int xsize, ysize, j; |
double negativeOneToI = -1; |
double rmse_vector; |
long currentLength; |
OSErr result; |
for (xsize=REAL_FFT2D_MIN_POWER_X; xsize <= REAL_FFT2D_MAX_POWER_X; xsize++) { |
for (ysize = REAL_FFT2D_MIN_POWER_Y; ysize <= REAL_FFT2D_MAX_POWER_Y; ysize++) { |
currentLength = (1 << xsize) * (1 << ysize); |
data1 = (float*)NewPtr(2*currentLength*sizeof(float)); |
data2 = (float*)NewPtr(2*currentLength*sizeof(float)); |
if (data1 == nil || data2 == nil ) { |
printf("\nerror allocating space for FFT 2D test"); |
} else { |
///////////////////////////////////////////////// |
// initialize signal |
///////////////////////////////////////////////// |
for (j=0; j<currentLength; j++) { |
negativeOneToI = -negativeOneToI; |
data1[j] = data2[j] = |
negativeOneToI + sin(3*PI*j/currentLength); |
} |
///////////////////////////////////////////////// |
// perform forward and inverse FFT and compare |
// result with original signal |
///////////////////////////////////////////////// |
result = FFT2DRealForward(data1, 1 << xsize, 1 << ysize); |
if (result != noErr) { |
printf("\nerror in fft 2d real forward"); |
} |
result = FFT2DRealInverse(data1, 1 << xsize, 1 << ysize); |
if (result != noErr) { |
printf("\nerror in fft 2d real inverse"); |
} |
rmse_vector = RMSE_real(data1, data2, currentLength); |
printf("\n real 2D %d X %d FFT rmse: %g", 1 << xsize, 1 << ysize, rmse_vector); |
} |
if (data1) DisposePtr((Ptr)data1); |
if (data2) DisposePtr((Ptr)data2); |
} |
} |
} |
#define COMPLEX_CONVOLVE2D_MIN_POWER_X 2 |
#define COMPLEX_CONVOLVE2D_MAX_POWER_X 5 |
#define COMPLEX_CONVOLVE2D_MIN_POWER_Y 2 |
#define COMPLEX_CONVOLVE2D_MAX_POWER_Y 5 |
#define REAL_CONVOLVE2D_MIN_POWER_X 3 |
#define REAL_CONVOLVE2D_MAX_POWER_X 8 |
#define REAL_CONVOLVE2D_MIN_POWER_Y 2 |
#define REAL_CONVOLVE2D_MAX_POWER_Y 8 |
static void TestConvolve2DReal() |
{ |
float *data1; |
float *data2; |
float *literalConvolve; |
int xsize, ysize, j; |
double negativeOneToI = -1; |
double rmse_vector; |
long currentLength; |
OSErr result; |
for (xsize=REAL_CONVOLVE2D_MIN_POWER_X; xsize <= REAL_CONVOLVE2D_MAX_POWER_X; xsize++) { |
for (ysize = REAL_CONVOLVE2D_MIN_POWER_Y; ysize <= REAL_CONVOLVE2D_MAX_POWER_Y; ysize++) { |
currentLength = (1 << xsize) * (1 << ysize); |
data1 = (float*)NewPtr(currentLength*sizeof(float)); |
data2 = (float*)NewPtr(currentLength*sizeof(float)); |
literalConvolve = (float*)NewPtr(currentLength*sizeof(float)); |
if (data1 == nil || data2 == nil || literalConvolve == nil) { |
printf("\nerror allocating space for convolve 2D real test"); |
} else { |
///////////////////////////////////////////////// |
// initialize signals |
///////////////////////////////////////////////// |
for (j=0; j<currentLength; j++) { |
negativeOneToI = -negativeOneToI; |
data1[j] = negativeOneToI + sin(0.05*PI*j/currentLength); |
data2[j] = -1.3 * cos(1.3*PI*j/currentLength); |
} |
///////////////////////////////////////////////// |
// perform literal convolution, and compare |
// with convolution performed through FFT |
///////////////////////////////////////////////// |
Convolve2DRealLiteral(data1, data2, literalConvolve, 1 << xsize, 1 << ysize); |
result = ConvolveRealAltivec2D(data1, data2, 1 << xsize, 1 << ysize); |
if (result != noErr) { |
printf("\nerror in fft 2d real"); |
} |
rmse_vector = RMSE_real(data2, literalConvolve, currentLength); |
printf("\n real 2D %d X %d convolve rmse: %g", 1 << xsize, 1 << ysize, rmse_vector); |
} |
if (data1) DisposePtr((Ptr)data1); |
if (data2) DisposePtr((Ptr)data2); |
if (literalConvolve) DisposePtr((Ptr)literalConvolve); |
} |
} |
} |
#pragma mark - |
int main() { |
//////////////////////////////////////////////////////////////////////////////// |
printf("\n\nTesting complex FFT:\n"); |
TestComplexFFT(); |
//////////////////////////////////////////////////////////////////////////////// |
printf("\n\nTesting real FFT:\n"); |
TestRealFFT(); |
//////////////////////////////////////////////////////////////////////////////// |
printf("\n\nTesting 2D complex FFT:\n"); |
TestFFT2DComplex(); |
//////////////////////////////////////////////////////////////////////////////// |
printf("\n\nTesting 2D real FFT:\n"); |
TestFFT2DReal(); |
//////////////////////////////////////////////////////////////////////////////// |
printf("\n\nTesting complex convolution:\n"); |
TestComplexConvolve(); |
//////////////////////////////////////////////////////////////////////////////// |
printf("\n\nTesting real convolution:\n"); |
TestRealConvolve(); |
//////////////////////////////////////////////////////////////////////////////// |
printf("\n\nTesting 2D real convolution:\n"); |
TestConvolve2DReal(); |
} |
Copyright © 2003 Apple Computer, Inc. All Rights Reserved. Terms of Use | Privacy Policy | Updated: 2003-01-14