THD+N measure without knowing the fundamental frequency

[Daniel Coelho]

Hi all,

I have a NI PCI 4461 and I want to measure the Vrms, THD and THD+N from an acquired signal.

I managed to measure de Vrms with CycleRMSAverage().

I'm using HarmonicAnalyzer() to obtain the THD and the THD+N but the fundamental frequency is asked and I want my code to find it out by itself, i. e., how do I find the fund. freq. and how do I measure the THD+N with it?

I'm trying to use the highest peak and its corresponding freq. on the AutoPowerSpectrum() to determine the fundamental freq. but I'm haven't managed to do it yet.

I apreciate your time.

Thanks in advance.

Daniel Coelho

Portugal

[Daniel Coelho]

I managed to find the fundamental frequency using the highest peak on the AutoPowerSpectrum() function.

But HarmonicAnalyzer() returns bad values for THD and THD+N.

What could be the problem?

thanks in advance.

Daniel Coelho

[miniMe]

Hello Daniel,

You can use following function to find the peak frequency:

AnalysisLibErrType PowerFrequencyEstimate (double autoSpectrum[], int numberOfElements, double searchFreq, WindowStruct windowConstants, double df, int freqSpan, double *freqPeak, double *powerPeak);

If you set searchFreq to -1, the function will estimate the frequency peak and power.

You can then feed the peak frequency to the HarmonicAnalyzer Function.

Regards

[Daniel Coelho]

Thanks miniMe,

you're correct. That function is more acurate to find the peak.

But a problem remains, I use the HarmonicAnalyser() to find THD and THD+N but the values are absurd.

My sampling frequency is more than twice as big as the fundamental frequency and I use a big number for samples.

It's an almost perfect sinewave because I see it on my Scopemeter and the sinewave is generated by a Signal Generator.

Also, some of the Harmonics on the harmonic_amplitudes array have NaN values.

Hope someone can help me.

Thanks for your help.

Daniel Coelho

Portugal

[miniMe]

Could you post some of the code ?

[Daniel Coelho]

...

DAQmxReadAnalogF64(taskHandle,samples,10.0,DAQmx_Val_GroupByChannel,data,samples*numChannels,&numRead,NULL);

...

CycleRMSAverage (&(data[i*numRead]), numRead, 1, 100, 50, 0, 1, STATE_LEVELS_AUTO_SELECT, 1, &CycleAverage, &Vrms1, &StartTime,  &EndTime, &lowRefLevel, &midRefLevel, &highRefLevel);

(not very sure about ", 100, 50, 0,")

...

auto_power_spectrum = (double *)malloc((samples/2)*sizeof(double));
  
   AutoPowerSpectrum (data, samples/2, 1.0/(sampFreq), auto_power_spectrum, &df);
  
   harmonic_amplitudes = (double *)malloc(number_of_harmonics * sizeof(double));

   harmonic_frequencies = (double *)malloc(number_of_harmonics * sizeof(double));
                                    
   HarmonicAnalyzer(auto_power_spectrum, samples/2, samples, number_of_harmonics,
                             0, sampFreq, fundFreq, harmonic_amplitudes, harmonic_frequencies, &thd1, &thdnoise1);

...

My fundamental frequency is *usually* 1500 Hz, and I use 4000 samples and 55 000 S/s on the sampling frequency.

But THD and THD+N have values like 1000000000 and NaN and change a lot. I should be getting somethin like 0.2 for THD (value obtained from LabVIEW).

Hope this is enough.

Tnhks 😉

[miniMe]

Have you tried to visualize the waveform array after the sampling ? You should be able to see a clearly visible sinewave - if input signal has enough amplitude.

I tested Harmonic Analyzer function with following code:

#define HARM_NO 9

#define POINTS 8192

double thd, thd_noise,df;

double harm_frq_array[HARM_NO];

double harm_amplitude_array[HARM_NO];

WindowConst constants;

double waveform_array[POINTS];

double spectrum_array[POINTS/2];

double converted_spectrum_array[POINTS/2];

double waveform_amplitude, frequency=10;

double amplitude_peak,frequency_peak;

//...

SineWave (POINTS, waveform_amplitude, frequency/(double)POINTS, &phase, waveform_array); // Generate a sinewave

ScaledWindow (waveform_array, POINTS, 0, &constants);

AutoPowerSpectrum (waveform_array, POINTS, 1 /(double) POINTS, spectrum_array,&df);

PowerFrequencyEstimate (spectrum_array, POINTS/2, -1.0, constants, df, 7, &frequency_peak, &amplitude_peak);

error_code = HarmonicAnalyzer (spectrum_array, POINTS/2, 0, HARM_NO, 2, POINTS, frequency_peak, harm_amplitude_array, harm_frq_array, &thd, &thd_noise);

//...

Regards

[Daniel Coelho]

Thank you for your quick reply.

I checked your code and, except on one or two arguments, everything looks OK.

I'm starting to believe that my Digital Function Generator is broken. LabVIEW is getting abnormal results too.

I'll try to rewrite the code and change the Generator.

I hope I can give good news soon.

Thnks for your help.

If I manage to correct this I'll post my conclusions, ifnot, I'll ask for your help again 😛

Daniel Coelho

[Daniel Coelho]

status = DAQmxCreateTask("",&taskHandle);
 if (status != 0)
 {
  strcpy(error, "Erro no Init\0");
  return status; 
 }
 strcpy(chan1, DeviceName);       
 strcat(chan1, Port);
 strcat(chan1, "/ai0\0");
 strcpy(chan2, DeviceName);
 strcat(chan2, Port);
 strcat(chan2, "/ai1\0");             // the purpose of this is to write Dev2/ai0, for instance
 
 DAQmxCreateAIVoltageChan(taskHandle,chan1,"",DAQmx_Val_Cfg_Default,min,max,DAQmx_Val_Volts,NULL);
 DAQmxCreateAIVoltageChan(taskHandle,chan2,"",DAQmx_Val_Cfg_Default,min,max,DAQmx_Val_Volts,NULL);
 DAQmxGetTaskAttribute(taskHandle,DAQmx_Task_NumChans,&numChannels);
 if (numChannels <2)
 {
  strcpy(error, "Erro no Init\0");
  return -1; 
 }

....

 DAQmxCfgSampClkTiming(taskHandle,"",sampFreq,DAQmx_Val_Rising,DAQmx_Val_FiniteSamps,samples);
 if( (data=malloc(samples*numChannels*sizeof(double)))==NULL )
 {
  MessagePopup("Error","Memoria insuficiente para guardar dados\0");
  return (-1);
 }
 DAQmxStartTask(taskHandle);
 DAQmxReadAnalogF64(taskHandle,samples,10.0,DAQmx_Val_GroupByChannel,data,samples*numChannels,&numRead,NULL);
 
 // ********************Vrms
 if( numRead>0 )
 {
  for(i=0;i<numChannels;i++)
  {
   if(i == 0) //channel 0
    CycleRMSAverage (&(data[i*numRead]), numRead, 1, 100, 50, 0, 1, STATE_LEVELS_AUTO_SELECT, 1, &CycleAverage, &Vrms1, &StartTime,
      &EndTime, &lowRefLevel, &midRefLevel, &highRefLevel);
   if(i == 1) // channel 1
    CycleRMSAverage (&(data[i*numRead]), numRead, 1, 100, 50, 0, 1, STATE_LEVELS_AUTO_SELECT, 1, &CycleAverage, &Vrms2, &StartTime,
      &EndTime, &lowRefLevel, &midRefLevel, &highRefLevel);
  }

 }
 else
 {
  strcpy(error, "Erro no numero de pontos\0");
  return 0;
 }

 ScaledWindow (data, samples, 0, &WinConst);
 auto_power_spectrum = (double *)malloc((samples/2)*sizeof(double));
 AutoPowerSpectrum (data, samples/2, 1.0/(sampFreq),auto_power_spectrum, &df);
 // ********************FundFreq
 PowerFrequencyEstimate (auto_power_spectrum, samples/2, -1, WinConst, df, 7, &Freq_Peak, &Power_Peak);
 Converted_Spec = calloc (samples/2, sizeof(double));
 SpectrumUnitConversion (auto_power_spectrum, samples/2, 0, 0, 0, df, WinConst, Converted_Spec, Unit);
 number_of_harmonics = sampFreq/(2 * Freq_Peak);
 if ( number_of_harmonics == 0)
  number_of_harmonics = 5;
 harmonic_amplitudes = (double *)malloc(number_of_harmonics * sizeof(double));
    harmonic_frequencies = (double *)malloc(number_of_harmonics * sizeof(double));
 // ********************THD e THD+N
    HarmonicAnalyzer(auto_power_spectrum, samples/2, samples, number_of_harmonics,
                     2  , sampFreq, Freq_Peak, harmonic_amplitudes,
      harmonic_frequencies, &thd1, &thdnoise1);
 
 if (strcmp(KeyCode, "GETVRMS") == 0)
  sprintf(Response, "%f", Vrms1);
 if (strcmp(KeyCode, "GETTHD") == 0)
  sprintf(Response, "%f", thd1);
 if (strcmp(KeyCode, "GETTHDNOISE") == 0)
  sprintf(Response, "%f", thdnoise1);
 if (strcmp(KeyCode, "GETFUNDFREQ") == 0)
  sprintf(Response, "%f", Freq_Peak);
 if (strcmp(KeyCode, "GETFAMPL") == 0)
  sprintf(Response, "%f", amplitude);
 free(data);
 free(auto_power_spectrum);
 free(harmonic_amplitudes);
 free(harmonic_frequencies);
 free(Converted_Spec);

...

if( taskHandle!=0 )
 {
  DAQmxStopTask(taskHandle);
  DAQmxClearTask(taskHandle);
 }

This is the code I m using, more or less.

I've seen on a Graph the sine I'm generating and it seems fine.

I've changed my Digital Function Generator and THD (on LabVIEW) is nicer but LabWindowsCVI still gives me bad values.

Is there a problem with LabWindows or am I doing something wrong?

I apreciate all the help.

Daniel Coelho

[miniMe]

Daniel,

I think your problem with Harmonic Analyzer is due to a DC component on the sampled waveform.

Please check if your sinewave is simetrical in respect to 0.

If the sinewave has a DC offset, Harmonic Analyzer might return a THD value of 1000000000.00 (this value I get for a 0.25 amplitude sinewave with a 1.0 DC offset).

If your waveform has a mean value other than 0,  try to remove the DC component before ScaledWindow call.

Regards