# Estimate frequency and peak value of a signal's fundamental

How can I find out frequency and peak value of a signal's fundamental? Please do consider the following:

• The peak value of some harmonics might be higher than that one of the fundamental
• The frequency of the harmonic is not an integer (e.g. 3.45 Hz)
• I have information about where to search for the frequency (3.45 +- 0.5 Hz).
• I need to estimate frequency and peak value as accurate as possible

Any response or sample code (MATLAB) is highly appreciated.

• The answer will depend on whether you are looking for the frequency of only the fundamental spectrum peak, or the f0 of the fundamental pitch or periodicity, which can be different depending on whether the harmonics are exactly harmonic or slightly inharmonic as in many real world signals. Nov 7, 2013 at 20:04

There are a few Matlab implementations of various frequency estimators here. However, they will probably be "tricked" by the higher harmonics.

That you know a range, perhaps the best technique is one of Eric Jacobsen's estimators. Have a look at his page.

In essence, Eric's estimators use the FFT bin at $X(k)$, and on either side of ($X(k+1), X(k-1)$) the peak. If you choose your FFT length to match your known range (so that $k \rightarrow 3.45$ and $k+1 \rightarrow 3.50$ and $k-1 \rightarrow 3.4$, then you should be good to go.

Harmonic product Spectrum (HPS) seems be useful for you, a feature of the HPS algorithm is multiply a number of harmonics, this give you a peak of your fundamental frequency, I believe that you will need a big FFT points like 12288 for signals sampled at 44100hz, for it the first bin are in 3.5889hz. HPS source code example here: https://stackoverflow.com/questions/19765486/matlab-code-for-harmonic-product-spectrum

• Thank you for your answer. Unfortunately, using the code example without the correctFactor, the result was not accurate enough for my requirements. I need to estimate the frequency with a tolerance of only +- 0.01 Hz. Nov 6, 2013 at 20:21

If the harmonics are exactly harmonic (and not slightly inharmonic as produced by many physical systems), and you have a good frequency estimate of a strong harmonic, such that that harmonic frequency estimate fits within a unique integer multiple of the search range of the fundamental frequency, then you can just divide by that integer multiple to get a frequency estimate of the fundamental.

Once you have a frequency estimate, you can then compute a magnitude estimate using 1 bin of a DFT or a Goertzel filter with a length that is exactly an integer multiple of the period of that frequency estimate. (... or 3 very close lengths followed by parabolic peak interpolation.)

• Thanks a lot for this comment. I have adapted the code below which leads to good results concerning frequency estimation. I do not fully understand what is the most accurate way to estimate the fundamentals peak value. Can I just obtain it at the bin corresponding to the estimated frequency or would your recommend another method which leads to better results? Nov 8, 2013 at 10:33
• If the estimated frequency is between bin centers, then the nearest bin value may not be accurate. However a different length of FFT or DFT may put the estimate closer to a bin center. Nov 8, 2013 at 18:30
• Thank you for your response. This is enough for this question. I will create a new one dealing with peak estimation in more detail. Nov 9, 2013 at 21:19

Thanks a lot, that looks promising.

I have used Erics code (Macleod's algorithm) to implement the following snipped:

% Amplitude of signal.
ACW=1;

% Specify signal.
Fs = 8000;                   % samples per second
dt = 1/Fs;                   % seconds per sample
nHarmonic = 7;               % order of harmonic
Fc = 21.5;                   % Hz
StopTime = 0.25;             % seconds
t = (0:dt:StopTime-dt)';     % seconds

% Specify evaluation tolerance
Tol = 0.01;

% Effective length of test vector.
tstlen=length(t);

% Generate signal
cw = ACW * sin(2*pi*(1:tstlen)*Fc/tstlen) + 1.5 * ACW * sin(2*pi*(1:tstlen)*(Fc/tstlen) * 3) ...
+ 2.5 * ACW * sin(2*pi*(1:tstlen)*(Fc/tstlen) * 5)+ 0.3 * ACW * sin(2*pi*(1:tstlen)*(Fc/tstlen) * 7);

% Number of trials in test.
N=10000;

% Desired bin number of tone relative to long test length.
bin=Fc;

% Allocate vector for results.
macldest=zeros(size(1:N));

% Calculate desired target result for comparison.
target=bin;

fprintf('Peak is at %f.\n',bin);

for I = 1:N,        % Run trials.
% Calculate signal power.
sigp=sum(abs(cw(1:tstlen).^2));

% DFT.
dftshrt(1:tstlen)=fft(cw(1:tstlen));

% DFT magnitude.
magshrt(1:tstlen)=abs(dftshrt);

% Find raw peak magnitude and location.
[rawmag,rawind]=max(magshrt);

% Isolated 3 samples around peak.
pk3vect(1:3)=dftshrt(rawind-1:rawind+1);

%Do Macleod's estimation.
ref = pk3vect(2);                               % Isolate phase reference.
R = real(pk3vect.*conj(ref));                   % Generate phase corrected coefficient vector.
gamma = (R(1)-R(3))/((2*R(2))+R(1)+R(3));       % Calculate offset.
delta = (sqrt(1 + 8*gamma*gamma)-1)/(4*gamma);  % Final estimate.
macldest=rawind-1+delta;
end

% Calculate average result.
macldr=mean(macldest)

% Evaluate estimated frequency
Factor = round(macldr/Fc)

if (macldr > (Fc-Tol) & macldr<(Fc+Tol) )
result=macldr
elseif ((macldr/Factor)>(Fc-Tol) & (macldr/Factor)<(Fc+Tol))
result=macldr/Factor
end


But what is the most accurate method now to estimate the fundamentals peak value?