luke
/
ngspice
mirror of https://git.code.sf.net/p/ngspice/ngspice
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

convert float to double

This commit is contained in:
h_vogt 2011-08-08 19:39:15 +00:00
parent 324b27a4d4
commit e7ce26c118
9 changed files with 3111 additions and 3064 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
ChangeLog
View File

@ -1,3 +1,8 @@
2011-08-08 Holger Vogt
* src/frontend/com_fft.c, src/maths/fft/fftext.c, fftext.h,
fftlib.c, fftlib.h, matlib.c, matlib.h: convert float to double
* winmain.c: increase text output buffer
2011-08-08 Robert Larice
* src/xspice/cm/cmutil.c :
ngspice.h always must be the first included file

40
src/frontend/com_fft.c
View File

@ -257,14 +257,14 @@ com_psd(wordlist *wl)
double **tdvec;
double *freq, *win, *time, *ave;
double delta_t, span, noipower;
int ngood, mm;
unsigned long fpts, i, j, tlen, jj, smooth, hsmooth;
int mm;
unsigned long size, ngood, fpts, i, j, tlen, jj, smooth, hsmooth;
char *s;
struct dvec *f, *vlist, *lv = NULL, *vec;
struct pnode *names, *first_name;
float *reald, *imagd;
int size, sign, isreal;
double *reald, *imagd;
int sign, isreal;
double scaling, sum;
int order;
double scale, sigma;
@ -403,7 +403,7 @@ com_psd(wordlist *wl)
vec = ft_evaluate(names);
names = names->pn_next;
while (vec) {
if (vec->v_length != tlen) {
if (vec->v_length != (int)tlen) {
fprintf(cp_err, "Error: lengths of %s vectors don't match: %d, %d\n",
vec->v_name, vec->v_length, tlen);
vec = vec->v_link2;
@ -428,7 +428,7 @@ com_psd(wordlist *wl)
ngood++;
}
}
free_pnode(first_name);
free_pnode_o(first_name);
if (!ngood) {
return;
}
@ -440,7 +440,7 @@ com_psd(wordlist *wl)
plot_cur->pl_name = copy("PSD");
plot_cur->pl_date = copy(datestring( ));
freq = (double *) tmalloc(fpts * sizeof(double));
freq = TMALLOC(double, fpts);
f = alloc(struct dvec);
ZERO(f, struct dvec);
f->v_name = copy("frequency");
@ -474,14 +474,14 @@ com_psd(wordlist *wl)
sign = 1;
isreal = 1;
reald = TMALLOC(float, size);
imagd = TMALLOC(float, size);
reald = TMALLOC(double, size);
imagd = TMALLOC(double, size);
// scale = 0.66;
for (i = 0; i<ngood; i++) {
for (j = 0; j < tlen; j++){
reald[j] = (float)(tdvec[i][j]*win[j]);
reald[j] = (tdvec[i][j]*win[j]);
imagd[j] = 0.;
}
for (j = tlen; j < size; j++){
@ -496,12 +496,12 @@ com_psd(wordlist *wl)
scaling = size*0.3;
/* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
noipower = fdvec[i][0].cx_real = (double)reald[0]*(double)reald[0];
fdvec[i][fpts-1].cx_real = (double)reald[1]*(double)reald[1];
noipower = fdvec[i][0].cx_real = reald[0]*reald[0];
fdvec[i][fpts-1].cx_real = reald[1]*reald[1];
noipower += fdvec[i][fpts-1].cx_real;
for (j=1; j<(fpts - 1); j++){
jj = j<<1;
fdvec[i][j].cx_real = ((double)reald[jj]*(double)reald[jj] + (double)reald[jj + 1]*(double)reald[jj + 1]);
fdvec[i][j].cx_real = reald[jj]*reald[jj] + reald[jj + 1]*reald[jj + 1];
fdvec[i][j].cx_imag = 0;
noipower += fdvec[i][j].cx_real;
}
@ -517,22 +517,22 @@ com_psd(wordlist *wl)
sum = 0.;
for (jj = 0; jj < hsmooth + j; jj++)
sum += fdvec[i][jj].cx_real;
sum /= (double)(hsmooth + j);
reald[j] = (float)(sqrt(sum)/scaling);
sum /= (hsmooth + j);
reald[j] = (sqrt(sum)/scaling);
}
for (j=hsmooth; j<fpts-hsmooth; j++){
sum = 0.;
for (jj = 0; jj < smooth; jj++)
sum += fdvec[i][j-hsmooth+jj].cx_real;
sum /= (double)smooth;
reald[j] = (float)(sqrt(sum)/scaling);
sum /= smooth;
reald[j] = (sqrt(sum)/scaling);
}
for (j=fpts-hsmooth; j<fpts; j++){
sum = 0.;
for (jj = 0; jj < smooth; jj++)
sum += fdvec[i][j-hsmooth+jj].cx_real;
sum /= (double)(fpts - j + hsmooth - 1);
reald[j] = (float)(sqrt(sum)/scaling);
sum /= (fpts - j + hsmooth - 1);
reald[j] = (sqrt(sum)/scaling);
}
for (j=0; j<fpts; j++)
fdvec[i][j].cx_real = reald[j];
@ -626,7 +626,7 @@ static void fftext(double* x, double* y, long int n, long int nn, int dir)
/* Scaling for forward transform */
if (dir == 1) {
double scale = 1.0 / (double) nn;
double scale = 1.0 / nn;
for (i=0;i<n;i++) {
x[i] *= scale; /* don't consider zero padded values */
y[i] *= scale;

239
src/maths/fft/fftext.c
View File

@ -1,8 +1,8 @@
/*******************************************************************
This file extends the fftlib with calls to maintain the cosine and bit reversed tables
for you (including mallocs and free's). Call the init routine for each fft size you will
be using. Then you can call the fft routines below which will make the fftlib library
call with the appropriate tables passed. When you are done with all fft's you can call
be using. Then you can call the fft routines below which will make the fftlib library
call with the appropriate tables passed. When you are done with all fft's you can call
fftfree to release the storage for the tables. Note that you can call fftinit repeatedly
with the same size, the extra calls will be ignored. So, you could make a macro to call
fftInit every time you call ffts. For example you could have someting like:
@ -14,69 +14,72 @@
#include "fftext.h"
// pointers to storage of Utbl's and BRLow's
static float *UtblArray[8*sizeof(long)] = {0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0};
static double *UtblArray[8*sizeof(long)] =
{0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0};
static short *BRLowArray[8*sizeof(long)/2] = {0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0};
int fftInit(long M){
int fftInit(long M)
{
// malloc and init cosine and bit reversed tables for a given size fft, ifft, rfft, rifft
/* INPUTS */
/* M = log2 of fft size (ex M=10 for 1024 point fft) */
/* OUTPUTS */
/* private cosine and bit reversed tables */
/* INPUTS */
/* M = log2 of fft size (ex M=10 for 1024 point fft) */
/* OUTPUTS */
/* private cosine and bit reversed tables */
int theError = 1;
/*** I did NOT test cases with M>27 ***/
if ((M >= 0) && (M < 8*sizeof(long))){
theError = 0;
if (UtblArray[M] == 0){ // have we not inited this size fft yet?
// init cos table
UtblArray[M] = (float *) malloc( (POW2(M)/4+1)*sizeof(float) );
if (UtblArray[M] == 0)
theError = 2;
else{
fftCosInit(M, UtblArray[M]);
}
if (M > 1){
if (BRLowArray[M/2] == 0){ // init bit reversed table for cmplx fft
BRLowArray[M/2] = (short *) malloc( POW2(M/2-1)*sizeof(short) );
if (BRLowArray[M/2] == 0)
theError = 2;
else{
fftBRInit(M, BRLowArray[M/2]);
}
}
}
if (M > 2){
if (BRLowArray[(M-1)/2] == 0){ // init bit reversed table for real fft
BRLowArray[(M-1)/2] = (short *) malloc( POW2((M-1)/2-1)*sizeof(short) );
if (BRLowArray[(M-1)/2] == 0)
theError = 2;
else{
fftBRInit(M-1, BRLowArray[(M-1)/2]);
}
}
}
}
};
return theError;
int theError = 1;
/*** I did NOT test cases with M>27 ***/
if ((M >= 0) && (M < 8*sizeof(long))) {
theError = 0;
if (UtblArray[M] == 0) { // have we not inited this size fft yet?
// init cos table
UtblArray[M] = (double *) malloc( (POW2(M)/4+1)*sizeof(double) );
if (UtblArray[M] == 0)
theError = 2;
else {
fftCosInit(M, UtblArray[M]);
}
if (M > 1) {
if (BRLowArray[M/2] == 0) { // init bit reversed table for cmplx fft
BRLowArray[M/2] = (short *) malloc( POW2(M/2-1)*sizeof(short) );
if (BRLowArray[M/2] == 0)
theError = 2;
else {
fftBRInit(M, BRLowArray[M/2]);
}
}
}
if (M > 2) {
if (BRLowArray[(M-1)/2] == 0) { // init bit reversed table for real fft
BRLowArray[(M-1)/2] = (short *) malloc( POW2((M-1)/2-1)*sizeof(short) );
if (BRLowArray[(M-1)/2] == 0)
theError = 2;
else {
fftBRInit(M-1, BRLowArray[(M-1)/2]);
}
}
}
}
};
return theError;
}
void fftFree(void){
void fftFree(void)
{
// release storage for all private cosine and bit reversed tables
long i1;
for (i1=8*sizeof(long)/2-1; i1>=0; i1--){
if (BRLowArray[i1] != 0){
free(BRLowArray[i1]);
BRLowArray[i1] = 0;
};
};
for (i1=8*sizeof(long)-1; i1>=0; i1--){
if (UtblArray[i1] != 0){
free(UtblArray[i1]);
UtblArray[i1] = 0;
};
};
long i1;
for (i1=8*sizeof(long)/2-1; i1>=0; i1--) {
if (BRLowArray[i1] != 0) {
free(BRLowArray[i1]);
BRLowArray[i1] = 0;
};
};
for (i1=8*sizeof(long)-1; i1>=0; i1--) {
if (UtblArray[i1] != 0) {
free(UtblArray[i1]);
UtblArray[i1] = 0;
};
};
}
/*************************************************
@ -84,73 +87,77 @@ for (i1=8*sizeof(long)-1; i1>=0; i1--){
Just make sure fftlib has been called for each M first.
**************************************************/
void ffts(float *data, long M, long Rows){
/* Compute in-place complex fft on the rows of the input array */
/* INPUTS */
/* *ioptr = input data array */
/* M = log2 of fft size (ex M=10 for 1024 point fft) */
/* Rows = number of rows in ioptr array (use 1 for Rows for a single fft) */
/* OUTPUTS */
/* *ioptr = output data array */
ffts1(data, M, Rows, UtblArray[M], BRLowArray[M/2]);
void ffts(double *data, long M, long Rows)
{
/* Compute in-place complex fft on the rows of the input array */
/* INPUTS */
/* *ioptr = input data array */
/* M = log2 of fft size (ex M=10 for 1024 point fft) */
/* Rows = number of rows in ioptr array (use 1 for Rows for a single fft) */
/* OUTPUTS */
/* *ioptr = output data array */
ffts1(data, M, Rows, UtblArray[M], BRLowArray[M/2]);
}
void iffts(float *data, long M, long Rows){
/* Compute in-place inverse complex fft on the rows of the input array */
/* INPUTS */
/* *ioptr = input data array */
/* M = log2 of fft size (ex M=10 for 1024 point fft) */
/* Rows = number of rows in ioptr array (use 1 for Rows for a single fft) */
/* OUTPUTS */
/* *ioptr = output data array */
iffts1(data, M, Rows, UtblArray[M], BRLowArray[M/2]);
void iffts(double *data, long M, long Rows)
{
/* Compute in-place inverse complex fft on the rows of the input array */
/* INPUTS */
/* *ioptr = input data array */
/* M = log2 of fft size (ex M=10 for 1024 point fft) */
/* Rows = number of rows in ioptr array (use 1 for Rows for a single fft) */
/* OUTPUTS */
/* *ioptr = output data array */
iffts1(data, M, Rows, UtblArray[M], BRLowArray[M/2]);
}
void rffts(float *data, long M, long Rows){
/* Compute in-place real fft on the rows of the input array */
/* The result is the complex spectra of the positive frequencies */
/* except the location for the first complex number contains the real */
/* values for DC and Nyquest */
/* See rspectprod for multiplying two of these spectra together- ex. for fast convolution */
/* INPUTS */
/* *ioptr = real input data array */
/* M = log2 of fft size (ex M=10 for 1024 point fft) */
/* Rows = number of rows in ioptr array (use 1 for Rows for a single fft) */
/* OUTPUTS */
/* *ioptr = output data array in the following order */
/* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
rffts1(data, M, Rows, UtblArray[M], BRLowArray[(M-1)/2]);
void rffts(double *data, long M, long Rows)
{
/* Compute in-place real fft on the rows of the input array */
/* The result is the complex spectra of the positive frequencies */
/* except the location for the first complex number contains the real */
/* values for DC and Nyquest */
/* See rspectprod for multiplying two of these spectra together- ex. for fast convolution */
/* INPUTS */
/* *ioptr = real input data array */
/* M = log2 of fft size (ex M=10 for 1024 point fft) */
/* Rows = number of rows in ioptr array (use 1 for Rows for a single fft) */
/* OUTPUTS */
/* *ioptr = output data array in the following order */
/* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
rffts1(data, M, Rows, UtblArray[M], BRLowArray[(M-1)/2]);
}
void riffts(float *data, long M, long Rows){
/* Compute in-place real ifft on the rows of the input array */
/* data order as from rffts */
/* INPUTS */
/* *ioptr = input data array in the following order */
/* M = log2 of fft size (ex M=10 for 1024 point fft) */
/* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
/* Rows = number of rows in ioptr array (use 1 for Rows for a single fft) */
/* OUTPUTS */
/* *ioptr = real output data array */
riffts1(data, M, Rows, UtblArray[M], BRLowArray[(M-1)/2]);
void riffts(double *data, long M, long Rows)
{
/* Compute in-place real ifft on the rows of the input array */
/* data order as from rffts */
/* INPUTS */
/* *ioptr = input data array in the following order */
/* M = log2 of fft size (ex M=10 for 1024 point fft) */
/* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
/* Rows = number of rows in ioptr array (use 1 for Rows for a single fft) */
/* OUTPUTS */
/* *ioptr = real output data array */
riffts1(data, M, Rows, UtblArray[M], BRLowArray[(M-1)/2]);
}
void rspectprod(float *data1, float *data2, float *outdata, long N){
void rspectprod(double *data1, double *data2, double *outdata, long N)
{
// When multiplying a pair of spectra from rfft care must be taken to multiply the
// two real values seperately from the complex ones. This routine does it correctly.
// the result can be stored in-place over one of the inputs
/* INPUTS */
/* *data1 = input data array first spectra */
/* *data2 = input data array second spectra */
/* N = fft input size for both data1 and data2 */
/* OUTPUTS */
/* *outdata = output data array spectra */
if(N>1){
outdata[0] = data1[0] * data2[0]; // multiply the zero freq values
outdata[1] = data1[1] * data2[1]; // multiply the nyquest freq values
cvprod(data1 + 2, data2 + 2, outdata + 2, N/2-1); // multiply the other positive freq values
}
else{
outdata[0] = data1[0] * data2[0];
}
/* INPUTS */
/* *data1 = input data array first spectra */
/* *data2 = input data array second spectra */
/* N = fft input size for both data1 and data2 */
/* OUTPUTS */
/* *outdata = output data array spectra */
if(N>1) {
outdata[0] = data1[0] * data2[0]; // multiply the zero freq values
outdata[1] = data1[1] * data2[1]; // multiply the nyquest freq values
cvprod(data1 + 2, data2 + 2, outdata + 2, N/2-1); // multiply the other positive freq values
} else {
outdata[0] = data1[0] * data2[0];
}
}

14
src/maths/fft/fftext.h
View File

@ -19,7 +19,7 @@ int fftInit(long M);
void fftFree(void);
// release storage for all private cosine and bit reversed tables
void ffts(float *data, long M, long Rows);
void ffts(double *data, long M, long Rows);
/* Compute in-place complex fft on the rows of the input array */
/* INPUTS */
/* *ioptr = input data array */
@ -28,7 +28,7 @@ void ffts(float *data, long M, long Rows);
/* OUTPUTS */
/* *ioptr = output data array */
void iffts(float *data, long M, long Rows);
void iffts(double *data, long M, long Rows);
/* Compute in-place inverse complex fft on the rows of the input array */
/* INPUTS */
/* *ioptr = input data array */
@ -37,7 +37,7 @@ void iffts(float *data, long M, long Rows);
/* OUTPUTS */
/* *ioptr = output data array */
void rffts(float *data, long M, long Rows);
void rffts(double *data, long M, long Rows);
/* Compute in-place real fft on the rows of the input array */
/* The result is the complex spectra of the positive frequencies */
/* except the location for the first complex number contains the real */
@ -51,7 +51,7 @@ void rffts(float *data, long M, long Rows);
/* *ioptr = output data array in the following order */
/* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
void riffts(float *data, long M, long Rows);
void riffts(double *data, long M, long Rows);
/* Compute in-place real ifft on the rows of the input array */
/* data order as from rffts */
/* INPUTS */
@ -62,7 +62,7 @@ void riffts(float *data, long M, long Rows);
/* OUTPUTS */
/* *ioptr = real output data array */
void rspectprod(float *data1, float *data2, float *outdata, long N);
void rspectprod(double *data1, double *data2, double *outdata, long N);
// When multiplying a pair of spectra from rfft care must be taken to multiply the
// two real values seperately from the complex ones. This routine does it correctly.
// the result can be stored in-place over one of the inputs
@ -80,10 +80,10 @@ void rspectprod(float *data1, float *data2, float *outdata, long N);
//Note that most of the fft routines require full matrices, ie Rsiz==Ncols
//This is how I like to define a real matrix:
//struct matrix { // real matrix
// float *d; // pointer to data
// double *d; // pointer to data
// long Nrows; // number of rows in the matrix
// long Ncols; // number of columns in the matrix (can be less than Rsiz)
// long Rsiz; // number of floats from one row to the next
// long Rsiz; // number of doubles from one row to the next
//};
//typedef struct matrix matrix;

5288
src/maths/fft/fftlib.c
View File

File diff suppressed because it is too large Load Diff

12
src/maths/fft/fftlib.h
View File

@ -1,14 +1,14 @@
#define MYRECIPLN2 1.442695040888963407359924681001892137426 // 1.0/log(2)
/* some useful conversions between a number and its power of 2 */
#define LOG2(a) (MYRECIPLN2*log(a)) // floating point logarithm base 2
#define LOG2(a) (MYRECIPLN2*log(a)) // doubleing point logarithm base 2
#define POW2(m) ((unsigned long) 1 << (m)) // integer power of 2 for m<32
/*******************************************************************
lower level fft stuff called by routines in fftext.c and fft2d.c
*******************************************************************/
void fftCosInit(long M, float *Utbl);
void fftCosInit(long M, double *Utbl);
/* Compute Utbl, the cosine table for ffts */
/* of size (pow(2,M)/4 +1) */
/* INPUTS */
@ -24,7 +24,7 @@ void fftBRInit(long M, short *BRLow);
/* OUTPUTS */
/* *BRLow = bit reversed counter table */
void ffts1(float *ioptr, long M, long Rows, float *Utbl, short *BRLow);
void ffts1(double *ioptr, long M, long Rows, double *Utbl, short *BRLow);
/* Compute in-place complex fft on the rows of the input array */
/* INPUTS */
/* *ioptr = input data array */
@ -35,7 +35,7 @@ void ffts1(float *ioptr, long M, long Rows, float *Utbl, short *BRLow);
/* OUTPUTS */
/* *ioptr = output data array */
void iffts1(float *ioptr, long M, long Rows, float *Utbl, short *BRLow);
void iffts1(double *ioptr, long M, long Rows, double *Utbl, short *BRLow);
/* Compute in-place inverse complex fft on the rows of the input array */
/* INPUTS */
/* *ioptr = input data array */
@ -46,7 +46,7 @@ void iffts1(float *ioptr, long M, long Rows, float *Utbl, short *BRLow);
/* OUTPUTS */
/* *ioptr = output data array */
void rffts1(float *ioptr, long M, long Rows, float *Utbl, short *BRLow);
void rffts1(double *ioptr, long M, long Rows, double *Utbl, short *BRLow);
/* Compute in-place real fft on the rows of the input array */
/* The result is the complex spectra of the positive frequencies */
/* except the location for the first complex number contains the real */
@ -62,7 +62,7 @@ void rffts1(float *ioptr, long M, long Rows, float *Utbl, short *BRLow);
/* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
void riffts1(float *ioptr, long M, long Rows, float *Utbl, short *BRLow);
void riffts1(double *ioptr, long M, long Rows, double *Utbl, short *BRLow);
/* Compute in-place real ifft on the rows of the input array */
/* data order as from rffts1 */
/* INPUTS */

569
src/maths/fft/matlib.c
View File

@ -1,297 +1,300 @@
/* a few routines from a vector/matrix library */
#include "matlib.h"
void xpose(float *indata, long iRsiz, float *outdata, long oRsiz, long Nrows, long Ncols){
/* not in-place matrix transpose */
/* INPUTS */
/* *indata = input data array */
/* iRsiz = offset to between rows of input data array */
/* oRsiz = offset to between rows of output data array */
/* Nrows = number of rows in input data array */
/* Ncols = number of columns in input data array */
/* OUTPUTS */
/* *outdata = output data array */
void xpose(double *indata, long iRsiz, double *outdata, long oRsiz, long Nrows, long Ncols)
{
/* not in-place matrix transpose */
/* INPUTS */
/* *indata = input data array */
/* iRsiz = offset to between rows of input data array */
/* oRsiz = offset to between rows of output data array */
/* Nrows = number of rows in input data array */
/* Ncols = number of columns in input data array */
/* OUTPUTS */
/* *outdata = output data array */
float *irow; /* pointer to input row start */
float *ocol; /* pointer to output col start */
float *idata; /* pointer to input data */
float *odata; /* pointer to output data */
long RowCnt; /* row counter */
long ColCnt; /* col counter */
float T0; /* data storage */
float T1; /* data storage */
float T2; /* data storage */
float T3; /* data storage */
float T4; /* data storage */
float T5; /* data storage */
float T6; /* data storage */
float T7; /* data storage */
const long inRsizd1 = iRsiz;
const long inRsizd2 = 2*iRsiz;
const long inRsizd3 = inRsizd2+iRsiz;
const long inRsizd4 = 4*iRsiz;
const long inRsizd5 = inRsizd3+inRsizd2;
const long inRsizd6 = inRsizd4+inRsizd2;
const long inRsizd7 = inRsizd4+inRsizd3;
const long inRsizd8 = 8*iRsiz;
double *irow; /* pointer to input row start */
double *ocol; /* pointer to output col start */
double *idata; /* pointer to input data */
double *odata; /* pointer to output data */
long RowCnt; /* row counter */
long ColCnt; /* col counter */
double T0; /* data storage */
double T1; /* data storage */
double T2; /* data storage */
double T3; /* data storage */
double T4; /* data storage */
double T5; /* data storage */
double T6; /* data storage */
double T7; /* data storage */
const long inRsizd1 = iRsiz;
const long inRsizd2 = 2*iRsiz;
const long inRsizd3 = inRsizd2+iRsiz;
const long inRsizd4 = 4*iRsiz;
const long inRsizd5 = inRsizd3+inRsizd2;
const long inRsizd6 = inRsizd4+inRsizd2;
const long inRsizd7 = inRsizd4+inRsizd3;
const long inRsizd8 = 8*iRsiz;
ocol = outdata;
irow = indata;
for (RowCnt=Nrows/8; RowCnt>0; RowCnt--){
idata = irow;
odata = ocol;
for (ColCnt=Ncols; ColCnt>0; ColCnt--){
T0 = *idata;
T1 = *(idata+inRsizd1);
T2 = *(idata+inRsizd2);
T3 = *(idata+inRsizd3);
T4 = *(idata+inRsizd4);
T5 = *(idata+inRsizd5);
T6 = *(idata+inRsizd6);
T7 = *(idata+inRsizd7);
*odata = T0;
*(odata+1) = T1;
*(odata+2) = T2;
*(odata+3) = T3;
*(odata+4) = T4;
*(odata+5) = T5;
*(odata+6) = T6;
*(odata+7) = T7;
idata++;
odata += oRsiz;
}
irow += inRsizd8;
ocol += 8;
}
if (Nrows%8 != 0){
for (ColCnt=Ncols; ColCnt>0; ColCnt--){
idata = irow++;
odata = ocol;
ocol += oRsiz;
for (RowCnt=Nrows%8; RowCnt>0; RowCnt--){
T0 = *idata;
*odata++ = T0;
idata += iRsiz;
}
}
}
ocol = outdata;
irow = indata;
for (RowCnt=Nrows/8; RowCnt>0; RowCnt--) {
idata = irow;
odata = ocol;
for (ColCnt=Ncols; ColCnt>0; ColCnt--) {
T0 = *idata;
T1 = *(idata+inRsizd1);
T2 = *(idata+inRsizd2);
T3 = *(idata+inRsizd3);
T4 = *(idata+inRsizd4);
T5 = *(idata+inRsizd5);
T6 = *(idata+inRsizd6);
T7 = *(idata+inRsizd7);
*odata = T0;
*(odata+1) = T1;
*(odata+2) = T2;
*(odata+3) = T3;
*(odata+4) = T4;
*(odata+5) = T5;
*(odata+6) = T6;
*(odata+7) = T7;
idata++;
odata += oRsiz;
}
irow += inRsizd8;
ocol += 8;
}
if (Nrows%8 != 0) {
for (ColCnt=Ncols; ColCnt>0; ColCnt--) {
idata = irow++;
odata = ocol;
ocol += oRsiz;
for (RowCnt=Nrows%8; RowCnt>0; RowCnt--) {
T0 = *idata;
*odata++ = T0;
idata += iRsiz;
}
}
}
}
void cxpose(float *indata, long iRsiz, float *outdata, long oRsiz, long Nrows, long Ncols){
/* not in-place complex float matrix transpose */
/* INPUTS */
/* *indata = input data array */
/* iRsiz = offset to between rows of input data array */
/* oRsiz = offset to between rows of output data array */
/* Nrows = number of rows in input data array */
/* Ncols = number of columns in input data array */
/* OUTPUTS */
/* *outdata = output data array */
void cxpose(double *indata, long iRsiz, double *outdata, long oRsiz, long Nrows, long Ncols)
{
/* not in-place complex double matrix transpose */
/* INPUTS */
/* *indata = input data array */
/* iRsiz = offset to between rows of input data array */
/* oRsiz = offset to between rows of output data array */
/* Nrows = number of rows in input data array */
/* Ncols = number of columns in input data array */
/* OUTPUTS */
/* *outdata = output data array */
float *irow; /* pointer to input row start */
float *ocol; /* pointer to output col start */
float *idata; /* pointer to input data */
float *odata; /* pointer to output data */
long RowCnt; /* row counter */
long ColCnt; /* col counter */
float T0r; /* data storage */
float T0i; /* data storage */
float T1r; /* data storage */
float T1i; /* data storage */
float T2r; /* data storage */
float T2i; /* data storage */
float T3r; /* data storage */
float T3i; /* data storage */
const long inRsizd1 = 2*iRsiz;
const long inRsizd1i = 2*iRsiz + 1;
const long inRsizd2 = 4*iRsiz;
const long inRsizd2i = 4*iRsiz + 1;
const long inRsizd3 = inRsizd2+inRsizd1;
const long inRsizd3i = inRsizd2+inRsizd1 + 1;
const long inRsizd4 = 8*iRsiz;
double *irow; /* pointer to input row start */
double *ocol; /* pointer to output col start */
double *idata; /* pointer to input data */
double *odata; /* pointer to output data */
long RowCnt; /* row counter */
long ColCnt; /* col counter */
double T0r; /* data storage */
double T0i; /* data storage */
double T1r; /* data storage */
double T1i; /* data storage */
double T2r; /* data storage */
double T2i; /* data storage */
double T3r; /* data storage */
double T3i; /* data storage */
const long inRsizd1 = 2*iRsiz;
const long inRsizd1i = 2*iRsiz + 1;
const long inRsizd2 = 4*iRsiz;
const long inRsizd2i = 4*iRsiz + 1;
const long inRsizd3 = inRsizd2+inRsizd1;
const long inRsizd3i = inRsizd2+inRsizd1 + 1;
const long inRsizd4 = 8*iRsiz;
ocol = outdata;
irow = indata;
for (RowCnt=Nrows/4; RowCnt>0; RowCnt--){
idata = irow;
odata = ocol;
for (ColCnt=Ncols; ColCnt>0; ColCnt--){
T0r = *idata;
T0i = *(idata +1);
T1r = *(idata+inRsizd1);
T1i = *(idata+inRsizd1i);
T2r = *(idata+inRsizd2);
T2i = *(idata+inRsizd2i);
T3r = *(idata+inRsizd3);
T3i = *(idata+inRsizd3i);
*odata = T0r;
*(odata+1) = T0i;
*(odata+2) = T1r;
*(odata+3) = T1i;
*(odata+4) = T2r;
*(odata+5) = T2i;
*(odata+6) = T3r;
*(odata+7) = T3i;
idata+=2;
odata += 2*oRsiz;
}
irow += inRsizd4;
ocol += 8;
}
if (Nrows%4 != 0){
for (ColCnt=Ncols; ColCnt>0; ColCnt--){
idata = irow;
odata = ocol;
for (RowCnt=Nrows%4; RowCnt>0; RowCnt--){
T0r = *idata;
T0i = *(idata+1);
*odata = T0r;
*(odata+1) = T0i;
odata+=2;
idata += 2*iRsiz;
}
irow+=2;
ocol += 2*oRsiz;
}
}
ocol = outdata;
irow = indata;
for (RowCnt=Nrows/4; RowCnt>0; RowCnt--) {
idata = irow;
odata = ocol;
for (ColCnt=Ncols; ColCnt>0; ColCnt--) {
T0r = *idata;
T0i = *(idata +1);
T1r = *(idata+inRsizd1);
T1i = *(idata+inRsizd1i);
T2r = *(idata+inRsizd2);
T2i = *(idata+inRsizd2i);
T3r = *(idata+inRsizd3);
T3i = *(idata+inRsizd3i);
*odata = T0r;
*(odata+1) = T0i;
*(odata+2) = T1r;
*(odata+3) = T1i;
*(odata+4) = T2r;
*(odata+5) = T2i;
*(odata+6) = T3r;
*(odata+7) = T3i;
idata+=2;
odata += 2*oRsiz;
}
irow += inRsizd4;
ocol += 8;
}
if (Nrows%4 != 0) {
for (ColCnt=Ncols; ColCnt>0; ColCnt--) {
idata = irow;
odata = ocol;
for (RowCnt=Nrows%4; RowCnt>0; RowCnt--) {
T0r = *idata;
T0i = *(idata+1);
*odata = T0r;
*(odata+1) = T0i;
odata+=2;
idata += 2*iRsiz;
}
irow+=2;
ocol += 2*oRsiz;
}
}
}
void cvprod(float *a, float *b, float *out, long N){
/* complex vector product, can be in-place */
/* product of complex vector *a times complex vector *b */
/* INPUTS */
/* N vector length */
/* *a complex vector length N complex numbers */
/* *b complex vector length N complex numbers */
/* OUTPUTS */
/* *out complex vector length N */
void cvprod(double *a, double *b, double *out, long N)
{
/* complex vector product, can be in-place */
/* product of complex vector *a times complex vector *b */
/* INPUTS */
/* N vector length */
/* *a complex vector length N complex numbers */
/* *b complex vector length N complex numbers */
/* OUTPUTS */
/* *out complex vector length N */
long OutCnt; /* output counter */
float A0R; /* A storage */
float A0I; /* A storage */
float A1R; /* A storage */
float A1I; /* A storage */
float A2R; /* A storage */
float A2I; /* A storage */
float A3R; /* A storage */
float A3I; /* A storage */
float B0R; /* B storage */
float B0I; /* B storage */
float B1R; /* B storage */
float B1I; /* B storage */
float B2R; /* B storage */
float B2I; /* B storage */
float B3R; /* B storage */
float B3I; /* B storage */
float T0R; /* TMP storage */
float T0I; /* TMP storage */
float T1R; /* TMP storage */
float T1I; /* TMP storage */
float T2R; /* TMP storage */
float T2I; /* TMP storage */
float T3R; /* TMP storage */
float T3I; /* TMP storage */
long OutCnt; /* output counter */
double A0R; /* A storage */
double A0I; /* A storage */
double A1R; /* A storage */
double A1I; /* A storage */
double A2R; /* A storage */
double A2I; /* A storage */
double A3R; /* A storage */
double A3I; /* A storage */
double B0R; /* B storage */
double B0I; /* B storage */
double B1R; /* B storage */
double B1I; /* B storage */
double B2R; /* B storage */
double B2I; /* B storage */
double B3R; /* B storage */
double B3I; /* B storage */
double T0R; /* TMP storage */
double T0I; /* TMP storage */
double T1R; /* TMP storage */
double T1I; /* TMP storage */
double T2R; /* TMP storage */
double T2I; /* TMP storage */
double T3R; /* TMP storage */
double T3I; /* TMP storage */
if (N>=4){
A0R = *a;
B0R = *b;
A0I = *(a +1);
B0I = *(b +1);
A1R = *(a +2);
B1R = *(b +2);
A1I = *(a +3);
B1I = *(b +3);
A2R = *(a +4);
B2R = *(b +4);
A2I = *(a +5);
B2I = *(b +5);
A3R = *(a +6);
B3R = *(b +6);
A3I = *(a +7);
B3I = *(b +7);
T0R = A0R * B0R;
T0I = (A0R * B0I);
T1R = A1R * B1R;
T1I = (A1R * B1I);
T2R = A2R * B2R;
T2I = (A2R * B2I);
T3R = A3R * B3R;
T3I = (A3R * B3I);
T0R -= (A0I * B0I);
T0I = A0I * B0R + T0I;
T1R -= (A1I * B1I);
T1I = A1I * B1R + T1I;
T2R -= (A2I * B2I);
T2I = A2I * B2R + T2I;
T3R -= (A3I * B3I);
T3I = A3I * B3R + T3I;
for (OutCnt=N/4-1; OutCnt > 0; OutCnt--){
a += 8;
b += 8;
A0R = *a;
B0R = *b;
A0I = *(a +1);
B0I = *(b +1);
A1R = *(a +2);
B1R = *(b +2);
A1I = *(a +3);
B1I = *(b +3);
A2R = *(a +4);
B2R = *(b +4);
A2I = *(a +5);
B2I = *(b +5);
A3R = *(a +6);
B3R = *(b +6);
A3I = *(a +7);
B3I = *(b +7);
*out = T0R;
*(out +1) = T0I;
*(out +2) = T1R;
*(out +3) = T1I;
*(out +4) = T2R;
*(out +5) = T2I;
*(out +6) = T3R;
*(out +7) = T3I;
T0R = A0R * B0R;
T0I = (A0R * B0I);
T1R = A1R * B1R;
T1I = (A1R * B1I);
T2R = A2R * B2R;
T2I = (A2R * B2I);
T3R = A3R * B3R;
T3I = (A3R * B3I);
T0R -= (A0I * B0I);
T0I = A0I * B0R + T0I;
T1R -= (A1I * B1I);
T1I = A1I * B1R + T1I;
T2R -= (A2I * B2I);
T2I = A2I * B2R + T2I;
T3R -= (A3I * B3I);
T3I = A3I * B3R + T3I;
out += 8;
}
a += 8;
b += 8;
*out = T0R;
*(out +1) = T0I;
*(out +2) = T1R;
*(out +3) = T1I;
*(out +4) = T2R;
*(out +5) = T2I;
*(out +6) = T3R;
*(out +7) = T3I;
out += 8;
}
for (OutCnt=N%4; OutCnt > 0; OutCnt--){
A0R = *a++;
B0R = *b++;
A0I = *a++;
B0I = *b++;
T0R = A0R * B0R;
T0I = (A0R * B0I);
T0R -= (A0I * B0I);
T0I = A0I * B0R + T0I;
*out++ = T0R;
*out++ = T0I;
}
if (N>=4) {
A0R = *a;
B0R = *b;
A0I = *(a +1);
B0I = *(b +1);
A1R = *(a +2);
B1R = *(b +2);
A1I = *(a +3);
B1I = *(b +3);
A2R = *(a +4);
B2R = *(b +4);
A2I = *(a +5);
B2I = *(b +5);
A3R = *(a +6);
B3R = *(b +6);
A3I = *(a +7);
B3I = *(b +7);
T0R = A0R * B0R;
T0I = (A0R * B0I);
T1R = A1R * B1R;
T1I = (A1R * B1I);
T2R = A2R * B2R;
T2I = (A2R * B2I);
T3R = A3R * B3R;
T3I = (A3R * B3I);
T0R -= (A0I * B0I);
T0I = A0I * B0R + T0I;
T1R -= (A1I * B1I);
T1I = A1I * B1R + T1I;
T2R -= (A2I * B2I);
T2I = A2I * B2R + T2I;
T3R -= (A3I * B3I);
T3I = A3I * B3R + T3I;
for (OutCnt=N/4-1; OutCnt > 0; OutCnt--) {
a += 8;
b += 8;
A0R = *a;
B0R = *b;
A0I = *(a +1);
B0I = *(b +1);
A1R = *(a +2);
B1R = *(b +2);
A1I = *(a +3);
B1I = *(b +3);
A2R = *(a +4);
B2R = *(b +4);
A2I = *(a +5);
B2I = *(b +5);
A3R = *(a +6);
B3R = *(b +6);
A3I = *(a +7);
B3I = *(b +7);
*out = T0R;
*(out +1) = T0I;
*(out +2) = T1R;
*(out +3) = T1I;
*(out +4) = T2R;
*(out +5) = T2I;
*(out +6) = T3R;
*(out +7) = T3I;
T0R = A0R * B0R;
T0I = (A0R * B0I);
T1R = A1R * B1R;
T1I = (A1R * B1I);
T2R = A2R * B2R;
T2I = (A2R * B2I);
T3R = A3R * B3R;
T3I = (A3R * B3I);
T0R -= (A0I * B0I);
T0I = A0I * B0R + T0I;
T1R -= (A1I * B1I);
T1I = A1I * B1R + T1I;
T2R -= (A2I * B2I);
T2I = A2I * B2R + T2I;
T3R -= (A3I * B3I);
T3I = A3I * B3R + T3I;
out += 8;
}
a += 8;
b += 8;
*out = T0R;
*(out +1) = T0I;
*(out +2) = T1R;
*(out +3) = T1I;
*(out +4) = T2R;
*(out +5) = T2I;
*(out +6) = T3R;
*(out +7) = T3I;
out += 8;
}
for (OutCnt=N%4; OutCnt > 0; OutCnt--) {
A0R = *a++;
B0R = *b++;
A0I = *a++;
B0I = *b++;
T0R = A0R * B0R;
T0I = (A0R * B0I);
T0R -= (A0I * B0I);
T0I = A0I * B0R + T0I;
*out++ = T0R;
*out++ = T0I;
}
}

6
src/maths/fft/matlib.h
View File

@ -1,6 +1,6 @@
/* a few routines from a vector/matrix library */
void xpose(float *indata, long iRsiz, float *outdata, long oRsiz, long Nrows, long Ncols);
void xpose(double *indata, long iRsiz, double *outdata, long oRsiz, long Nrows, long Ncols);
/* not in-place matrix transpose */
/* INPUTS */
/* *indata = input data array */
@ -11,7 +11,7 @@ void xpose(float *indata, long iRsiz, float *outdata, long oRsiz, long Nrows, lo
/* OUTPUTS */
/* *outdata = output data array */
void cxpose(float *indata, long iRsiz, float *outdata, long oRsiz, long Nrows, long Ncols);
void cxpose(double *indata, long iRsiz, double *outdata, long oRsiz, long Nrows, long Ncols);
/* not in-place complex matrix transpose */
/* INPUTS */
/* *indata = input data array */
@ -22,7 +22,7 @@ void cxpose(float *indata, long iRsiz, float *outdata, long oRsiz, long Nrows, l
/* OUTPUTS */
/* *outdata = output data array */
void cvprod(float *a, float *b, float *out, long N);
void cvprod(double *a, double *b, double *out, long N);
/* complex vector product, can be in-place */
/* product of complex vector *a times complex vector *b */
/* INPUTS */

2
src/winmain.c
View File

@ -29,7 +29,7 @@
#include "misc/misc_time.h" /* timediff */
/* Constants */
#define TBufSize 8192 // size of text buffer
#define TBufSize 65536 // size of text buffer
#define CR VK_RETURN // Carriage Return
#define LF 10 // Line Feed
#define SE 0 // String termination