4
$\begingroup$

Please let me ask you about a phenomenon in DFT. Below DFT program outputs excellent results for low sample rates, but for high results are very bad.

Specifically: for sample rate 50 samples/sec, segment 1000 samples, existing frequency 36.7634 Hz, grid from 35.2364 Hz to 40.7896 divided by 40, outputs maximum magnitude squared 1.22629e+15 corresponding to testing frequency 36.7635.

But for sample rate 50,000 samples/sec, segment also 1000 samples, and same existing frequency and grid, outputs maximum magnitude squared 8.46881e+14 corresponding to testing frequency 38.2907 Hz. Error is 5 significant digits.

Are you aware of this phenomenon? If yes please explain me that.

// Grid is created by variables lowTestFreq and highTestFreq
// both of double type, and integer type M.

#include <iostream>
#include "math.h"

int main()
{  
//Unchanged variables.
constexpr double pi = 3.141592653589793; //pure number.
int n = 0; // sample.
double xn = 0; //pure number.
double Re = 0; //pure number.
double Im = 0; //pure number.
double xn_cos = 0; //pure number.
double xn_sin = 0; //pure number.
double PowerSpectrum =0; //pure number. Name given by my dear dsp Teacher Mr. Lyons.

//Sampling frequency and samples of segment.
int samplFreq = 50000; // sample/second.
int N = 1000; //sample.

// 1st existing sinusoidal.
double ampl_1 = 70000; //pure number.
double existFreq_1 = 36.7634; // Hz.
double phase_1 = 0.67 * pi; // rad.
double unitAngle_1 = (existFreq_1 / samplFreq) * (2 * pi); // rad/sample.

// 2nd existing sinusoidal.
double ampl_2 = 60000; // pure number.
double existFreq_2 = 4505.75; // Hz.
double phase_2 = -0.37 * pi; // rad.
double unitAngle_2 = (existFreq_2 / samplFreq) * (2 * pi); // rad/sample.

// Testing grid.
double testFreq = 0; // Hz.
double testUnitAngle = 0; // rad/sample.
double testAngle = 0; // rad.
double lowTestFreq = 35.2364; //Hz.
double highTestFreq = 40.7896; //Hz.
int M = 40; // pure number.
int m = 0; // pure number.
double testFreqStep = (highTestFreq - lowTestFreq) / M; // Hz.

for (m = 0; m <= M; ++m)
{
testFreq = lowTestFreq + m * testFreqStep; // Hz.
testUnitAngle = (testFreq / samplFreq) * (2 * pi); // rad/sample.

Re = 0; Im = 0;

for (n = 1; n <= N; ++n)
{
xn = ampl_1 * sin(n * unitAngle_1 + phase_1)
   + ampl_2 * sin(n * unitAngle_2 + phase_2);

testAngle = n * testUnitAngle;
xn_cos = xn *  cos(testAngle);
xn_sin = xn * -sin(testAngle);

Re += xn_cos;
Im += xn_sin;
}  
 std::cout << "Test Freq " << testFreq; // Hz.
 PowerSpectrum = Re * Re + Im * Im; //pure number.       
 std::cout << " MagnSquared " << PowerSpectrum << std::endl;
}       
return 0;
}
$\endgroup$
  • 1
    $\begingroup$ Have you verified your implementation with a more flexible language such as python/numpy or matlab using the builtin fft function? Can you reproduce the problem there? If so, I think we would be able to help better, if we get the code in these more concise languages. $\endgroup$ – Maximilian Matthé Mar 25 '17 at 15:43
  • $\begingroup$ @ Maximilian Matthé: Unfortunately I have not, and I can not reproduce the problem there. Many thanks for you are interested in helping me. $\endgroup$ – George Theodosiou Mar 25 '17 at 15:48
  • $\begingroup$ in this case I would suggest you get a reference implementation and compare your intermediate results step by step wtih the reference to see where the problem is. $\endgroup$ – Maximilian Matthé Mar 26 '17 at 8:24
  • $\begingroup$ @ Maximilian Matthé: I tested it with sample rate 12,000 samples / sec, existing frequency 51 Hz, grid from 50 Hz to 52 divided by 1000 (M = 1000). Maximum magnitude squared corresponds to 51.01 Hz. I regard it sufficient. In general I can trade sample rate off for lowest existing / testing frequencies. My question was theoretical, just to learn the cause of this phenomenon. But never mind! Regards. $\endgroup$ – George Theodosiou Mar 28 '17 at 12:04
  • $\begingroup$ This is most likely an artifact of (the rectangular) windowing that results in amplitude distortion and can cause interference between the two sine waves. If I am right, it isn't high frequency versus low frequency that matters as much as it is the misalignment of the frequencies with the center frequencies of the DFT bins. Using a longer window should ameliorate this problem. Another way to avoid it would be to use coherent sampling (make sure the sampling period is commensurate with the sinusoid periods). $\endgroup$ – hops May 2 '17 at 16:07
1
$\begingroup$

If you are sampling 50 samples/sec your nyquist frequency is only 25Hz, so your test signal is being aliased. You should expect different results if your sampling rate is less than double the highest frequency in your signal.

$\endgroup$
  • $\begingroup$ @ Bumzur. Many thanks for you answered my question. I tested with sample rate 100 samples / sec and results are same. Regards. $\endgroup$ – George Theodosiou Mar 28 '17 at 8:32
  • 2
    $\begingroup$ @GeorgeTheodosiou, why don't you bump up the score on Bum's answer and also select it as the accepted answer? show the Bum a little appreciation. $\endgroup$ – robert bristow-johnson Apr 25 '17 at 6:36
  • $\begingroup$ @robert bristow-johnson: Gentleman, I appreciate Mr. Bumzur's answer, but it does not answer my question. As I commented his answer, results are same when sample rate is 100 samples / sec, so existing frequency (36.7634 Hz) is lower than Nyquist's. Many thanks for your comment. Regards. $\endgroup$ – George Theodosiou Apr 26 '17 at 7:09
0
$\begingroup$

You take 1000 samples at 50000 samples/sec. That's a sampling period of 20ms. For a frequency of 36.7634Hz (which has a period of 27ms more or less). So you sample only a partial period and expect useful statistics?

Note that your slow sampling rate more or less has its sampling points distributed randomly across the different phases. It's not really random, but at least you get a reasonably even coverage of different phases as long as your sampling interval is not too closely related to the periodicity of the signal.

All of that has very little to do with DFT, just with what you are actually sampling. Basically it is the garbage-in/garbage-out principle at work.

$\endgroup$
  • $\begingroup$ 20ms = 0.02s => 50 samples/sec. $\endgroup$ – Peter K. Apr 25 '17 at 11:26
  • $\begingroup$ @user28163: Gentleman, please accept my many thanks for your answer. As Mr. Peter K. commented, you did a common human error. For sample rate 50 ksample / sec, sampling period is 20 μsec. By trial and error I found that with 12 ksamples / sec, and 800 samples segment, results are acceptable even for 25 Hz frequency. My question was theoretical. If anybody knows the reason. Regards. $\endgroup$ – George Theodosiou Apr 26 '17 at 7:26
0
$\begingroup$

First thing you should do is change all your floats to doubles if you are having precision problems. If that doesn't fix it, I'll take a closer look at your code.

For example, running numbers in the 1e14 range is definitely going to cause precision loss in a float. It'll just barely fit in a double.

Ced

$\endgroup$
  • $\begingroup$ @ Cedron Dawg, Gentleman please accept my many thanks for you answered my questions. I did your instruction. Problem is the ratio testing frequencies to sampling rate. It demands compromise. Regards. $\endgroup$ – George Theodosiou Feb 7 '18 at 10:24
0
$\begingroup$

Gentlemen, answerers and commentators,

Please let me express my sincere gratitude for all your help. I'm interested in create program (project) catching, with the help of dsp kit, in real time, fundamental frequency (first harmonic, characterizing sound's pitch) of chanter's or singer's voice, a periodic signal. Range of voice signal I'm interested for, is 25 Hz - 5 kHz, that is piano (keyboard) range and little more. I expect catching it by "nude" DFT, performing operator ">" on magnitude square of DFT complex numbers output. First (local) maximum should be first harmonic's.

Because electronic (analog) filter distorts signal severely, I expect, instead, perform 96 ksamples/sec sample rate. Then, only frequencies in the range 91 kHz - 96 kHz can duplicate in interesting range. But such frequencies are negligible in voice. On the other hand with so high sample rate results from DFT are BAD, apparently due to high ratio (higher than 1000) sample rate/frequency. However by trial and error I found, for sample rate 12 ksample/sec and segment 800 samples, results are GOOD for frequencies higher than 100 Hz, and with 4 ksample/sec are GOOD for 25 Hz - 100 Hz and more. I asked this question for learn the theoretically right way (algorithm or formula or table) for get optimum sample rate and segment. One month has passed but did not learn that. Then I conclude perform what found by trial and error.

For implement that, I expect perform digital low pass FIR filter with cut-off zone 5 kHz - 6 kHz. Then, divide sample rate by 8, and perform DFT with 800 samples segment. If fundFreq is lower than 400 Hz, perform second low pass FIR filter with cut-off zone 500 Hz - 2000 Hz, divide sample rate further more by 3, and perform DFT. Example C++ program is following.

/* Outputs 1st (fundamental) frequency.
1. You can change values.
2. You can add more sinusoids, but They should be multiples of first existing.
3. Low testing frequency be lower than first existing, and high be midway 1st and 2nd existing or lower then midway but higher than first.
Grid is created by variables lowTestFreq and highTestFreq both of double type, and integer type M.*/

#include <iostream>
#include "math.h"

int main()
{
//Unchanged variables.
constexpr double pi = 3.141592653589793; //pure number.
int n = 0; // sample.
double xn = 0; //pure number.
double Re = 0; //pure number.
double Im = 0; //pure number.
double xn_cos = 0; //pure number.
double xn_sin = 0; //pure number.
double testPowerSpectrum =0; //pure number.
double firstHarmonic = 0;
double maxPowerSpectrum = 0;

//Sampling frequency and samples of segment.
int samplFreq = 8000; // sample/second.
int N = 1200; //sample.

// 1st existing sinusoidal.
double ampl_1 = 70000; //pure number.
double existFreq_1 = 256; // Hz.
double phase_1 = 0.67 * pi; // rad.
double unitAngle_1 = existFreq_1 * 2 * pi / samplFreq; // rad/sample.

// 2nd existing sinusoidal.
double ampl_2 = 600000; // pure number.
double existFreq_2 = 2 * existFreq_1; // Hz.
double phase_2 = -0.37 * pi; // rad.
double unitAngle_2 = existFreq_2 * 2 * pi / samplFreq; // rad/sample.

// 3rd existing sinusoidal.
double ampl_3 = 100000; // pure number.
double existFreq_3 = 3 * existFreq_1; // Hz.
double phase_3 = -0.73 * pi; // rad.
double unitAngle_3 = existFreq_3 * 2 * pi / samplFreq; // rad/sample.

// Testing grid.
double testFreq = 0; // Hz.
double testUnitAngle = 0; // rad/sample.
double testAngle = 0; // rad.
double lowTestFreq = 254; //Hz.
double highTestFreq = 258; //Hz.
int M = 400; // pure number.
int m = 0; // pure number.
double testFreqStep = (highTestFreq - lowTestFreq) / M; // Hz.

for (m = 0; m <= M; ++m)
{
testFreq = lowTestFreq + m * testFreqStep; // Hz.
testUnitAngle = testFreq * 2 * pi / samplFreq; // rad/sample.

Re = 0; Im = 0;

for (n = 1; n <= N; ++n)
{
xn = ampl_1 * sin(n * unitAngle_1 + phase_1)
+ ampl_2 * sin(n * unitAngle_2 + phase_2)
+ ampl_3 * sin(n * unitAngle_3 + phase_3);

testAngle = n * testUnitAngle;
xn_cos = xn *  cos(testAngle);
xn_sin = xn * -sin(testAngle);

Re += xn_cos;
Im += xn_sin;
}
testPowerSpectrum = Re * Re + Im * Im; //pure number.
if (testPowerSpectrum > maxPowerSpectrum)
{
maxPowerSpectrum = testPowerSpectrum;
firstHarmonic = testFreq;
}
}
std::cout << "First harmonic: " << firstHarmonic << " Hz" << "\n";
return 0;
}

I remain with regards and friendship, remembering that in other questions of mine, you answered helpfully.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.