Update of /cvsroot/audacity/lib-src/libnyquist/nyquist/ffts/Matlab-testing
In directory 23jxhf1.ch3.sourceforge.com:/tmp/cvs-serv8506/ffts/Matlab-testing
Added Files:
conv2dTest.c conv2dtest.m convTest.c convtest.m rfft2dTestML.c
rfft2dTestML.m
Log Message:
Updating to Nyquist v3.03.
--- NEW FILE: conv2dtest.m ---
% program to test 2d real fast conv
% let user select file then open it
[fname, pname] = uigetfile('*.c2d', 'select conv file');
cd(pname);
fidout=fopen(fname,'r');
% read header info
aN=fread(fidout,1,'long');
aM=fread(fidout,1,'long');
bN=fread(fidout,1,'long');
bM=fread(fidout,1,'long');
% read in data
%status=fseek(fidout,Nheader,'bof');
a=fread(fidout,aN*aM,'float');
a=reshape(a,aN,aM);
b=fread(fidout,bN*bM,'float');
b=reshape(b,bN,bM);
c=fread(fidout,(aN+bN-1)*(aM+bM-1),'float');
c=reshape(c,(aN+bN-1),(aM+bM-1));
fclose(fidout);
c2=conv2(a,b);
max(max(abs(c2-c)))
--- NEW FILE: rfft2dTestML.m ---
% program to test 2d real fft
% let user select file then open it
[fname, pname] = uigetfile('*.dr2', 'select conv file');
cd(pname);
fidout=fopen(fname,'r');
% read header info
N=fread(fidout,1,'long');
M=fread(fidout,1,'long');
% read in data
%status=fseek(fidout,Nheader,'bof');
a=fread(fidout,N*M,'float');
a=reshape(a,N,M);
c=fread(fidout,N*M,'float');
c=reshape(c,N,M);
c=c(1:2:N,:)+j*c(2:2:N,:);
fclose(fidout);
c2=fft2(a);
% Remember Matlab is column major order, C is row major order
% Matlab starts with index of 1, f and k start at 0
%%%% check the real transforms of DC and Nyquest frequencies
% check the four real values
maxerr=abs(real(c(1,1))-c2(1,1)); % 0 f,
0 k
maxerr=max(maxerr,abs(imag(c(1,1))-c2(1,M/2+1))); % 0 f, M/2 k
maxerr=max(maxerr,abs(real(c(1,M/2+1))-c2(N/2+1,1)));% N/2 f, 0 k
maxerr=max(maxerr,abs(imag(c(1,M/2+1))-c2(N/2+1,M/2+1)));% N/2 f, M/2 k
%check the rest of the pos wavenumbers of DC and Nyquest frequencies
maxerr=max(maxerr,max(abs(c(1,2:M/2)-c2(1,2:M/2)))); %f = 0, k=1 to M/2-1
maxerr=max(maxerr,max(abs(c(1,M/2+2:M)-c2(N/2+1,2:M/2))));%f = N/2, k=1 to M/2-1
%%%%
% check all the other positive frequencies at all wavenumbers
% f from 1 to N/2-1, k from 0 to M-1 (wraps around through negative k)
maxerr=max(maxerr,max(max(abs(c(2:N/2,:)-c2(2:N/2,:)))))
--- NEW FILE: conv2dTest.c ---
/* A program to test fast 2d real convolution */
#include <stdio.h>
#include <stdlib.h>
#include <fp.h>
#include <math.h>
#include "fftlib.h"
#include "fftext.h"
#include "fft2d.h"
#if macintosh
#include <timer.h>
#endif
#define BIPRAND(a) (2.0/(RAND_MAX+1.0)*a-1.0)
//#define BIPRAND(a) round(100*(2.0/(RAND_MAX+1.0)*a-1.0))
void main(){
long N = 256; /* the number of cols in 2d ffts, must be power of 2 */
long N2 = 64; /* the number of rows in 2d ffts, must be power of 2 */
long kernSize = 53; /* kernal cols must be less than N */
long kernSize2 = 29; /* kernal rows must be less than N2*/
long dataSize = N-kernSize+1; /* data cols */
long dataSize2 = N2-kernSize2+1; /* data rows */
float *a;
float *b;
long i1;
long i2;
long TheErr;
long M;
long M2;
FILE *fdataout; /* output file */
unsigned int randseed = 777;
int rannum;
#if macintosh
UnsignedWide TheTime1;
Microseconds(&TheTime1);
randseed = TheTime1.lo;
#endif
printf(" %6d Byte Floats \n", sizeof(a[0]));
printf(" randseed = %10u\n", randseed);
srand(randseed);
M = roundtol(LOG2(N));
N = POW2(M);
M2 = roundtol(LOG2(N2));
N2 = POW2(M2);
printf("fft size = %6d X%6d, ", N2, N);
if ((dataSize <= 0)||(dataSize2 <= 0)) TheErr = 22;
else TheErr = 0;
if(!TheErr){
TheErr = fft2dInit(M2, M);
}
a = (float *) calloc(N2*N,sizeof(float) ); // calloc to zero pad data to
fill N to N2
if (a == 0) TheErr = 2;
if(!TheErr){
b = (float *) calloc(N2*N,sizeof(float) ); // calloc to zero pad
data to fill N to N2
if (b == 0) TheErr = 2;
}
if(!TheErr){
fdataout = fopen("conv2ddat.c2d", "wb");
if (fdataout == NULL) TheErr = -50;
}
if(!TheErr){
/* write sizes to fdataout */
fwrite(&dataSize, sizeof(dataSize), 1, fdataout);
fwrite(&dataSize2, sizeof(dataSize2), 1, fdataout);
fwrite(&kernSize, sizeof(kernSize), 1, fdataout);
fwrite(&kernSize2, sizeof(kernSize2), 1, fdataout);
/* set up a simple test case and write to fdataout */
for (i2=0; i2<dataSize2; i2++){
for (i1=0; i1<dataSize; i1++){
rannum = rand();
a[i2*N+i1] = BIPRAND(rannum);
}
fwrite(&a[i2*N], dataSize*sizeof(float), 1, fdataout);
}
for (i2=0; i2<kernSize2; i2++){
for (i1=0; i1<kernSize; i1++){
rannum = rand();
b[i2*N+i1] = BIPRAND(rannum);
}
fwrite(&b[i2*N], kernSize*sizeof(float), 1, fdataout);
}
/* fast 2d convolution of zero padded sequences */
rfft2d(a, M2, M);
rfft2d(b, M2, M);
rspect2dprod(a, b, a, N2, N);
rifft2d(a, M2, M);
/* write out answer */
fwrite(a, N2*N*sizeof(float), 1, fdataout);
fclose(fdataout);
free(b);
free(a);
fft2dFree();
}
else{
if(TheErr==2) printf(" out of memory \n");
else printf(" error \n");
fft2dFree();
}
printf(" Done. \n");
return;
}
--- NEW FILE: rfft2dTestML.c ---
/* A program to test real 2d forward and inverse fast fourier transform
routines */
#include <stdio.h>
#include <stdlib.h>
#include <fp.h>
#include <math.h>
#include "fftlib.h"
#include "fftext.h"
#include "fft2d.h"
#if macintosh
#include <timer.h>
#endif
#define BIPRAND(a) (2.0/(RAND_MAX+1.0)*a-1.0)
void main(){
long N2 = 64; /* the number of rows in 2d ffts, must be power of 2 */
long N = 256; /* the number of cols in 2d ffts, must be power of 2 */
float *a;
float maxerrfft;
long i1;
long TheErr;
long M;
long M2;
FILE *fdataout; /* output file */
unsigned int randseed = 777;
int rannum;
#if macintosh
UnsignedWide TheTime1;
Microseconds(&TheTime1);
randseed = TheTime1.lo;
#endif
printf(" %6d Byte Floats \n", sizeof(a[0]));
printf(" randseed = %10u\n", randseed);
srand(randseed);
M = roundtol(LOG2(N));
N = POW2(M);
M2 = roundtol(LOG2(N2));
N2 = POW2(M2);
printf("fft size = %6d X%6d, ", N2, N);
TheErr = 0;
if(!TheErr){
TheErr = fft2dInit(M2, M);
}
a = (float *) malloc(N2*N*sizeof(float) );
if (a == 0) TheErr = 2;
if(!TheErr){
fdataout = fopen("fftdat.dr2", "wb");
if (fdataout == NULL) TheErr = -50;
}
if(!TheErr){
/* write sizes to fdataout */
fwrite(&N, sizeof(N), 1, fdataout);
fwrite(&N2, sizeof(N2), 1, fdataout);
/* set up a simple test case and write to fdataout */
for (i1=0; i1<N2*N; i1++){
rannum = rand();
a[i1] = BIPRAND(rannum);
}
fwrite(a, N2*N*sizeof(float), 1, fdataout);
/* real, 2d fast fourier transform */
rfft2d(a, M2, M);
/* write out answer */
fwrite(a, N2*N*sizeof(float), 1, fdataout);
fclose(fdataout);
/* compute and check inverse transform */
rifft2d(a, M2, M);
maxerrfft = 0;
srand(randseed);
for (i1=0; i1<N2*N; i1++){
rannum = rand();
maxerrfft = fmax(maxerrfft, fabs(BIPRAND(rannum)-a[i1]));
}
printf("maxerr rfft = %6.4e\n", maxerrfft);
free(a);
fft2dFree();
}
else{
if(TheErr==2) printf(" out of memory \n");
else printf(" error \n");
fft2dFree();
}
printf(" Done. \n");
return;
}
--- NEW FILE: convTest.c ---
/* A program to test fast 1d real convolution */
#include <stdio.h>
#include <stdlib.h>
#include <fp.h>
#include <math.h>
#include "fftlib.h"
#include "fftext.h"
#if macintosh
#include <timer.h>
#endif
#define BIPRAND(a) (2.0/(RAND_MAX+1.0)*a-1.0)
//#define BIPRAND(a) round(100*(2.0/(RAND_MAX+1.0)*a-1.0))
void main(){
const long N2 = 2; /* the number ffts to test */
long N = 2048; /* size of FFTs, must be power of 2 */
long kernSize = 1003; /* kernal size must be less than N */
long dataSize = N-kernSize+1; /* data size */
float *a;
float *b;
long i1;
long i2;
long TheErr;
long M;
FILE *fdataout; /* output file */
unsigned int randseed = 777;
int rannum;
#if macintosh
UnsignedWide TheTime1;
Microseconds(&TheTime1);
randseed = TheTime1.lo;
#endif
printf(" %6d Byte Floats \n", sizeof(a[0]));
printf(" randseed = %10u\n", randseed);
srand(randseed);
M = roundtol(LOG2(N));
N = POW2(M);
printf("fft size = %6d, ", N);
if (dataSize <= 0) TheErr = 22;
else TheErr = 0;
if(!TheErr){
TheErr = fftInit(M);
}
a = (float *) calloc(N2*N,sizeof(float) ); // calloc to zero pad data to
fill N
if (a == 0) TheErr = 2;
if(!TheErr){
b = (float *) calloc(N2*N,sizeof(float) ); // calloc to zero pad
data to fill N
if (b == 0) TheErr = 2;
}
if(!TheErr){
fdataout = fopen("convdat.cnv", "wb");
if (fdataout == NULL) TheErr = -50;
}
if(!TheErr){
/* write sizes to fdataout */
fwrite(&dataSize, sizeof(dataSize), 1, fdataout);
fwrite(&kernSize, sizeof(kernSize), 1, fdataout);
fwrite(&N2, sizeof(N2), 1, fdataout);
/* set up a simple test case and write to fdataout */
for (i2=0; i2<N2; i2++){
for (i1=0; i1<dataSize; i1++){
rannum = rand();
a[i2*N+i1] = BIPRAND(rannum);
}
fwrite(&a[i2*N], dataSize*sizeof(float), 1, fdataout);
}
for (i2=0; i2<N2; i2++){
for (i1=0; i1<kernSize; i1++){
rannum = rand();
b[i2*N+i1] = BIPRAND(rannum);
}
fwrite(&b[i2*N], kernSize*sizeof(float), 1, fdataout);
}
/* fast convolution of zero padded sequences */
rffts(a, M, N2);
rffts(b, M, N2);
for (i2=0; i2<N2*N; i2+=N){
rspectprod(&a[i2], &b[i2], &a[i2], N);
}
riffts(a, M, N2);
/* write out answer */
fwrite(a, N2*N*sizeof(float), 1, fdataout);
fclose(fdataout);
free(b);
free(a);
fftFree();
}
else{
if(TheErr==2) printf(" out of memory \n");
else printf(" error \n");
fftFree();
}
printf(" Done. \n");
return;
}
--- NEW FILE: convtest.m ---
% program to test 1d real fast conv
clear c2
% let user select file then open it
[fname, pname] = uigetfile('*.cnv', 'select conv file');
cd(pname);
fidout=fopen(fname,'r');
% read header info
aN=fread(fidout,1,'long');
bN=fread(fidout,1,'long');
M=fread(fidout,1,'long');
% read in data
%status=fseek(fidout,Nheader,'bof');
a=fread(fidout,aN*M,'float');
a=reshape(a,aN,M);
b=fread(fidout,bN*M,'float');
b=reshape(b,bN,M);
c=fread(fidout,(aN+bN-1)*M,'float');
c=reshape(c,(aN+bN-1),M);
fclose(fidout);
for i1=1:M;
c2(:,i1)=conv(a(:,i1),b(:,i1));
end;
max(max(abs(c2-c)))
------------------------------------------------------------------------------
Open Source Business Conference (OSBC), March 24-25, 2009, San Francisco, CA
-OSBC tackles the biggest issue in open source: Open Sourcing the Enterprise
-Strategies to boost innovation and cut costs with open source participation
-Receive a $600 discount off the registration fee with the source code: SFAD
http://p.sf.net/sfu/XcvMzF8H
_______________________________________________
Audacity-cvs mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/audacity-cvs
