This is gck_fft.c in view mode; [Download] [Up]
/* ******************************************************************** */
/* * GCK * */
/* * A Circuit Simulation Program * */
/* * by Tanju Cataltepe * */
/* * (c) Copyright 1989 * */
/* ******************************************************************** */
/* (c) Copyright 1989, Tanju Cataltepe */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <time.h>
#ifdef Macintosh
#include <console.h>
#include <Files.h>
#else
#include <sys/types.h>
#include <sys/stat.h>
#endif
#include <math.h>
#include <float.h>
#include "gck.h"
#include "gck_vars.h"
#include "gck_protos.h"
#include "gck_errors.h"
fftMRec *get_fft_mrec(fftRec *tmpP)
{
fftMRec *tmpM;
fftMRec *newptr;
tmpM = baseOfMFFT;
while (tmpM != NULL) {
if ((tmpP->vo1 == tmpM->vo1) && (tmpP->vo2 == tmpM->vo2))
break;
tmpM = tmpM->next;
}
if (tmpM != NULL)
return(tmpM);
else {
newptr = (fftMRec *) c_malloc(sizeof(fftMRec));
newptr->vo1 = tmpP->vo1;
newptr->vo2 = tmpP->vo2;
newptr->next = baseOfMFFT;
baseOfMFFT = newptr;
newptr->fft_length = 0;
newptr->base = (fftClist *) c_malloc(sizeof(fftClist));
newptr->current = newptr->base;
newptr->current->depth = 0;
newptr->current->next = NULL;
return(newptr);
}
}
void readFFT()
{
/*.FFT [Window windowname] V(n1, n2) I(name) ... [> filename] */
fftRec *tmpP, *lastP;
aWord tmpW;
nameString tmpN;
long i, j, k, m;
long counter;
i = 2;
GetWORD(tmpW,i);
while (tmpW[0] != NULL_CHAR) {
switch (tmpW[0]) {
case 'W':
i++;
GetWORD(tmpW,i);
if (strncmp(tmpW,"RECT",4) == 0)
fft_window = WINDOW_RECTANGULAR;
else if (strncmp(tmpW,"BART",4) == 0)
fft_window = WINDOW_TRIANGULAR;
else if (strncmp(tmpW,"TRIA",4) == 0)
fft_window = WINDOW_TRIANGULAR;
else if (strncmp(tmpW,"HANN",4) == 0)
fft_window = WINDOW_HANN;
else if (strncmp(tmpW,"HAMM",4) == 0)
fft_window = WINDOW_HAMMING;
else if (strncmp(tmpW,"BLAC",4) == 0)
fft_window = WINDOW_BLACKMAN;
else
errorReport(ERROR_WINDOW_NOT_FOUND,tmpW,ERR_VOID,0);
break;
case '>': /* new option */
fft_to_file = TRUE;
if (tmpW[1] == NULL_CHAR) {
i++;
GetWord(fft_filename,i);
}
else {
GetWord(tmpW,i);
strcpy(fft_filename,&tmpW[1]);
}
break;
/* Windowing On Another Card ! */
case 'V':
tmpP = (fftRec *)c_malloc(sizeof(fftRec));
tmpP->outType = branchVoltage;
strcpy(tmpN,"\0");
if (tmpW[1] != '(') {
m = 1;
switch(tmpW[1]) {
case 'M': tmpP->fptype = FP_MAGN; break;
case 'R': tmpP->fptype = FP_REAL; break;
case 'I': tmpP->fptype = FP_IMAG; break;
case 'P': tmpP->fptype = FP_PHASE; break;
case 'D': tmpP->fptype = FP_DB; m = 2; break;
}
}
else {
tmpP->fptype = FP_MAGN;
m = 0;
}
j = 3;
while ((tmpW[j - 1 + m] != ',') && (tmpW[j - 1 + m] != ')')) {
tmpN[j - 3] = tmpW[j - 1 + m];
j++;
}
tmpN[j - 3] = NULL_CHAR;
tmpP->n1 = nodeRecordFromList(tmpN,REQUEST_NODE);
if (tmpW[j - 1 + m] == ')') /* Remembering Ayse! */
tmpP->n2 = nodeRecordFromList("0",REQUEST_NODE);
else {
strcpy(tmpN,"\0");
k = j + 1;
while (tmpW[k - 1 + m] != ')') {
tmpN[k - j - 1] = tmpW[k - 1 + m];
k++;
}
tmpN[k - j - 1] = NULL_CHAR;
tmpP->n2 = nodeRecordFromList(tmpN,REQUEST_NODE);
}
if (baseOfFFT == NULL) {
baseOfFFT = tmpP;
lastP = tmpP;
} else {
lastP->next = tmpP;
lastP = tmpP;
}
break;
case 'I':
tmpP = (fftRec *)c_malloc(sizeof(fftRec));
tmpP->outType = branchCurrent;
strcpy(tmpN,"\0");
if (tmpW[1] != '(') {
m = 1;
switch(toupper(tmpW[1])) {
case 'M': tmpP->fptype = FP_MAGN; break;
case 'R': tmpP->fptype = FP_REAL; break;
case 'I': tmpP->fptype = FP_IMAG; break;
case 'P': tmpP->fptype = FP_PHASE; break;
case 'D': tmpP->fptype = FP_DB; m = 2; break;
}
}
else {
tmpP->fptype = FP_MAGN;
m = 0;
}
j = 3;
while (tmpW[j - 1 + m] != ')') {
tmpN[j - 3] = tmpW[j - 1 + m];
j++;
}
tmpN[j - 3] = NULL_CHAR;
strcpy(tmpP->name,tmpN);
if (baseOfFFT == NULL) {
baseOfFFT = tmpP;
lastP = tmpP;
} else {
lastP->next = tmpP;
lastP = tmpP;
}
break;
default:
errorReport(ERROR_BAD_FFTSPEC, progName, ERR_VOID, 0);
break;
}
i++;
GetWORD(tmpW,i);
}
}
void SetUpFFT()
{
fftRec *tmp;
branchRec *b;
/* the following may be a conflict, since it is now performed "twice." */
tmp = baseOfFFT;
while (tmp != NULL) {
switch (tmp->outType) {
case branchVoltage:
if (tmp->n1 == 0)
tmp->vo1 = GND;
else
tmp->vo1 = &vv_x[ID_nindex[tmp->n1]];
if (tmp->n2 == 0)
tmp->vo2 = GND;
else
tmp->vo2 = &vv_x[ID_nindex[tmp->n2]];
break;
case branchCurrent:
GetBranchByName(&b, tmp->name);
tmp->vo1 = &vv_x[b->ID.UU.n[0]];
tmp->vo2 = GND;
break;
}
/* create some space for the list of numbers ... */
/* must determine whether data space has been created for this entry yet */
tmp->data = get_fft_mrec(tmp);
tmp = tmp->next;
}
}
void RunFFT()
{
fftMRec *tmp;
fftClist *holder;
if (baseOfMFFT == NULL)
return;
tmp = baseOfMFFT;
while (tmp != NULL) {
holder = tmp->current;
/* should place parameters in appropriate bins... */
holder->x[holder->depth].re = (*tmp->vo1 - *tmp->vo2);
holder->x[holder->depth].im = (double) 0.0;
tmp->fft_length++;
holder->depth++;
if (holder->depth == FFT_ARRAY_SIZE) {
tmp->current = (fftClist *) c_malloc(sizeof(fftClist));
holder->next = tmp->current;
tmp->current->depth = 0;
tmp->current->next = NULL;
}
tmp = tmp->next;
}
}
void process_fft()
{
fftMRec *tmp;
gck_complex *x;
long length;
long position;
fftClist *clistptr;
fftClist *c2;
if (baseOfMFFT == NULL)
return;
tmp = baseOfMFFT;
while (tmp != NULL) {
x = (gck_complex *) c_malloc(sizeof(gck_complex) * tmp->fft_length);
/* here rdh stress */
clistptr = tmp->base;
position = 0;
for (length = 0; length < tmp->fft_length; length++) {
x[length].re = clistptr->x[position].re;
x[length].im = clistptr->x[position].im;
position++;
if (position == FFT_ARRAY_SIZE) {
c2 = clistptr->next;
free(clistptr);
clistptr = c2;
position = 0;
}
}
window_data(x,tmp->fft_length);
compute_fft(x,tmp->fft_length);
tmp->result = x;
tmp = tmp->next;
}
}
void print_fft()
{
fftRec *tmp;
long length;
long stages;
long counter;
long index;
register double frequency_step;
double z;
double value;
FILE *fft_file;
if (fft_to_file) {
fft_file = fopen(fft_filename,"w");
if (fft_file == NULL)
errorReport(ERROR_FILE_FFT_CANNOT_OPEN,fft_filename,ERR_VOID,0);
}
else {
fft_file = stdout;
}
if (baseOfFFT == NULL)
return;
length = baseOfFFT->data->fft_length;
/* frequency step fixed April 13, 1992 RdH */
z = ((double) deltaT * (double) NoOfTicks / (double) samplesOn);
frequency_step = (double) (1.0 / (z * (double) length));
stages = get_max_stages(length);
for (counter = 0; counter < (length / 2); counter++) {
fprintf(fft_file,"% .5E%c",(counter * frequency_step),TAB_CHAR);
tmp = baseOfFFT;
if (fft_must_reorder)
index = bit_reversed_index(counter,stages);
else
index = counter;
while (tmp != NULL) {
switch (tmp->fptype) {
case FP_REAL:
fprintf(fft_file,"% .5E",tmp->data->result[index].re);
break;
case FP_IMAG:
fprintf(fft_file,"% .5E",tmp->data->result[index].im);
break;
case FP_MAGN:
value = (double)
(sqrt(pow(tmp->data->result[index].re,2.0) + pow(tmp->data->result[index].im,2.0)));
fprintf(fft_file,"% .5E",value);
break;
case FP_PHASE:
value = atan2(tmp->data->result[index].im,tmp->data->result[index].re);
fprintf(fft_file,"% .5E",value);
break;
case FP_DB:
value = (double)
(20.0 * log10(sqrt(pow(tmp->data->result[index].re,2.0)
+ pow(tmp->data->result[index].im,2.0))));
fprintf(fft_file,"% .5E",value);
break;
}
if (tmp->next != NULL)
fprintf(fft_file,"%c",TAB_CHAR);
tmp = tmp->next;
}
fprintf(fft_file,"\n");
}
}
long bit_reversed_index(index,stages)
long index;
long stages;
{
long result;
long stage;
long hold;
result = 0;
for (stage = 0; stage < stages; stage++) {
hold = (index >> stage) & 0x01;
result += (hold << (stages - stage - 1));
}
return(result);
}
long get_max_stages(N)
long N;
{
long i;
for (i = 0; i < 30; i++) {
if (((N >> i) & 0x01) == 1)
break;
}
return((long) i);
}
void compute_fft(fft_x,N)
gck_complex *fft_x;
long N;
{
long butterfly;
long butterfly_first_position;
long butterfly_index_spacing;
long max_butterflies;
long max_stages;
long stage;
long element_position;
register double ti;
register double tj;
long inner_loop;
long inner_loop_upper_limit;
double alpha;
long outer_loop;
long step_factor;
long twiddle;
long loop;
double twi[2],twj[2];
gck_complex *twiddle_array;
gck_complex *xn,*hold;
boolean can_fft;
double junk1;
long junk2;
register long k,n;
if (PRINT_PAINFUL)
fprintf(stderr,"N is %ld \n",N);
k = 0;
for (n = 0; n < 30; n++) { /* 30 because a long has 32 bits */
if (((N >> n) & 0x01) == (long) 1)
k++;
}
if (k == 1)
can_fft = TRUE;
else
can_fft = FALSE;
if (can_fft) {
if (PRINT_PAINFUL)
fprintf(stderr,"Performing FFT ... \n");
twiddle_array = (gck_complex *) c_malloc(sizeof(gck_complex) * N);
max_stages = get_max_stages(N);
max_butterflies = (N >> 1);
butterfly_index_spacing = (N >> 1); /* initial spacing */
for (stage = 0; stage < max_stages; stage++) {
if (stage != 0) {
element_position = 0;
/* compute twiddle factor... */
step_factor = (1 << (stage - 1) * 1);
for (outer_loop = 0; outer_loop < (1 << (stage * 1)); outer_loop++) {
inner_loop_upper_limit = (N / (step_factor << 1));
for (inner_loop = 0; inner_loop < inner_loop_upper_limit; inner_loop++) {
twiddle = (step_factor * (outer_loop % 2) * inner_loop);
alpha = (-2.0 * PI * (double) twiddle) / ((double) N);
twiddle_array[element_position].re = cos(alpha);
twiddle_array[element_position].im = sin(alpha);
element_position++;
}
}
}
/* twiddle_all_elements() */
if (stage != 0) {
for (element_position = 0; element_position < N; element_position++) {
ti = fft_x[element_position].re;
tj = fft_x[element_position].im;
fft_x[element_position].re = (ti * twiddle_array[element_position].re)
- (tj * twiddle_array[element_position].im);
fft_x[element_position].im = (ti * twiddle_array[element_position].im)
+ (tj * twiddle_array[element_position].re);
}
}
/* must butterfly all next */
for (butterfly = 0; butterfly < max_butterflies; butterfly++) {
butterfly_first_position = butterfly / butterfly_index_spacing;
butterfly_first_position *= (butterfly_index_spacing << 1);
butterfly_first_position += (butterfly % butterfly_index_spacing);
/* radix_2_butterfly(butterfly_first_position,butterfly_index_spacing); */
for (loop = 0; loop < 2; loop++) {
element_position = butterfly_first_position + (loop * butterfly_index_spacing);
twi[loop] = fft_x[element_position].re;
twj[loop] = fft_x[element_position].im;
}
/* Oppenheim and Schafer, Discrete Time Signal Processing, Figure 9.6, page 591 */
element_position = butterfly_first_position;
fft_x[element_position].re = twi[0] + twi[1];
fft_x[element_position].im = twj[0] + twj[1];
element_position += butterfly_index_spacing;
fft_x[element_position].re = twi[0] - twi[1];
fft_x[element_position].im = twj[0] - twj[1];
}
butterfly_index_spacing >>= 1;
}
/* bit reversal performed during printout...*/
free(twiddle_array);
}
else {
if (PRINT_NORMAL) {
fprintf(stderr,"*Warning: Performing DFT ... This may take a while.\n");
}
fft_must_reorder = FALSE;
/* now, a DFT */
/* Equation 8.61, Oppenheim and Schafer, Discrete-Time Signal Processing */
twiddle_array = (gck_complex *) c_malloc(sizeof(gck_complex) * N);
xn = (gck_complex *) c_malloc(sizeof(gck_complex) * N);
for (k = 0; k < N; k++) {
xn[k].re = fft_x[k].re;
xn[k].im = fft_x[k].im;
alpha = (-2.0 * PI * (double) k) / ((double) N);
twiddle_array[k].re = cos(alpha);
twiddle_array[k].im = sin(alpha);
}
for (k = 0; k < N; k++) {
ti = 0.0; tj = 0.0;
for (n = 0; n < N; n++) {
element_position = (k * n) % N;
ti += xn[n].re * twiddle_array[element_position].re - xn[n].im * twiddle_array[element_position].im;
tj += xn[n].re * twiddle_array[element_position].im + xn[n].im * twiddle_array[element_position].re;
}
fft_x[k].re = ti;
fft_x[k].im = tj;
}
free(xn);
free(twiddle_array);
}
}
void window_data(fft_x,N)
gck_complex *fft_x;
long N;
{
long n;
long spare;
double scalar;
if ((PRINT_PAINFUL) && (fft_window != WINDOW_RECTANGULAR))
fprintf(stderr,"Windowing Data ... \n");
switch(fft_window) {
case WINDOW_RECTANGULAR:
/* Whew, that was easy. */
break;
case WINDOW_TRIANGULAR:
if (((double) N / 2.0) - (double) (long) (N / 2) == 0.5)
spare = 1;
else
spare = 0;
for (n = 0; n <= (N/2); n++) {
scalar = 2.0 * (double) n / N;
fft_x[n].re *= scalar;
}
for (n = (N/2) + spare; n < N; n++) {
scalar = 2.0 - 2.0 * (double) n / N;
fft_x[n].re *= scalar;
}
break;
case WINDOW_HANN:
for (n = 0; n < N; n++) {
scalar = 0.50 - 0.50 * cos(TWOPI * (double) n / (double) N);
fft_x[n].re *= scalar;
}
break;
case WINDOW_HAMMING:
for (n = 0; n < N; n++) {
scalar = 0.54 - 0.46 * cos(TWOPI * (double) n / (double) N);
fft_x[n].re *= scalar;
}
break;
case WINDOW_BLACKMAN:
for (n = 0; n < N; n++) {
scalar = 0.42 - 0.50 * cos(TWOPI * (double) n / (double) N) + 0.08 * cos(2.0 * TWOPI * (double) n / (double) N);
fft_x[n].re *= scalar;
}
break;
}
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.