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 sampling 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(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.

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.

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)
r = thd(x)

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.

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.

Obtain the periodogram PSD estimate of the signal and use the
PSD estimate as the input to thd. 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;
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);
[pxx,f] = periodogram(x,rectwin(length(x)),length(x),fs);
r = thd(pxx,f,3,'psd');

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. Specify
the number of harmonics to 7. Determine the THD.

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);
[sxx,f] = periodogram(x,hamming(length(x)),length(x),fs,'power');
rbw = enbw(hamming(length(x)),fs);
r = thd(sxx,f,rbw,7,'power');

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];

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

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

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.

Sampling frequency, specified as a positive scalar. The sampling
frequency is the number of samples per unit time. If the unit of time
is seconds, the sampling frequency has units of hertz.

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.

Power of the harmonics specified as a nonnegative scalar or
vector. Whether harmpow is a scalar or a vector
depends on the number of harmonics you specify as the input argument n.

Frequencies of the harmonics specified 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.

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.