Detecting frequencies in a signal

Contents

You need to download an m-file

You first have to download the m-file minabs.m and put it in the same folder as your other m-files.

Sample signal

Our sample signal contains two frequencies $a_1$ and $a_2$, a nonperiodic part, and noise

rng default;    % reset random number generator so that we get the same results each time

h = .01;
x = -5:.01:5;
N = length(x);

a(1) = sqrt(2); a(2) = exp(1);      % the two frequencies 1.4142 and 2.7183 we want to detect
f = exp(-2*abs(x)) + cos(2*pi*a(1)*x-3) + cos(2*pi*a(2)*x-2);

plot(x,f)
title('signal without noise')

Add random noise

We add some normally distributed white noise.

How accurately can we approximate $a_1$ and $a_2$ from the noisy signal?

f = f + 0.5*randn(1,N);
plot(x,f)
title('signal with noise')

Method 1: Use autocorrelation

The first relative maximum corresponds to the period $T_1 = 1/\xi_1$, the second relative maximum to the period $T_2 = 1/\xi_2$. Note that we will also relative maxima near integer multiples of $T_1$ and $T_2$.

We find the first three relative maxima of the autocorrelation vector using findpeaks. Note that the first and third value seem to correspond to "true relative maxima" whereas the second value seems to be spurious due to the noise in the data.

Therefore we detect the frequencies 2.7027 and 1.3514 (the true values are 2.7183 and 1.4142).

N2 = 2*N;
% F = [f,zeros(1,N)];      % pad with zeros
F = f;                     % without padding (also works)
Fh = four(F);
Gh = abs(Fh).^2;
G = ifour(Gh);             % autocorrelation
xv = (0:N-1)*h;            % entries of G correspond to autocorrelation for 0,h,2h,...

[pks,ind] = findpeaks(G);  % ind are indices of local maxima of vector
lm = ind(1:3);             % first 3 local maxima
periods = xv(lm)           % corresponding x-values give periods
freq = 1./xv(lm)           % corresponding frequencies

error1 = freq(3) - a(1)    % error for freq. a1
error2 = freq(1) - a(2)    % error for freq. a2

figure(1)
plot(xv(1:500),G(1:500),'.-'); grid on
title('autocorrelation')
figure(2)
plot(xv(1:100),G(1:100),'.-',xv(lm),G(lm),'ro'); grid on
title('autocorrelation with 3 relative maxima')
periods =
                      0.37                      0.39                      0.74
freq =
           2.7027027027027          2.56410256410256          1.35135135135135
error1 =
       -0.0628622110217438
error2 =
       -0.0155791257563429

Method 2: Use Discrete Fourier Transform to detect frequencies

We take the Discrete Fourier Transform (DFT) of the data vector (without any padding).

We observe sharp peaks at $|\hat f_{k_1}|$, $|\hat f_{k_2}|$, $|\hat f_{N-k_1}|$, $|\hat f_{N-k_2}|$.

In Matlab we can locate the "peaks" in a vector v using findpeaks:

[pks,ind] = findpeaks(v,'MinPeakProminence',.1,'MinPeakHeight',.2)

finds the values and indices of relative maxima which stand out from their "reference level" by at least .1 and have a height of at least .2. We have to pick the values for MinPeakProminence and MinPeakHeigth carefully so that we capture only the "true peaks" and not the relative maxima from the "noise".

For the DFT the data values correspond to the discrete times 0,1/N,...,(N-1)/N. The k-values correspond to this spacing.

However, our original data $f_j$ correspond to a spacing of h rather than 1/h. Therefore we have to divide the obtained k-values by N*h to obtain the frequencies.

Here we obtain the frequencies 1.3986 and 2.6973 (the true values are 1.4142 and 2.7183).

fh = four(f);              % Discrete Fourier Transform of signal
v = abs(fh(1:100));
kv = 0:99;
[pks,ind] = findpeaks(v,'MinPeakProminence',.1);  % find relative maxima
                     % which stand out from their "reference level" by at least .1
                     % (use "doc findpeaks" to see the definition of "reference level")
                     % You can also specify 'MinPeakHeight'

kp = kv(ind)         % k-values are frequencies for spacing 1/N
freq = kp/(N*h)      % convert to frequencies for spacing h
errors = freq - a    % errors for frequencies a(1), a(2)

plot(kv,v,'.-',kv(ind),v(ind),'ro')  % plot |fhat_k|, mark detected peaks with red circles
xlabel('k'); title('|fhat_k|')
kp =
    14    27
freq =
           1.3986013986014           2.6973026973027
errors =
       -0.0156121637716966       -0.0209791311563481

More accurate approximations for 1st and 2nd frequency

The first peak for $|\hat f_k|$ occurs at $k=14$, which gave us the frequency estimate 1.399 with an error of .015.

If we connect the points $z_k:=1/\hat f_k$ for $k=13,14,15$ in the complex plane by two line segments, we can find the closest point to the origin. On this polygonal line the closest point to the origin corresponds to kpeak = 14.14, corresponding the frequency xipeak=kpeak/(N*h)=1.4147. The error for the 1st frequency is now only 4.6e-4 (before it was .015).

We proceed analogously for the 2nd frequency: The second peak for $|\hat f_k|$ occurs at $k=27$, which gave us the frequency estimate 2.697 with an error of .02. Using minabs we obtain the improved frequency estimate xipeak=2.719. The error for the 1st frequency is now only 1.1e-3 (before it was .02).

for j=1:2
  K = kp(j)                   % location of j-th peak: K=14 for j=1, K=27 for j=1
  absvals = abs(fh(K:K+2))    % [ |fhat_{K-1}|, |fhat_K|, |fhat_{K+1}| ]
  t = minabs(1./fh(K:K+2))
  kpeak = K + t               % best estimate for location of peak
  xipeak = kpeak/(N*h)        % convert to frequencies for original data at x_j = j*h
  fprintf('error = %e\n',xipeak-a(j))
end
K =
    14
absvals =
        0.0666782093240122         0.476717377011674        0.0916745410923054
t =
         0.160852557431261
kpeak =
          14.1608525574313
xipeak =
          1.41467058515797
error = 4.570228e-04
K =
    27
absvals =
        0.0773635422576332         0.439933703779493         0.125015261384777
t =
         0.221285136319633
kpeak =
          27.2212851363196
xipeak =
          2.71940910452744
error = 1.127276e-03