Skip to content

Instantly share code, notes, and snippets.

@MartinBloedorn
Created August 7, 2018 21:07
Show Gist options
  • Save MartinBloedorn/4935cd3f923e4282212ff6ee9b19020f to your computer and use it in GitHub Desktop.
Save MartinBloedorn/4935cd3f923e4282212ff6ee9b19020f to your computer and use it in GitHub Desktop.
function [thd, p0, l0] = thdf( P1, f, varargin )
%thdf Computes the total harmonic distortion to the fundamental.
% Fundamental is the first harmonic peak that's larger than 2% than
% the largest harmonic in the supplied frequency spectrum.
%
% thd = thdf( P1, f )
%
% P1 Single-sided amplitude spectrum of the signal.
% f Frequency span of P1.
% thd Total harmonic distortion.
%
% Options:
% 'NFundamentals' [1] Amount of fundamentals assumed in computation.
% 'FThres' [Inf] Upper limit for accounted frequencies.
p = inputParser;
p.addOptional('NFundamentals', 1, @isnumeric);
p.addOptional('FThres', Inf, @isnumeric);
p.parse(varargin{:});
% Amount of fundamental harmonics
nf = p.Results.NFundamentals;
% Frequency threshold
ft = p.Results.FThres;
% Fundamental amplitude threshold: >2% than the largest harmonic in P1
athres = 0.02*max(P1);
% Find all harmonics above threshold
[p, l] = findpeaks(P1, f);
% Filter out irrelevant harmonics
p0 = p(p > athres & l < ft);
l0 = l(p > athres & l < ft);
% Fundamentals
p1 = p0(1:nf);
% Compute distortion
thd = sqrt(sum(p0((nf+1):end).^2))/sqrt(sum(p1.^2));
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment