Estimate frequency and peak value of a signal's fundamental

Question

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.

Answer 1

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.

Answer 2

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

Answer 3

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

Answer 4

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?