I am not an expert in DSP, in fact I know some of the theory but I have never really written a DSP program. I must write one in order to detect the presence of a carrier signal (either 2222Hz or 2475Hz) The signal is present for 16 to 20 ms. I am writing a program for the Atmel UC3A3 (AVR32). So far, I have programmed the ADC and the DMA to write conversion results in a double buffer. The DMA writes the result in memory as soon as the conversion is finished. When N samples have been converted, an interrupt is raised and the buffers are switched. I intend to use the idle buffer, where the samples have been stored, in order to perform the detection. I imagine that I could do a FFT of the buffer and look for the results around the frequencies of interest. Is that valid or are there problems with that approach? What should the acquisition frequency be? (5kHz that is 2*2475Hz ?) Should I first filter the signal? How can I determine the presence of the carrier signal when the computed frequencies don't match exactly the looked for frequencies? You see that those are basic questions, I hope you will be patient.
carrier detection
Started by ●October 15, 2014
Reply by ●October 15, 20142014-10-15
On Wed, 15 Oct 2014 17:51:01 +0200, Lanarcam wrote:> I am not an expert in DSP, in fact I know some of the theory but I have > never really written a DSP program. > > I must write one in order to detect the presence of a carrier signal > (either 2222Hz or 2475Hz) > > The signal is present for 16 to 20 ms. > > I am writing a program for the Atmel UC3A3 (AVR32). > > So far, I have programmed the ADC and the DMA to write conversion > results in a double buffer. The DMA writes the result in memory as soon > as the conversion is finished. > > When N samples have been converted, an interrupt is raised and the > buffers are switched. > > I intend to use the idle buffer, where the samples have been stored, in > order to perform the detection.This will fall down if you get part of your carrier in one buffer and part in another. I would do this with some sort of a filter (more on that later) that acts on all the received samples with no seams at the buffer boundaries.> I imagine that I could do a FFT of the buffer and look for the results > around the frequencies of interest. Is that valid or are there problems > with that approach?It is complex, error ridden (because you're looking in too large of chunks), hard to interpret the results, and very processor intensive for the task you're trying to accomplish -- but other than that it's fine.> What should the acquisition frequency be? (5kHz that is 2*2475Hz ?)AAAGH! Noooooo! Please see <http://wescottdesign.com/articles/Sampling/ sampling.pdf> -- it should explain why that's a bad idea. IF you have enough horsepower (which you should) and IF you don't have lots of noise to contend with, then sampling at 10kHz should reduce your level of difficulty dramatically. Sampling at 20 or even 40kHz is not to be sneezed at.> Should I first filter the signal?Before you FFT, or before you sample? Filtering before sampling depends on what you're trying to accomplish, and the source of the signal -- basically, if you don't have any out-of-band noise to alias into your intended signal, you don't need to filter.> How can I determine the presence of the carrier signal when the computed > frequencies don't match exactly the looked for frequencies? > > You see that those are basic questions, I hope you will be patient.The easiest and most robust (and fairly computationally intensive) way to do this is to use a pair of FIR bandpass filters which run on the data. Monitor their outputs, and if one of them goes above a threshold (which you must determine, either theoretically or experimentally) declare that a tone is present. Second easiest and robust, and far less computationally intensive, is to use a pair of IIR bandpass filters which run on the data. Everything else is the same as above. There's various other things that you can try, but neither of the above two should be using a significant portion of your overall processing power: if you need to quibble over fractions of a percent of the total available processor cycles, then you've got way bigger problems than your stated one. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com
Reply by ●October 15, 20142014-10-15
Hi,>>I must write one in order to detect the presence of a carrier signal(either 2222Hz or 2475Hz) Detect over a common multiple of cycles, which is n*202 and n*225 cycles, or n*90.909090_ ms. (To get this number, divide both frequencies by the common factor 11) - multiply input signal with both a sine wave and a cosine wave of the detection frequency - add the sine products into one accumulator, the cosine products into another - after the detection interval, square both accumulators and add - the result is the power detected at the frequency (take square root for magnitude) - Compare threshold for any compromise between sensitivity and number of "false alarms". - do this in parallel for both frequencies. A computationally efficient solution is the so-called "Goertzel" algorithm. It may be more complex than it looks... See here: http://www.mathworks.se/products/demos/signaltlbx/dtmf/dtmfdemo.html http://en.wikipedia.org/wiki/Goertzel_algorithm _____________________________ Posted through www.DSPRelated.com
Reply by ●October 15, 20142014-10-15
Le 15/10/2014 18:44, Tim Wescott a écrit :> On Wed, 15 Oct 2014 17:51:01 +0200, Lanarcam wrote: > >> I am not an expert in DSP, in fact I know some of the theory but I have >> never really written a DSP program. >> >> I must write one in order to detect the presence of a carrier signal >> (either 2222Hz or 2475Hz) >> >> The signal is present for 16 to 20 ms. >> >> I am writing a program for the Atmel UC3A3 (AVR32). >> >> So far, I have programmed the ADC and the DMA to write conversion >> results in a double buffer. The DMA writes the result in memory as soon >> as the conversion is finished. >> >> When N samples have been converted, an interrupt is raised and the >> buffers are switched. >> >> I intend to use the idle buffer, where the samples have been stored, in >> order to perform the detection. > > This will fall down if you get part of your carrier in one buffer and part > in another. I would do this with some sort of a filter (more on that > later) that acts on all the received samples with no seams at the buffer > boundaries. > >> I imagine that I could do a FFT of the buffer and look for the results >> around the frequencies of interest. Is that valid or are there problems >> with that approach? > > It is complex, error ridden (because you're looking in too large of > chunks), hard to interpret the results, and very processor intensive for > the task you're trying to accomplish -- but other than that it's fine. > >> What should the acquisition frequency be? (5kHz that is 2*2475Hz ?) > > AAAGH! Noooooo! Please see <http://wescottdesign.com/articles/Sampling/ > sampling.pdf> -- it should explain why that's a bad idea. > > IF you have enough horsepower (which you should) and IF you don't have > lots of noise to contend with, then sampling at 10kHz should reduce your > level of difficulty dramatically. Sampling at 20 or even 40kHz is not to > be sneezed at. > >> Should I first filter the signal? > > Before you FFT, or before you sample? Filtering before sampling depends > on what you're trying to accomplish, and the source of the signal -- > basically, if you don't have any out-of-band noise to alias into your > intended signal, you don't need to filter. > >> How can I determine the presence of the carrier signal when the computed >> frequencies don't match exactly the looked for frequencies? >> >> You see that those are basic questions, I hope you will be patient. > > The easiest and most robust (and fairly computationally intensive) way to > do this is to use a pair of FIR bandpass filters which run on the data. > Monitor their outputs, and if one of them goes above a threshold (which > you must determine, either theoretically or experimentally) declare that a > tone is present. > > Second easiest and robust, and far less computationally intensive, is to > use a pair of IIR bandpass filters which run on the data. Everything else > is the same as above. > > There's various other things that you can try, but neither of the above > two should be using a significant portion of your overall processing > power: if you need to quibble over fractions of a percent of the total > available processor cycles, then you've got way bigger problems than your > stated one. >Thank you for your answer, I must "absorb" it, I am no expert. Let's suppose I sample at 10 kHz, it will be difficult to perform the computation (two IIR filters) every 100 microseconds. I don't see how one can avoid buffering the inputs and computing later on. Is there a trick?
Reply by ●October 15, 20142014-10-15
Le 15/10/2014 19:12, mnentwig a �crit :> Hi, > >>> I must write one in order to detect the presence of a carrier signal > (either 2222Hz or 2475Hz) > > Detect over a common multiple of cycles, which is n*202 and n*225 cycles, > or n*90.909090_ ms. > (To get this number, divide both frequencies by the common factor 11) > > - multiply input signal with both a sine wave and a cosine wave of the > detection frequency > - add the sine products into one accumulator, the cosine products into > another > - after the detection interval, square both accumulators and add > - the result is the power detected at the frequency (take square root for > magnitude) > - Compare threshold for any compromise between sensitivity and number of > "false alarms". > - do this in parallel for both frequencies. > > A computationally efficient solution is the so-called "Goertzel" algorithm. > It may be more complex than it looks... > > See here: > http://www.mathworks.se/products/demos/signaltlbx/dtmf/dtmfdemo.html > http://en.wikipedia.org/wiki/Goertzel_algorithm >Thank you very much. I don't quite get: "Detect over a common multiple of cycles, which is n*202 and n*225 cycles, or n*90.909090_ ms." I can see that 1s/11 = 90.909090_ ms and that 11 is a common factor of both frequencies. Does that mean to find n such that: n/11 = k/2222 = m/2475 where n, k and m are integers?
Reply by ●October 15, 20142014-10-15
Hi, what I meant is that the number of cycles in 90.9 ms is exactly 202 (2222 Hz) and 225 (2475 Hz). Other candidates are n*90.9 ms (with n integer), and the number of cycles in that length scales accordingly with n. There is no shorter solution. The integer number of cycles is preferable, because the waveforms are orthogonal over an integer number of cycles. That is, a 2222 Hz input signal results in exactly zero output from the 2475 Hz detection, and vice versa. Most likely, it will work well enough with other lengths for "yes/no" decisions, but I'd start with 90.9 ms for simplicity. _____________________________ Posted through www.DSPRelated.com
Reply by ●October 15, 20142014-10-15
On Wed, 15 Oct 2014 21:41:23 +0200, Lanarcam wrote:> Le 15/10/2014 18:44, Tim Wescott a écrit : >> On Wed, 15 Oct 2014 17:51:01 +0200, Lanarcam wrote: >> >>> I am not an expert in DSP, in fact I know some of the theory but I >>> have never really written a DSP program. >>> >>> I must write one in order to detect the presence of a carrier signal >>> (either 2222Hz or 2475Hz) >>> >>> The signal is present for 16 to 20 ms. >>> >>> I am writing a program for the Atmel UC3A3 (AVR32). >>> >>> So far, I have programmed the ADC and the DMA to write conversion >>> results in a double buffer. The DMA writes the result in memory as >>> soon as the conversion is finished. >>> >>> When N samples have been converted, an interrupt is raised and the >>> buffers are switched. >>> >>> I intend to use the idle buffer, where the samples have been stored, >>> in order to perform the detection. >> >> This will fall down if you get part of your carrier in one buffer and >> part in another. I would do this with some sort of a filter (more on >> that later) that acts on all the received samples with no seams at the >> buffer boundaries. >> >>> I imagine that I could do a FFT of the buffer and look for the results >>> around the frequencies of interest. Is that valid or are there >>> problems with that approach? >> >> It is complex, error ridden (because you're looking in too large of >> chunks), hard to interpret the results, and very processor intensive >> for the task you're trying to accomplish -- but other than that it's >> fine. >> >>> What should the acquisition frequency be? (5kHz that is 2*2475Hz ?) >> >> AAAGH! Noooooo! Please see >> <http://wescottdesign.com/articles/Sampling/ sampling.pdf> -- it should >> explain why that's a bad idea. >> >> IF you have enough horsepower (which you should) and IF you don't have >> lots of noise to contend with, then sampling at 10kHz should reduce >> your level of difficulty dramatically. Sampling at 20 or even 40kHz is >> not to be sneezed at. >> >>> Should I first filter the signal? >> >> Before you FFT, or before you sample? Filtering before sampling >> depends on what you're trying to accomplish, and the source of the >> signal -- basically, if you don't have any out-of-band noise to alias >> into your intended signal, you don't need to filter. >> >>> How can I determine the presence of the carrier signal when the >>> computed frequencies don't match exactly the looked for frequencies? >>> >>> You see that those are basic questions, I hope you will be patient. >> >> The easiest and most robust (and fairly computationally intensive) way >> to do this is to use a pair of FIR bandpass filters which run on the >> data. Monitor their outputs, and if one of them goes above a threshold >> (which you must determine, either theoretically or experimentally) >> declare that a tone is present. >> >> Second easiest and robust, and far less computationally intensive, is >> to use a pair of IIR bandpass filters which run on the data. >> Everything else is the same as above. >> >> There's various other things that you can try, but neither of the above >> two should be using a significant portion of your overall processing >> power: if you need to quibble over fractions of a percent of the total >> available processor cycles, then you've got way bigger problems than >> your stated one. >> > Thank you for your answer, I must "absorb" it, I am no expert. > > Let's suppose I sample at 10 kHz, it will be difficult to perform the > computation (two IIR filters) every 100 microseconds. I don't see how > one can avoid buffering the inputs and computing later on. Is there a > trick?What clock rate are you running your 32-bit processor at? In my experience with ARM processors, if you're clocking above 25Mhz, this shouldn't present a huge load unless you want to do everything with double- precision floating point. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com
Reply by ●October 16, 20142014-10-16
On Wed, 15 Oct 2014 17:51:01 +0200, Lanarcam wrote: (snip) First, I agree with the other comments that a 5K sample rate is too low. 10K might be more reasonable. If the actual carrier frequency is slightly different than what you expect, then the first thing may want to do is estimate it. You could use Gaussian ratio estimation to do so, as described here: https://www.edn.com/design/analog/4352503/EDN--03-03-94-Ratio-detection-precisely-characterizes-signals-amplitude-and-frequenc http://www.ptodirect.com/Results/Patents?query=PN/5214708 (note: the patent was issued in 1993, so it's more than 20 years old, and besides, it's for speech) Basically, you window your input data with a Gaussian window, then do an FFT, after which you take the log of the FFT amplitudes. As it turns out, when you use a window, you're creating a series of band-pass filters in the frequency domain, where each filter is centered at an FFT bin. Now if you have a sinusoid in the input that lies between two consecutive bins of an FFT, then it will contribute to (at least) 2 of the FFT bins - mostly to the one that's higher or lower in frequency. Well, it turns out that you can write 2 equations with 2 unknowns (amplitude and frequency). By using the log values from the FFT, and taking the ratio of the 2 equations, the amplitude drops out, and you can solve for frequency. There is another, slightly different Gaussian method that's used at CERN, but it's a little more complicated, and you probably don't need that. Incidentally, CERN has done a number of detailed accuracy studies on their method - after all, they use it for estimating beam frequencies in a $10 billion supercollider. Now, given your numbers, and presuming a sample rate (s_r) of 10K and a 32 point FFT (N = 32), the 2 bin numbers bracketing your frequencies are: k1 = s_r/N = f1 = 2187.5 --> k1 = 7 k2 = s_r/N = f2 = 2500.0 --> k2 = 8 So you'd collect 32 samples, take the FFT, compute the log of bins 7 and 8, then put those numbers into an equation which will give you the frequency estimate. Then you can use a second equation to compute amplitude. I originally wrote this post as a step by step sequence, but, since you're a beginner, you might need something a bit more detailed. So I cobbled together a C++ program to demonstrate (I used Ariel 10 format, so I hope it prints OK): #include <cstdio> #include <cstdlib> #include <iostream> #include <cmath> using namespace std; void fft_recur(long, double* r, double* i); // prototype declaration int main (int nNumberofArgs, char* pszArgs[ ] ) { const long N = 32; double sample_rate = 10000., r[N] = {0.}, i[N] = {0.}, w[N] = {0.} ; long n; double t, twopi = 6.2831853071795865; double F, fbin, amp, power, c = 3.0; for (n = 0; n < N; n++) { // generate test data t = twopi*n/sample_rate; r[n] = 1.*cos(2222.0*t + twopi/5.) ; // amplitude = 1. frequency = 2222., arbitrary phase } // end for // Gaussian window as computed in [MCEAC94] program double z; for (n = 0; n < N/2; n++) { z = (.5+(N/2)-1.-n)*(twopi/2.)/N; w[n] = exp(-c*z*z); w[N-1-n] = w[n]; } //multiply input data by window points; for (n = 0; n < N; n++) { r[n] = w[n]*r[n]; // we're overwriting r[n] with the windowed r[n] } //end for fft_recur(N, r, i); // note: fft is not scaled // compute amplitude squared of bins 7 and 8 r[7] = r[7]*r[7] + i[7]*i[7]; // r[7] contains amplitude squared r[8] = r[8]*r[8] + i[8]*i[8]; // r[8] contains amplitude squared // compute natural log of FFT amplitude output bins 7 and 8 r[7] = .5*log( r[7] ) ; // r[7] is now the natural log of the amplitude of FFT bin 7 r[8] = .5*log( r[8] ) ; // .5*log(A squared) = log(A) cout<<"\n\nfrequency and amplitude estimation results\n\n frequency amplitude\n\n"; fbin = 7+.5+.5*c*(r[8]-r[7]); F = (sample_rate/N)*fbin; amp = sqrt(c)*exp(r[7]+((fbin-7)*(fbin-7))/c)*(2.0*sqrt(twopi/2.)/N) ; printf("%10.6f\t%10.6f\n",F,amp); system ("PAUSE"); return 0; } // end main //******************** fft_recur *********************** void fft_recur(long N, double *r, double *i) { long h, i1, j = 0, k, i2 = N/2, l1, l2 = 1; double c = -1.0, s = 0.0, t1, t2, u1, u2; for (h = 0; h < N-2; h++) { // ***** bit reverse starts here *** if (h < j) { t1 = r[h]; r[h] = r[j]; r[j] = t1; t2 = i[h]; i[h] = i[j]; i[j] = t2; } k = i2; while (k <= j) { j = j - k; k = k/2; } j = j + k; } //****** bit reverse done ****** for (k = 1; k < N; k = k*2) { l1 = l2; l2 = l2*2; u1 = 1.0; u2 = 0.0; for (j = 0; j < l1; j++) { for (h = j; h < N; h = h + l2) { i1 = h + l1; t2 = (r[i1] - i[i1])*u2 ; t1 = t2 + r[i1]*(u1 - u2) ; t2 = t2 + i[i1]*(u1 + u2) ; r[i1] = r[h] - t1; i[i1] = i[h] - t2; r[h] = r[h] + t1; i[h] = i[h] + t2; } // end for over h t1 = u1 * c - u2 * s; u2 = u1 * s + u2 * c; u1 = t1; //x = u1 - u2; y = u1 + u2; } // end for over j s = - sqrt((1.0 - c) / 2.0); c = sqrt((1.0 + c) / 2.0); } // end for over k } // end fft The above program does the following: 1) collect 32 points 2) use a Gaussian window on the data (note 1: the variable 'c' in the program determines the beamwidth of the Gaussian window - it is a critical and sometimes sensitive value, and it can range from about 2.5 to 4 - sometimes higher or lower depending on what you're trying to do and what you're estimating - a value of about 3 is used above). (note 2: the Gaussian window shown above is actually not a correct window for use with an FFT - which would require that the zero point be unique - I have derived a corrected version of this window, but the one shown here is adequate as long as N isn't too small) 3) take 32 point FFT of the windowed data (not the most efficient way, since you only need bins 7 and 8, but it gives you flexibility if you have to change any of your numbers, such as sample rate or N, or you're estimating different frequencies than the ones you've mentioned). 4) compute amplitude squared of bins 7 and 8 (r squared plus i squared) 5) compute the natural log of bins 7 and 8: 6) estimate the frequency 7) estimate the amplitude The results of the above program are: ******************************************* frequency and amplitude estimation results frequency amplitude 2222.122610 0.999918 ****************************************** Change the frequency in the program to 2475, and the results are: ****************************************** frequency and amplitude estimation results frequency amplitude 2474.863271 0.999900 ***************************************** These results are ... well, terrible. If I tweaked things a little, and placed my bins more optimally, I'd get better results - at least 6+ decimal digits of accuracy - and that's only if I'm not trying too hard. With no added noise and careful computation and programming, the accuracy of the results could be limited by the noise floor created by the FFT. I was going to do some tweaking but then I thought: what if my understanding of what you're doing is not correct? Do those frequencies of yours run immediately one after another (ie: 2222 goes on for 16 ms and then abruptly changes to 2475, which goes on for 16 ms, and so on?). That's going to create a problem, since the estimates are going to be stable for awhile, then get mixed up, and back to stable again. You'd probably need to establish a threshold for when the frequency or amplitude estimates get too far out of range. And how good is the frequency control of whatever is generating them? (can the frequencies vary by .1%, 1%, 5% ??). Or might they be stable, but you're not certain of their exact frequency. I've used this technique for over 20 years now, and I know it well. And I also know that's it's very easy for a beginner to mess up, because they're not aware of the pitfalls (eg: operating too near the 0 or N/2 bins, not having enough cycles of the input waveform, not using the optimum Gaussian beamwidth, or the wrong equation for the window, improper sampling rate or wrong N for the problem at hand, etc, etc). But it works well for those who know what they're doing and have some experience with it. Kevin McGee
Reply by ●October 17, 20142014-10-17
Thank you very much for your time and your explanations. Le 17/10/2014 02:25, kevinjmcee a �crit :> On Wed, 15 Oct 2014 17:51:01 +0200, Lanarcam wrote: > (snip) > > First, I agree with the other comments that a 5K sample rate is too low. 10K might be more reasonable. > > If the actual carrier frequency is slightly different than what you expect, then the first thing may want to do is estimate it. You could use Gaussian ratio estimation to do so, as described here: > > https://www.edn.com/design/analog/4352503/EDN--03-03-94-Ratio-detection-precisely-characterizes-signals-amplitude-and-frequenc > > http://www.ptodirect.com/Results/Patents?query=PN/5214708 > > (note: the patent was issued in 1993, so it's more than 20 years old, and besides, it's for speech) > > Basically, you window your input data with a Gaussian window, then do an FFT, after which you take the log of the FFT amplitudes. As it turns out, when you use a window, you're creating a series of band-pass filters in the frequency domain, where each filter is centered at an FFT bin. > > Now if you have a sinusoid in the input that lies between two consecutive bins of an FFT, then it will contribute to (at least) 2 of the FFT bins - mostly to the one that's higher or lower in frequency. Well, it turns out that you can write 2 equations with 2 unknowns (amplitude and frequency). By using the log values from the FFT, and taking the ratio of the 2 equations, the amplitude drops out, and you can solve for frequency. > > There is another, slightly different Gaussian method that's used at CERN, but it's a little more complicated, and you probably don't need that. > > Incidentally, CERN has done a number of detailed accuracy studies on their method - after all, they use it for estimating beam frequencies in a $10 billion supercollider. > > Now, given your numbers, and presuming a sample rate (s_r) of 10K and a 32 point FFT (N = 32), the 2 bin numbers bracketing your frequencies are: > > k1 = s_r/N = f1 = 2187.5 --> k1 = 7 > k2 = s_r/N = f2 = 2500.0 --> k2 = 8 > > So you'd collect 32 samples, take the FFT, compute the log of bins 7 and 8, then put those numbers into an equation which will give you the frequency estimate. Then you can use a second equation to compute amplitude. > > I originally wrote this post as a step by step sequence, but, since you're a beginner, you might need something a bit more detailed. So I cobbled together a C++ program to demonstrate (I used Ariel 10 format, so I hope it prints OK): > > #include <cstdio> > #include <cstdlib> > #include <iostream> > #include <cmath> > using namespace std; > void fft_recur(long, double* r, double* i); // prototype declaration > int main (int nNumberofArgs, char* pszArgs[ ] ) { > const long N = 32; double sample_rate = 10000., r[N] = {0.}, i[N] = {0.}, w[N] = {0.} ; > long n; double t, twopi = 6.2831853071795865; > double F, fbin, amp, power, c = 3.0; > > for (n = 0; n < N; n++) { // generate test data > t = twopi*n/sample_rate; > r[n] = 1.*cos(2222.0*t + twopi/5.) ; // amplitude = 1. frequency = 2222., arbitrary phase > } // end for > > // Gaussian window as computed in [MCEAC94] program > double z; > for (n = 0; n < N/2; n++) { > z = (.5+(N/2)-1.-n)*(twopi/2.)/N; > w[n] = exp(-c*z*z); > w[N-1-n] = w[n]; > } > > //multiply input data by window points; > for (n = 0; n < N; n++) { > r[n] = w[n]*r[n]; // we're overwriting r[n] with the windowed r[n] > } //end for > > fft_recur(N, r, i); // note: fft is not scaled > > // compute amplitude squared of bins 7 and 8 > r[7] = r[7]*r[7] + i[7]*i[7]; // r[7] contains amplitude squared > r[8] = r[8]*r[8] + i[8]*i[8]; // r[8] contains amplitude squared > > // compute natural log of FFT amplitude output bins 7 and 8 > r[7] = .5*log( r[7] ) ; // r[7] is now the natural log of the amplitude of FFT bin 7 > r[8] = .5*log( r[8] ) ; // .5*log(A squared) = log(A) > > cout<<"\n\nfrequency and amplitude estimation results\n\n frequency amplitude\n\n"; > > fbin = 7+.5+.5*c*(r[8]-r[7]); > F = (sample_rate/N)*fbin; > amp = sqrt(c)*exp(r[7]+((fbin-7)*(fbin-7))/c)*(2.0*sqrt(twopi/2.)/N) ; > printf("%10.6f\t%10.6f\n",F,amp); > > > system ("PAUSE"); > return 0; > } // end main > //******************** fft_recur *********************** > void fft_recur(long N, double *r, double *i) { > long h, i1, j = 0, k, i2 = N/2, l1, l2 = 1; > double c = -1.0, s = 0.0, t1, t2, u1, u2; > > for (h = 0; h < N-2; h++) { // ***** bit reverse starts here *** > if (h < j) { > t1 = r[h]; r[h] = r[j]; r[j] = t1; > t2 = i[h]; i[h] = i[j]; i[j] = t2; > } > k = i2; > while (k <= j) { > j = j - k; k = k/2; > } > j = j + k; > } //****** bit reverse done ****** > > for (k = 1; k < N; k = k*2) { > l1 = l2; l2 = l2*2; > u1 = 1.0; u2 = 0.0; > for (j = 0; j < l1; j++) { > for (h = j; h < N; h = h + l2) { > i1 = h + l1; > t2 = (r[i1] - i[i1])*u2 ; > t1 = t2 + r[i1]*(u1 - u2) ; > t2 = t2 + i[i1]*(u1 + u2) ; > r[i1] = r[h] - t1; > i[i1] = i[h] - t2; > r[h] = r[h] + t1; > i[h] = i[h] + t2; > } // end for over h > t1 = u1 * c - u2 * s; > u2 = u1 * s + u2 * c; > u1 = t1; //x = u1 - u2; y = u1 + u2; > } // end for over j > s = - sqrt((1.0 - c) / 2.0); > c = sqrt((1.0 + c) / 2.0); > } // end for over k > > } // end fft > > > The above program does the following: > > 1) collect 32 points > > 2) use a Gaussian window on the data > (note 1: the variable 'c' in the program determines the beamwidth of the Gaussian window - it is a critical and sometimes sensitive value, and it can range from about 2.5 to 4 - sometimes higher or lower depending on what you're trying to do and what you're estimating - a value of about 3 is used above). > (note 2: the Gaussian window shown above is actually not a correct window for use with an FFT - which would require that the zero point be unique - I have derived a corrected version of this window, but the one shown here is adequate as long as N isn't too small) > > 3) take 32 point FFT of the windowed data (not the most efficient way, since you only need bins 7 and 8, but it gives you flexibility if you have to change any of your numbers, such as sample rate or N, or you're estimating different frequencies than the ones you've mentioned). > > 4) compute amplitude squared of bins 7 and 8 (r squared plus i squared) > > 5) compute the natural log of bins 7 and 8: > > 6) estimate the frequency > > 7) estimate the amplitude > > The results of the above program are: > ******************************************* > frequency and amplitude estimation results > > frequency amplitude > > 2222.122610 0.999918 > ****************************************** > > Change the frequency in the program to 2475, and the results are: > ****************************************** > frequency and amplitude estimation results > > frequency amplitude > > 2474.863271 0.999900 > ***************************************** > > These results are ... well, terrible. If I tweaked things a little, and placed my bins more optimally, I'd get better results - at least 6+ decimal digits of accuracy - and that's only if I'm not trying too hard. With no added noise and careful computation and programming, the accuracy of the results could be limited by the noise floor created by the FFT. > > I was going to do some tweaking but then I thought: what if my understanding of what you're doing is not correct? Do those frequencies of yours run immediately one after another (ie: 2222 goes on for 16 ms and then abruptly changes to 2475, which goes on for 16 ms, and so on?). That's going to create a problem, since the estimates are going to be stable for awhile, then get mixed up, and back to stable again. You'd probably need to establish a threshold for when the frequency or amplitude estimates get too far out of range. > > And how good is the frequency control of whatever is generating them? (can the frequencies vary by .1%, 1%, 5% ??). Or might they be stable, but you're not certain of their exact frequency. > > I've used this technique for over 20 years now, and I know it well. And I also know that's it's very easy for a beginner to mess up, because they're not aware of the pitfalls (eg: operating too near the 0 or N/2 bins, not having enough cycles of the input waveform, not using the optimum Gaussian beamwidth, or the wrong equation for the window, improper sampling rate or wrong N for the problem at hand, etc, etc). But it works well for those who know what they're doing and have some experience with it. > > Kevin McGee >
Reply by ●October 27, 20142014-10-27
On Thu, 16 Oct 2014 17:25:23 -0700 (PDT), kevinjmcee <kevinjmcgee@netscape.net> wrote:>On Wed, 15 Oct 2014 17:51:01 +0200, Lanarcam wrote: >(snip) > >First, I agree with the other comments that a 5K sample rate is too low. 10K might be more reasonable. >[Snipped by Lyons> >Kevin McGeeHi Kevin, I sent you a private e-mail yesterday [on 10/26/14]. Did you receive it? Thanks, [-Rick-]
