thd

Total harmonic distortion

Syntax

r = thd(x)

r = thd(x,fs,n)

r = thd(pxx,f,"psd")

r = thd(pxx,f,n,"psd")

r = thd(sxx,f,rbw,"power")

r = thd(sxx,f,rbw,n,"power")

r = thd(___,"aliased")

[r,harmpow,harmfreq] = thd(___)

thd(___)

Description

r = thd(x) returns the total harmonic distortion (THD) in dBc of the real-valued sinusoidal signal x. The total harmonic distortion is determined from the fundamental frequency and the first five harmonics using a modified periodogram of the same length as the input signal. The modified periodogram uses a Kaiser window with β = 38.

r = thd(x,fs,n) specifies the sample rate fs and the number of harmonics (including the fundamental) to use in the THD calculation.

r = thd(pxx,f,"psd") specifies the input pxx as a one-sided power spectral density (PSD) estimate. f is a vector of frequencies corresponding to the PSD estimates in pxx.

r = thd(pxx,f,n,"psd") specifies the number of harmonics (including the fundamental) to use in the THD calculation.

r = thd(sxx,f,rbw,"power") specifies the input as a one-sided power spectrum. rbw is the resolution bandwidth over which each power estimate is integrated.

r = thd(sxx,f,rbw,n,"power") specifies the number of harmonics (including the fundamental) to use in the THD calculation.

r = thd(___,"aliased") reports harmonics of the fundamental that are aliased into the Nyquist range. Use this option when the input signal is undersampled. If you do not specify this option, or if you set it to "omitaliases", then the function ignores any harmonics of the fundamental frequency that lie beyond the Nyquist range.

[r,harmpow,harmfreq] = thd(___) returns the powers (in dB) and frequencies of the harmonics, including the fundamental.

thd(___) with no output arguments plots the spectrum of the signal and annotates the harmonics in the current figure window. It uses different colors to draw the fundamental component, the harmonics, and the DC level and noise. The THD appears above the plot. The fundamental and harmonics are labeled. The DC term is excluded from the measurement and is not labeled.

Examples

Determine THD for a Signal with Two Harmonics

This example shows explicitly how to calculate the total harmonic distortion in dBc for a signal consisting of the fundamental and two harmonics. The explicit calculation is checked against the result returned by thd.

Create a signal sampled at 1 kHz. The signal consists of a 100 Hz fundamental with amplitude 2 and two harmonics at 200 and 300 Hz with amplitudes 0.01 and 0.005. Obtain the total harmonic distortion explicitly and using thd.

t = 0:0.001:1-0.001;
x = 2*cos(2*pi*100*t)+0.01*cos(2*pi*200*t)+0.005*cos(2*pi*300*t);
tharmdist = 10*log10((0.01^2+0.005^2)/2^2)

tharmdist = 
-45.0515

r = thd(x)

r = 
-45.0515

Specify Number of Harmonics

Create a signal sampled at 1 kHz. The signal consists of a 100 Hz fundamental with amplitude 2 and three harmonics at 200, 300, and 400 Hz with amplitudes 0.01, 0.005, and 0.0025.

Set the number of harmonics to 3. This includes the fundamental. Accordingly, the power at 100, 200, and 300 Hz is used in the THD calculation.

t = 0:0.001:1-0.001;
x = 2*cos(2*pi*100*t)+0.01*cos(2*pi*200*t)+ ...
    0.005*cos(2*pi*300*t)+0.0025*sin(2*pi*400*t);
r = thd(x,1000,3)

r = 
-45.0515

Specifying the number of harmonics equal to 3 ignores the power at 400 Hz in the THD calculation.

Specify Number of Harmonics (PSD Input)

Create a signal sampled at 1 kHz. The signal consists of a 100 Hz fundamental with amplitude 2 and three harmonics at 200, 300, and 400 Hz with amplitudes 0.01, 0.005, and 0.0025, respectively.

t = 0:0.001:1-0.001;
fs = 1000;
x = 2*cos(2*pi*100*t) + 0.01*cos(2*pi*200*t)+ ...
    0.005*cos(2*pi*300*t) + 0.0025*sin(2*pi*400*t);

Obtain the periodogram PSD estimate of the signal and use the PSD estimate as the input to compute the THD. Set the number of harmonics to 3. This includes the fundamental frequency component. Accordingly, the powers at 100, 200, and 300 Hz are used in the THD calculation.

[pxx,f] = periodogram(x,rectwin(length(x)),length(x),fs);
r = thd(pxx,f,3,"psd")

r = 
-45.0515

THD from Power Spectrum

Determine the THD by inputting the power spectrum obtained with a Hamming window and the resolution bandwidth of the window.

Create a signal sampled at 10 kHz. The signal consists of a 100 Hz fundamental with amplitude 2 and three odd-numbered harmonics at 300, 500, and 700 Hz with amplitudes 0.01, 0.005, and 0.0025, respectively.

fs = 10000;
t = 0:1/fs:1-1/fs;
x = 2*cos(2*pi*100*t) + 0.01*cos(2*pi*300*t)+ ...
    0.005*cos(2*pi*500*t) + 0.0025*sin(2*pi*700*t);

Compute the THD from the power spectrum of the signal. Specify the number of harmonics to 7.

sigLength = length(x);
win = hamming(sigLength);
[sxx,f] = periodogram(x,win,sigLength,fs,"power");
rbw = enbw(win,fs);

r = thd(sxx,f,rbw,7,"power")

r = 
-44.8396

THD with and Without Aliased Harmonics

Generate a signal that resembles the output of a weakly nonlinear amplifier with a 2.1 kHz tone as input. The signal is sampled for 1 second at 10 kHz. Compute and plot the power spectrum of the signal. Use a Kaiser window with β = 38 for the computation.

Fs = 10000;
f = 2100;

t = 0:1/Fs:1; 
x = tanh(sin(2*pi*f*t)+0.1) + 0.001*randn(1,length(t));

periodogram(x,kaiser(length(x),38),[],Fs,'power')

Harmonics stick out from the noise at frequencies of 4.2 kHz, 6.3 kHz, 8.4 kHz, 10.5 kHz, 12.6 kHz, and 14.7 kHz. All frequencies except for the first one are greater than the Nyquist frequency. The harmonics are aliased respectively into 3.7 kHz, 1.6 kHz, 0.5 kHz, 2.6 kHz, and 4.7 kHz.

Compute the total harmonic distortion of the signal. By default, thd treats the aliased harmonics as part of the noise.

thd(x,Fs,7);

Repeat the computation, but now treat the aliased harmonics as part of the signal.

thd(x,Fs,7,'aliased');

Harmonic Powers and Corresponding Frequencies

Create a signal sampled at 10 kHz. The signal consists of a 100 Hz fundamental with amplitude 2 and three odd-numbered harmonics at 300, 500, and 700 Hz with amplitudes 0.01, 0.005, and 0.0025. Specify the number of harmonics to 7. Determine the THD, the power at the harmonics, and the corresponding frequencies.

fs = 10000;
t = 0:1/fs:1-1/fs;
x = 2*cos(2*pi*100*t)+0.01*cos(2*pi*300*t)+ ...
    0.005*cos(2*pi*500*t)+0.0025*sin(2*pi*700*t);
[r,harmpow,harmfreq] = thd(x,10000,7);
[harmfreq harmpow]

ans = 7×2

  100.0000    3.0103
  201.0000 -321.0983
  300.0000  -43.0103
  399.0000 -281.9259
  500.0000  -49.0309
  599.0000 -282.1066
  700.0000  -55.0515

The powers at the even-numbered harmonics are on the order of $$-300$$ dB, which corresponds to an amplitude of $$10^{-15}$$.

THD of an Amplified Signal

Generate a sinusoid of frequency 2.5 kHz sampled at 50 kHz. Add Gaussian white noise with standard deviation 0.00005 to the signal. Pass the result through a weakly nonlinear amplifier. Plot the THD.

fs = 5e4;
f0 = 2.5e3;
N = 1024;
t = (0:N-1)/fs;

ct = cos(2*pi*f0*t);
cd = ct + 0.00005*randn(size(ct));

amp = [1e-5 5e-6 -1e-3 6e-5 1 25e-3];
sgn = polyval(amp,cd);
thd(sgn,fs);

The plot shows the spectrum used to compute the ratio and the region treated as noise. The DC level is excluded from the computation. The fundamental and harmonics are labeled.

Input Arguments

x — Real-valued sinusoidal input signal

vector

Real-valued sinusoidal input signal, specified as a row or column vector.

Example: cos(pi/4*(0:159))+cos(pi/2*(0:159))

Data Types: single | double

fs — Sample rate

positive scalar

Sample rate, specified as a positive scalar. The sample rate is the number of samples per unit time. If the unit of time is seconds, then the sample rate has units of Hz.

n — Number of harmonics

positive integer

Number of harmonics, specified as a positive integer.

pxx — One-sided PSD estimate

vector

One-sided PSD estimate, specified as a real-valued, nonnegative column vector.

The power spectral density must be expressed in linear units, not decibels. Use db2pow to convert decibel values to power values.

Example: [pxx,f] = periodogram(cos(pi./[4;2]*(0:159))'+randn(160,2)) specifies the periodogram PSD estimate of a noisy two-channel sinusoid sampled at 2π Hz and the frequencies at which it is computed.

Data Types: single | double

f — Cyclical frequencies

vector

Cyclical frequencies corresponding to the one-sided PSD estimate, pxx, specified as a row or column vector. The first element of f must be 0.

Data Types: double | single

sxx — Power spectrum

nonnegative real-valued row or column vector

Power spectrum, specified as a real-valued nonnegative row or column vector.

The power spectrum must be expressed in linear units, not decibels. Use db2pow to convert decibel values to power values.

Example: [sxx,w] = periodogram(cos(pi./[4;2]*(0:159))'+randn(160,2),'power') specifies the periodogram power spectrum estimate of a two-channel sinusoid embedded in white Gaussian noise and the normalized frequencies at which it is computed.

rbw — Resolution bandwidth

positive scalar

Resolution bandwidth, specified as a positive scalar. The resolution bandwidth is the product of the frequency resolution of the discrete Fourier transform and the equivalent noise bandwidth of the window.

Output Arguments

r — Total harmonic distortion in dBc

real-valued scalar

Total harmonic distortion in dBc, returned as a real-valued scalar.

harmpow — Power of the harmonics

real-valued scalar or vector

Power of the harmonics, returned as a real-valued scalar or vector expressed in dB. Whether harmpow is a scalar or a vector depends on the number of harmonics you specify as the input argument n.

harmfreq — Frequencies of the harmonics

nonnegative scalar or vector

Frequencies of the harmonics, returned as a nonnegative scalar or vector. Whether harmfreq is a scalar or a vector depends on the number of harmonics you specify as the input argument n.

More About

Distortion Measurement Functions

The functions thd, sfdr, sinad, and snr measure the response of a weakly nonlinear system stimulated by a sinusoid.

When given time-domain input, thd performs a periodogram using a Kaiser window with large sidelobe attenuation. To find the fundamental frequency, the algorithm searches the periodogram for the largest nonzero spectral component. It then computes the central moment of all adjacent bins that decrease monotonically away from the maximum. To be detectable, the fundamental should be at least in the second frequency bin. Higher harmonics are at integer multiples of the fundamental frequency. If a harmonic lies within the monotonically decreasing region in the neighborhood of another, its power is considered to belong to the larger harmonic. This larger harmonic may or may not be the fundamental.

thd fails if the fundamental is not the highest spectral component in the signal.

Ensure that the frequency components are far enough apart to accommodate for the sidelobe width of the Kaiser window. If this is not feasible, you can use the "power" flag and compute a periodogram with a different window.

Extended Capabilities

C/C++ Code Generation

Generate C and C++ code using MATLAB® Coder™.

Usage notes and limitations:

• If supplied, input arguments "power", "psd", "aliased", and "omitaliases" must be compile-time constants.

GPU Arrays

Accelerate code by running on a graphics processing unit (GPU) using Parallel Computing Toolbox™.

This function fully supports GPU arrays. For more information, see Run MATLAB Functions on a GPU (Parallel Computing Toolbox).

Version History

Introduced in R2013b

R2023a: Use gpuArray objects

The thd function supports gpuArray objects. You must have Parallel Computing Toolbox™ to use this functionality.