-
December 21st, 2020, 11:01 AM
#1
Can I ask a question about FFT here?
Hello,
I am not sure I am allowed to ask FFT question here. I am using Visual Studio to program it.
First I've got a FFT function handy. But I will need to collect a certain amount of samples before passing them into it.
When the FFT is done processing that lot of samples, there are a vector of complex numbers, I understand they are from lowest to highest frequencies,
but they are not labelled. When I display values, I clear that lot of samples and that's the end of the window
The FFT function basically comes from the book called Toolchest, a very commonly placed book
Code:
/*--------------------------------------------------------------------*/
#include "FFT42.h"
double pi_ = 3.1416;
void matherr_(char*, int) { }
void fft42(Complex_Vector data, unsigned n, int isign)
{
/*
subroutine performs an n-point, radix "4+2" Cooley-Tukey fft on complex array
"data". This routine is about twice as fast as the standard radix-2 approach.
n must be a power of two
data[] is the complex input array, and, on exit, it is also the complex
output array (i.e. the input is destroyed)
forward transform ==> isgn=-1
inverse transform ==> isgn=+1
*/
unsigned nm1,ndv2,j,k,l,m,ipar,mmax,lmax,kmin,kdif,kstep,k1,k2,k3,k4;
double tempr,tempi,theta,wr,wi,wstpr,wstpi,w2r,w2i,w3r,w3i,fn;
double u1r,u1i,u2r,u2i,u3r,u3i,u4r,u4i,t2r,t2i,t3r,t3i,t4r,t4i;
unsigned lfft;
char *fname = "fft42";
if(data == NULL) {
matherr_("fname", E_NULLPTR);
return;
}
/* see if fft length, n, is a power of two */
m = (unsigned) ceil(log((double) n)/0.693147181);
lfft = 2;
for(j=1; j<m; j++) lfft += lfft;
if(lfft != n) {
matherr_(fname, E_FFTPOWER2);
return;
}
/* m = log2 (n) */
nm1=n-1;
wr = wi = wstpr = wstpi = 0.0;
ndv2 = n>>1;
j=0;
for (k=0; k < nm1; k++) {
if (k < j) { /* bit reversal stage */
tempr=data[j].r;
tempi=data[j].i;
data[j].r=data[k].r;
data[j].i=data[k].i;
data[k].r=(Real)tempr;
data[k].i=(Real)tempi;
}
m=ndv2;
while (m < j+1) { j -= m; m >>= 1; }
j += m;
}
ipar=n;
while (ipar > 2) ipar >>= 2;
if (ipar == 2) {
for (k1=0; k1<=nm1; k1+=2) {
k2=k1+1;
tempr=data[k2].r;
tempi=data[k2].i;
data[k2].r=data[k1].r-tempr;
data[k2].i=data[k1].i-tempi;
data[k1].r=data[k1].r+tempr;
data[k1].i=data[k1].i+tempi;
}
}
mmax=2;
while (mmax < n) {
lmax=MAX(4,(mmax >> 1));
if(mmax > 2) {
theta=-pi_/mmax;
if(isign >= 0) theta=-theta;
wr=cos(theta);
wi=sin(theta);
wstpr=-2.0*wi*wi;
wstpi=2.0*wr*wi;
}
for (l=2; l<=lmax; l+=4) {
m=l;
if(mmax > 2) {
zero:
w2r=wr*wr-wi*wi;
w2i=2.0*wr*wi;
w3r=w2r*wr-w2i*wi;
w3i=w2r*wi+w2i*wr;
}
kmin=ipar*m; kmin >>= 1;
if(mmax <= 2) kmin = 0;
kdif=ipar*mmax; kdif >>= 1;
five:
kstep = (kdif << 2);
for (k1=kmin; k1<=nm1; k1+=kstep) {
k2=k1+kdif;
k3=k2+kdif;
k4=k3+kdif;
if(mmax <= 2) {
u1r=data[k1].r+data[k2].r;
u1i=data[k1].i+data[k2].i;
u2r=data[k3].r+data[k4].r;
u2i=data[k3].i+data[k4].i;
u3r=data[k1].r-data[k2].r;
u3i=data[k1].i-data[k2].i;
if(isign < 0) {
u4r=data[k3].i-data[k4].i;
u4i=data[k4].r-data[k3].r;
goto ten;
}
u4r=data[k4].i-data[k3].i;
u4i=data[k3].r-data[k4].r;
goto ten;
}
t2r=w2r*data[k2].r-w2i*data[k2].i;
t2i=w2r*data[k2].i+w2i*data[k2].r;
t3r=wr*data[k3].r-wi*data[k3].i;
t3i=wr*data[k3].i+wi*data[k3].r;
t4r=w3r*data[k4].r-w3i*data[k4].i;
t4i=w3r*data[k4].i+w3i*data[k4].r;
u1r=data[k1].r+t2r;
u1i=data[k1].i+t2i;
u2r=t3r+t4r;
u2i=t3i+t4i;
u3r=data[k1].r-t2r;
u3i=data[k1].i-t2i;
if(isign < 0) {
u4r=t3i-t4i;
u4i=t4r-t3r;
goto ten;
}
u4r=t4i-t3i;
u4i=t3r-t4r;
ten:
data[k1].r=(Real)u1r+u2r;
data[k1].i=(Real)u1i+u2i;
data[k2].r=(Real)u3r+u4r;
data[k2].i=(Real)u3i+u4i;
data[k3].r=(Real)u1r-u2r;
data[k3].i=(Real)u1i-u2i;
data[k4].r=(Real)u3r-u4r;
data[k4].i=(Real)u3i-u4i;
}
kmin = (kmin << 2);
kdif = kstep;
if(kdif-n) goto five;
m=mmax-m;
if(isign < 0) {
tempr=wr;
wr=-wi;
wi=-tempr;
goto fifteen;
}
tempr=wr;
wr=wi;
wi=tempr;
fifteen:
if(m > lmax) goto zero;
tempr=wr;
wr=wr*wstpr-wi*wstpi+wr;
wi=wi*wstpr+tempr*wstpi+wi;
}
ipar=3-ipar;
mmax += mmax;
}
if(isign < 0) return;
fn = 1.0/n;
for (k=0; k < n; k++) { data[k].r *= fn; data[k].i *= fn;}
}
So the questions are:
1st) How to find out the names of the labels? (the frequencies), hence also finding out the intervals between the frequencies.
2nd) Is this windowing technique correct?
Thanks for helping
Jacky
Last edited by luckiejacky; December 21st, 2020 at 11:04 AM.
-
December 21st, 2020, 11:52 AM
#2
Re: Can I ask a question about FFT here?
A good example of the "bad fortran" code written in C.
What do you mean by "label" in your question?
Victor Nijegorodov
-
December 21st, 2020, 12:00 PM
#3
Re: Can I ask a question about FFT here?
The program wasn't written by me, however, when I said label, I meant the frequency "number", for example, 10KHz, 20KHz.
When the data comes out of the function, I receive a stack of vectors. But I do not know what does it represent.
-
December 21st, 2020, 04:58 PM
#4
Re: Can I ask a question about FFT here?
Don't worry! I didn't mean its your own code!
Victor Nijegorodov
Posting Permissions
- You may not post new threads
- You may not post replies
- You may not post attachments
- You may not edit your posts
-
Forum Rules
|
Click Here to Expand Forum to Full Width
|