# thd

Total harmonic distortion

## Description

example

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.

example

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.

example

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

example

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.

example

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.

example

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

example

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

collapse all

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

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')
r = -45.0515

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')
r = -44.8396

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');

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 $1{0}^{-15}$.

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

collapse all

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

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.

Number of harmonics, specified as a positive integer.

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

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

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.

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

collapse all

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

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.

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.

collapse all

### 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.

## Version History

Introduced in R2013b