Do you know this phenomenon in DFT? 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 &lt;iostream&gt;
#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 &lt;= M; ++m)
{
testFreq = lowTestFreq + m * testFreqStep; // Hz.
testUnitAngle = (testFreq / samplFreq) * (2 * pi); // rad/sample.

Re = 0; Im = 0;

for (n = 1; n &lt;= 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 &lt;&lt; "Test Freq " &lt;&lt; testFreq; // Hz.
 PowerSpectrum = Re * Re + Im * Im; //pure number.       
 std::cout &lt;&lt; " MagnSquared " &lt;&lt; PowerSpectrum &lt;&lt; std::endl;
}       
return 0;
}
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.