On a sunny day (Wed, 30 Jul 2014 20:40:49 -0400) it happened "Maynard A. Philbrook Jr." wrote in :
Octave is basicaly an open source version (clone?) of mathlab. That means mathlab files run in octave (at least the ones I tried).
formatting link
Some were not so lucky:
formatting link
Anyways, fft is very basic, just a few lines of C code really, I kept this from 1988 Wireless World:
#include #include #include
/* Fourier analysis program Ref. used: Electronics and Wireless World october 1988, article by Howard J.Hutchings "Interfacing and signal processing in C"
float rad(float angle) { return (angle * M_PI) / 180.0; } /* end function rad */
int main(int argc, char **argv) { FILE *fptr; float frequency;
//fptr = fopen("six_periods.dat", "r"); //fptr = fopen("66hz_16ms_1_period.dat", "r"); fptr = fopen("6600Hz_600us_4_periods.dat", "r"); /* 4.6875 us / step, 128 samples */ if(! fptr) { fprintf(stderr, "could not open file four_periods.dat for read, aborting.\n");
exit(1); }
int ival;
i = 0; while(1) { a = fscanf(fptr, "%x", &ival); if(a == EOF) break; if(a == 0) break;
z[i] = (float) ival;
fprintf(stderr, "i=%d %.2f\n", i, z[i]); i++; }
/* display shows 4 periods of 6600 Hz.
1 period of 6600 Hz = 1 / 6600 = 1.52. 10^-4 = 152 us
4 periods makes 4 x 152 = 600 uS
128 samples per screen, so per sample the time is 600 / 128 = 4.73 us.
/* Hamming window */ for(i = 0;i < samples;i++) { angle1 = (360.0 / (float)samples) * (float)i; window[i] = 0.5 * (1 - cos(rad(angle1))); /* Hamming window, 0 at start and end of samples, 1 at centre, sine function */ z[i] = z[i] * window[i]; /* Hamming */ }
for(i = 0; i < samples; i++) { fprintf(stderr, "i=%d window[i]=%.2f z[i]=%.2f\n", i, window[i], z[i]); }
fprintf(stderr, "WAS 1\n");
/* discrete Fourier transform Note: for time interval t, samples n, the spectral resolution df = 1 / (n.t) by multipliction with sine function or cosine function for real and imag part the present of a frequency component can be determined. (where the sine function is at maximum that frequency is expected, 90 deg. shifted for the cosine function imag part.)
*/
ic = 0; /* calculated component cnt */ for(k = 0; k < (samples / 2); k += 1.0) /* samples / 2 = niquist freq., above that mirror */ /* for each frequency */ { realsum = 0.0;/*reset*/ imagsum = 0.0; modulus = 0.0; for(i = 0; i < samples; i++) /* test how much we wiggle for this frequency */ { angle = (360.0 / (float)samples) * (float)i * (float)k; realsum += z[i] * cos(rad(angle)) / (float)samples; imagsum += z[i] * sin(rad(angle)) / (float)samples;
//fprintf(stderr, "===========modulus=%.6f\n", modulus); modulus = sqrt(modulus); /* back to normal values */
result[ic] = modulus; /* save result */
fprintf(stderr, "ic=%d ===========modulus_sqrt=%.6f result[ic]=%.6f\n", ic, modulus, result[ic]); ic++; /* the result is at maximum .5 * maximum input value, due to Hamming window */ }
for(i = 0; i < ic; i++) { frequency = (float)i * deltaf;
fprintf(stdout, "%d %.2f %.12f\n", i, frequency, result[i]); }
Anyways, people usually use libraries, and then sometimes start using them the wrong way. Open source fftw library is available and I have used that,
formatting link
formatting link
Still needs some work, will work on it once I need it... but it (libfftw) is huge, and many functions you will likely not need.
The above C code is easily rewritten in some other language, I have rewritten it in PIC 18 asm, used signed 32 bit integers, and use it here (used sine lookup table):
formatting link
just a few kB for the whole thing... Nice ASCII output it can do:
formatting link
I think by the look of it python interfaces with libfftw. As to python, I am not sure there is any advantage over using C, on the contrary. Ever so often somebody who is too lazy to learn C invents a new language, one of the worst examples is C++, a crime against humanity,
Except that the code you have posted is *NOT* an FFT! Your DFT code is O(N^2) an FFT is O(NlogN)
And if they roll their own they frequently don't know what they are doing and get into even more trouble. Most numerical library code is extremely well tested and characterised for numerical robustness.
If it does then that is a very good and fast portable FFT library.
Actually no. That WW code isn't an FFT at all! That is a long handed DFT done *THE* most inefficient way imaginable.
This is about the shortest C snippet for a canonical radix two FFT. It uses the original notation and is translated from a FORTRAN original. The lines commented out // would be the naive way to step the angle and considerably less accurate. The next step up in accuracy is to make up a precomputed lookup table containing the discrete 2^N roots of -1. Computing trig functions is still a little slow even on modern CPUs which is why FFTW and the like cache their computed twiddle factors.
Basically it splits down a 2^N FFT as 2x 2^N-1 FFTs combined together with suitable twiddle factors and so on down to the trivial 2 case.
/* FFT is a simple Fourier Transform based on the original Cooley-Tukey factorisation and notation, it takes a real array of length [0..2^N-1] Array length must be of the form 2^N, even index for reals and odd indexes are imaginary. The result is normalised after the -1 transform
*/
void FFT(float *x, unsigned int n, int ind) { unsigned int i, j, le, le2, ip, n2; double u1, u2, u3, t1, t2; double arg, s, s2;
On a sunny day (Thu, 31 Jul 2014 09:31:00 +0100) it happened Martin Brown wrote in :
Right!
I have this one:
/* fix_fft.c - Fixed-point in-place Fast Fourier Transform */ /* All data are fixed-point short integers, in which -32768 to +32768 represent -1.0 to +1.0 respectively. Integer arithmetic is used for speed, instead of the more natural floating-point.
For the forward FFT (time -> freq), fixed scaling is performed to prevent arithmetic overflow, and to map a 0dB sine/cosine wave (i.e. amplitude = 32767) to two -6dB freq coefficients. The return value is always 0.
For the inverse FFT (freq -> time), fixed scaling cannot be done, as two 0dB coefficients would sum to a peak amplitude of 64K, overflowing the 32k range of the fixed-point integers. Thus, the fix_fft() routine performs variable scaling, and returns a value which is the number of bits LEFT by which the output must be shifted to get the actual amplitude (i.e. if fix_fft() returns 3, each value of fr[] and fi[] must be multiplied by 8 (2**3) for proper scaling. Clearly, this cannot be done within fixed-point short integers. In practice, if the result is to be used as a filter, the scale_shift can usually be ignored, as the result will be approximately correctly normalized as is.
Written by: Tom Roberts 11/8/89 Made portable: Malcolm Slaney 12/15/94 snipped-for-privacy@interval.com Enhanced: Dimitrios P. Bouras 14 Jun 2006 snipped-for-privacy@ieee.org
*/
#define N_WAVE 1024 /* full length of Sinewave[] */ #define LOG2_N_WAVE 10 /* log2(N_WAVE) */
/* Henceforth "short" implies 16-bit word. If this is not the case in your architecture, please replace "short" with a type definition which *is* a 16-bit word.
*/
/* Since we only use 3/4 of N_WAVE, we define only this many samples, in order to conserve data space.
/* FIX_MPY() - fixed-point multiplication & scaling. Substitute inline assembly for hardware-specific optimization suited to a particluar DSP processor. Scaling ensures that result remains 16-bit.
*/ inline short FIX_MPY(short a, short b) { /* shift right one less bit (i.e. 15-1) */ int c = ((int)a * (int)b) >> 14; /* last bit shifted out = rounding-bit */ b = c & 0x01; /* last shift + rounding bit */ a = (c >> 1) + b; return a; }
/* fix_fft() - perform forward/inverse fast Fourier transform. fr[n],fi[n] are real and imaginary arrays, both INPUT AND RESULT (in-place FFT), with 0 = 1; } while (mr+l > nn); mr = (mr & (l-1)) + l;
if (mr 16383) { shift = 1; break; } } if (shift) ++scale; } else { /* fixed scaling, for proper normalization -- there will be log2(n) passes, so this results in an overall factor of 1/n, distributed to maximize arithmetic accuracy. */ shift = 1; } /* it may not be obvious, but the shift will be performed on each data point exactly once, during this pass. */ istep = l = 1; } for (i=m; i>= 1; qi >>= 1; } fr[j] = qr - tr; fi[j] = qi - ti; fr[i] = qr + tr; fi[i] = qi + ti; } } --k; l = istep; } return scale; }
/* fix_fftr() - forward/inverse FFT on array of real numbers. Real FFT/iFFT using half-size complex FFT by distributing even/odd samples into real/imaginary arrays respectively. In order to save data space (i.e. to avoid two arrays, one for real, one for imaginary samples), we proceed in the following two steps: a) samples are rearranged in the real array so that all even samples are in places 0-(N/2-1) and all imaginary samples in places (N/2)-(N-1), and b) fix_fft is called with fr and fi pointing to index 0 and index N/2 respectively in the original array. The above guarantees that fix_fft "sees" consecutive real samples as alternating real and imaginary samples in the complex array.
*/ int fix_fftr(short f[], int m, int inverse) { int i, N = 1
On a sunny day (Thu, 31 Jul 2014 09:31:00 +0100) it happened Martin Brown wrote in :
PS actually for a small number of samples that code is just fast enough, I used if for scope_pid, because tehre are more than just speed constraints.
1 the display is only 128 x 64, so 128 samples per bin or screen.
2 there is only a few kb code space.
3 there is only 256 bytes RAM plus some vars.
4 you have all the time in the world to do a transform, as soon as the samples are in. is a storage scope after all...
Since the real DFT is simpler to code, needs less space, it just works fine here on a 64 MHz PIC 18. I just use that memory several times and sometimes the bits twice. Although that last thing seems impossible, you see I display the FFT result in the wave form result. Why not:
formatting link
So 'speed' as such is not always an issue.
As to resolution (in amplitude) 64 pixels vertical I had to bit shift right several times to fit the result... We always need to be practical.
Ah Labview---finally we can have spaghetti code that looks like spaghetti.
Did you write a prototype for them, which they chose to ignore and reimplement---because Windows, maybe? It used to be hard to port Linux code but now with Qt and other multiplatform options like Python or JS/ HTML5, it's just reflexive lack of imagination.
The proto did need prettying up quite a bit--it was made on an optical breadboard with mostly stuff I had in my drawer, plus some hand-wired electronics, an RC airplane servo to turn the grating, a LabJack to to the A/D and control functions. It used JB Weld putty for a number of the mechanical bits, e.g. the holders for the bulb and the fibre bundles. (Great stuff, JB Weld.)
The s/w was a multithreaded console-mode program running on Windows. The commonality with Linux is pretty good, but some things would have had to be changed, e.g. the code to start a subprocess running Gnuplot and talk to it over a pipe. I already have a library that implements some little things you really can't get along without in a console mode control program, e.g. _kbhit(), which I implemented using select(), as well as portable mutexes and event semaphores.
My gizmo was ugly, took six weeks to build, cost a total of about $80k including a detailed photon budget, hardware, and software, and basically worked very well. There was some stray light, and the RC servo wasn't really repeatable enough, but those things would have been fairly easy to fix. Of course, you can't very well show the FDA something made from toy parts. ;)
Theirs took almost a year, cost many hundreds of thousands of dollars, looks much prettier, is somewhat better optically, but works much less well and is much harder to fix because every time I peel off one onion layer I find more buried treasure. The most recent thing is a "PID" control loop on the chopper that exhibits gross amounts of hunting and doesn't null out a steady-state disturbance, despite my varying the PID parameters over a wide range.
The firmware appears to run the PID control not from a timer interrupt, but in the main housekeeping loop. That's very poor practice but probably only a minor contribution in this case, since the irregular hunting is so bad (3-4 degrees or more at 1 kHz). Even a jittery loop should null out the phase error due to a constant drag force. (In this case the drag was from a steady stream of compressed air flowing tangentially to the rotor.) I measured this with nothing more than an HP 3325A synthesizer driving a white LED, making a strobe light tuned to match the chopper frequency (I did put a negative offset on the synth output and cranked up the amplitude, which narrows the conduction angle on the sinusoidal peaks.
I haven't ploughed all the way through the firmware source code yet, but I suspect that they're sampling the optointerrupter every time round the main loop, and occasionally miss the edge. The fact that the hunting has an amplitude of just about one chopper tooth is very suspicious. (The proto used a 4046 PLL, and was very stable. Twenty cents goes a long way sometimes.)
Plus I'm told that they broke a $100k FTIR lent to them by the client.
I'm _not_ looking forward to finding out what their digital lock-in code is doing with this signal. (The proto used an analogue lock-in that cost about $2 all up, and worked great.)
As I say, those folks are not a candidate for repeat business.
Cheers
Phil Hobbs
--
Dr Philip C D Hobbs
Principal Consultant
ElectroOptical Innovations LLC
Optics, Electro-optics, Photonics, Analog Electronics
160 North State Road #203
Briarcliff Manor NY 10510
hobbs at electrooptical dot net
http://electrooptical.net
I have two versions, one using getchar() rather like Jan's, and one using select(). They both use tcsetaddr()/tcgetaddr() instead of system() calls to get in and out of raw mode, but they basically do the same thing.
Cheers
Phil Hobbs
--
Dr Philip C D Hobbs
Principal Consultant
ElectroOptical Innovations LLC
Optics, Electro-optics, Photonics, Analog Electronics
160 North State Road #203
Briarcliff Manor NY 10510
hobbs at electrooptical dot net
http://electrooptical.net
ElectronDepot website is not affiliated with any of the manufacturers or service providers discussed here.
All logos and trade names are the property of their respective owners.